[3294] | 1 | !> @file ocean_mod.f90 |
---|
[4797] | 2 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 3 | ! This file is part of the PALM model system. |
---|
| 4 | ! |
---|
[4797] | 5 | ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General |
---|
| 6 | ! Public License as published by the Free Software Foundation, either version 3 of the License, or |
---|
| 7 | ! (at your option) any later version. |
---|
[3294] | 8 | ! |
---|
[4797] | 9 | ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the |
---|
| 10 | ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General |
---|
| 11 | ! Public License for more details. |
---|
[3294] | 12 | ! |
---|
[4797] | 13 | ! You should have received a copy of the GNU General Public License along with PALM. If not, see |
---|
| 14 | ! <http://www.gnu.org/licenses/>. |
---|
[3294] | 15 | ! |
---|
[4481] | 16 | ! Copyright 2017-2020 Leibniz Universitaet Hannover |
---|
[4797] | 17 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 18 | ! |
---|
| 19 | ! Current revisions: |
---|
| 20 | ! ----------------- |
---|
[4671] | 21 | ! |
---|
| 22 | ! |
---|
[3294] | 23 | ! Former revisions: |
---|
| 24 | ! ----------------- |
---|
| 25 | ! $Id: ocean_mod.f90 4797 2020-11-26 16:02:39Z eckhard $ |
---|
[4797] | 26 | ! file re-formatted to follow the PALM coding standard |
---|
| 27 | ! |
---|
| 28 | ! 4768 2020-11-02 19:11:23Z suehring |
---|
[4768] | 29 | ! Enable 3D data output also with 64-bit precision |
---|
| 30 | ! |
---|
| 31 | ! 4731 2020-10-07 13:25:11Z schwenkel |
---|
[4731] | 32 | ! Move exchange_horiz from time_integration to modules |
---|
| 33 | ! |
---|
| 34 | ! 4671 2020-09-09 20:27:58Z pavelkrc |
---|
[4671] | 35 | ! Implementation of downward facing USM and LSM surfaces |
---|
[4797] | 36 | ! |
---|
[4671] | 37 | ! 4535 2020-05-15 12:07:23Z raasch |
---|
[4535] | 38 | ! bugfix for restart data format query |
---|
[4797] | 39 | ! |
---|
[4535] | 40 | ! 4517 2020-05-03 14:29:30Z raasch |
---|
[4517] | 41 | ! added restart with MPI-IO for reading local arrays |
---|
[4797] | 42 | ! |
---|
[4517] | 43 | ! 4495 2020-04-13 20:11:20Z raasch |
---|
[4495] | 44 | ! restart data handling with MPI-IO added |
---|
[4797] | 45 | ! |
---|
[4495] | 46 | ! 4481 2020-03-31 18:55:54Z maronga |
---|
[4370] | 47 | ! vector directives added to force vectorization on Intel19 compiler |
---|
[4797] | 48 | ! |
---|
[4370] | 49 | ! 4346 2019-12-18 11:55:56Z motisi |
---|
[4797] | 50 | ! Introduction of wall_flags_total_0, which currently sets bits based on static topography |
---|
| 51 | ! information used in wall_flags_static_0 |
---|
| 52 | ! |
---|
[4346] | 53 | ! 4329 2019-12-10 15:46:36Z motisi |
---|
[4329] | 54 | ! Renamed wall_flags_0 to wall_flags_static_0 |
---|
[4797] | 55 | ! |
---|
[4329] | 56 | ! 4272 2019-10-23 15:18:57Z schwenkel |
---|
[4797] | 57 | ! Further modularization of boundary conditions: moved boundary conditions to respective modules |
---|
[4272] | 58 | ! |
---|
| 59 | ! 4196 2019-08-29 11:02:06Z gronemeier |
---|
[4196] | 60 | ! Consider rotation of model domain for calculating the Stokes drift |
---|
[4797] | 61 | ! |
---|
[4196] | 62 | ! 4182 2019-08-22 15:20:23Z scharf |
---|
[4182] | 63 | ! Corrected "Former revisions" section |
---|
[4797] | 64 | ! |
---|
[4182] | 65 | ! 4110 2019-07-22 17:05:21Z suehring |
---|
[4797] | 66 | ! Pass integer flag array as well as boundary flags to WS scalar advection routine |
---|
| 67 | ! |
---|
[4110] | 68 | ! 4109 2019-07-22 17:00:34Z suehring |
---|
[3873] | 69 | ! implemented ocean_actions |
---|
[4797] | 70 | ! |
---|
[3873] | 71 | ! 3767 2019-02-27 08:18:02Z raasch |
---|
[4797] | 72 | ! unused variable for file index and tmp_2d removed from rrd-subroutine parameter list |
---|
| 73 | ! |
---|
[3767] | 74 | ! 3719 2019-02-06 13:10:18Z kanani |
---|
[4797] | 75 | ! Changed log_point to log_point_s, otherwise this overlaps with 'all progn.equations' |
---|
| 76 | ! cpu measurement. |
---|
| 77 | ! |
---|
[3719] | 78 | ! 3684 2019-01-20 20:20:58Z knoop |
---|
[3636] | 79 | ! nopointer option removed |
---|
[4797] | 80 | ! |
---|
[4182] | 81 | ! 3294 2018-10-01 02:37:10Z raasch |
---|
| 82 | ! initial revision |
---|
[3294] | 83 | ! |
---|
[4182] | 84 | ! |
---|
| 85 | ! Authors: |
---|
| 86 | ! -------- |
---|
| 87 | ! @author Siegfried Raasch |
---|
| 88 | ! |
---|
[3294] | 89 | ! Description: |
---|
| 90 | ! ------------ |
---|
[4797] | 91 | !> This module contains relevant code for PALM's ocean mode, e.g. the prognostic equation for |
---|
| 92 | !> salinity, the equation of state for seawater, the Craik Leibovich force (Stokes force), and wave |
---|
| 93 | !> breaking effects |
---|
| 94 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 95 | MODULE ocean_mod |
---|
| 96 | |
---|
| 97 | |
---|
[4797] | 98 | USE arrays_3d, & |
---|
| 99 | ONLY: diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, prho, prho_1, rho_1, rho_ocean, sa, & |
---|
| 100 | sa_1, sa_2, sa_3, sa_init, sa_p, tsa_m |
---|
[3294] | 101 | |
---|
[4797] | 102 | USE control_parameters, & |
---|
| 103 | ONLY: atmos_ocean_sign, bottom_salinityflux, constant_top_salinityflux, loop_optimization,& |
---|
| 104 | ocean_mode, restart_data_format_output, top_salinityflux, wall_salinityflux, & |
---|
| 105 | ws_scheme_sca |
---|
| 106 | |
---|
[3294] | 107 | USE kinds |
---|
| 108 | |
---|
[4797] | 109 | USE pegrid, & |
---|
[3873] | 110 | ONLY: threads_per_task |
---|
[3294] | 111 | |
---|
[4797] | 112 | USE statistics, & |
---|
[3873] | 113 | ONLY: sums_wssas_ws_l |
---|
| 114 | |
---|
[4797] | 115 | USE indices, & |
---|
| 116 | ONLY: advc_flags_s, nbgp, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0 |
---|
[3873] | 117 | |
---|
[4495] | 118 | USE restart_data_mpi_io_mod, & |
---|
[4517] | 119 | ONLY: rd_mpi_io_check_array, rrd_mpi_io, rrd_mpi_io_global_array, wrd_mpi_io, & |
---|
| 120 | wrd_mpi_io_global_array |
---|
[4495] | 121 | |
---|
[4797] | 122 | USE surface_mod, & |
---|
| 123 | ONLY: bc_h, surf_def_v, surf_def_h, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v |
---|
[4272] | 124 | |
---|
[4731] | 125 | USE exchange_horiz_mod, & |
---|
| 126 | ONLY: exchange_horiz |
---|
| 127 | |
---|
| 128 | |
---|
[3294] | 129 | IMPLICIT NONE |
---|
| 130 | |
---|
| 131 | CHARACTER (LEN=20) :: bc_sa_t = 'neumann' !< namelist parameter |
---|
| 132 | |
---|
[4797] | 133 | INTEGER(iwp) :: ibc_sa_t !< integer flag for bc_sa_t |
---|
| 134 | INTEGER(iwp) :: iran_ocean = -1234567 !< random number used for wave breaking effects |
---|
[3294] | 135 | INTEGER(iwp) :: sa_vertical_gradient_level_ind(10) = -9999 !< grid index values of sa_vertical_gradient_level(s) |
---|
| 136 | |
---|
[3303] | 137 | LOGICAL :: salinity = .TRUE. !< switch for using salinity |
---|
[3302] | 138 | LOGICAL :: stokes_force = .FALSE. !< switch to switch on the Stokes force |
---|
[4797] | 139 | LOGICAL :: surface_cooling_switched_off = .FALSE. !< variable to check if surface heat flux has been switched off |
---|
[3302] | 140 | LOGICAL :: wave_breaking = .FALSE. !< switch to switch on wave breaking effects |
---|
[3294] | 141 | |
---|
[3302] | 142 | REAL(wp) :: alpha_wave_breaking = 3.0_wp !< coefficient for wave breaking generated turbulence from Noh et al. (2004), JPO |
---|
[3294] | 143 | REAL(wp) :: prho_reference !< reference state of potential density at ocean surface |
---|
| 144 | REAL(wp) :: sa_surface = 35.0_wp !< surface salinity, namelist parameter |
---|
| 145 | REAL(wp) :: sa_vertical_gradient(10) = 0.0_wp !< namelist parameter |
---|
| 146 | REAL(wp) :: sa_vertical_gradient_level(10) = -999999.9_wp !< namelist parameter |
---|
[3302] | 147 | REAL(wp) :: stokes_waveheight = 0.0_wp !< wave height assumed for Stokes drift velocity |
---|
| 148 | REAL(wp) :: stokes_wavelength = 0.0_wp !< wavelength assumed for Stokes drift velocity |
---|
[3381] | 149 | REAL(wp) :: surface_cooling_spinup_time = 999999.9_wp !< time after which surface heat flux is switched off |
---|
[3302] | 150 | REAL(wp) :: timescale_wave_breaking !< time scale of random forcing |
---|
| 151 | REAL(wp) :: u_star_wave_breaking !< to store the absolute value of friction velocity at the ocean surface |
---|
[3294] | 152 | |
---|
[4797] | 153 | REAL(wp), DIMENSION(12), PARAMETER :: nom = & |
---|
| 154 | (/ 9.9984085444849347D2, 7.3471625860981584D0, & |
---|
| 155 | -5.3211231792841769D-2, 3.6492439109814549D-4, & |
---|
| 156 | 2.5880571023991390D0, -6.7168282786692354D-3, & |
---|
| 157 | 1.9203202055760151D-3, 1.1798263740430364D-2, & |
---|
| 158 | 9.8920219266399117D-8, 4.6996642771754730D-6, & |
---|
| 159 | -2.5862187075154352D-8, -3.2921414007960662D-12 /) |
---|
| 160 | !< constants used in equation of state for seawater |
---|
[3294] | 161 | |
---|
[4797] | 162 | REAL(wp), DIMENSION(13), PARAMETER :: den = & |
---|
| 163 | (/ 1.0D0, 7.2815210113327091D-3, & |
---|
| 164 | -4.4787265461983921D-5, 3.3851002965802430D-7, & |
---|
| 165 | 1.3651202389758572D-10, 1.7632126669040377D-3, & |
---|
| 166 | -8.8066583251206474D-6, -1.8832689434804897D-10, & |
---|
| 167 | 5.7463776745432097D-6, 1.4716275472242334D-9, & |
---|
| 168 | 6.7103246285651894D-6, -2.4461698007024582D-17, & |
---|
| 169 | -9.1534417604289062D-18 /) |
---|
| 170 | !< constants used in equation of state for seawater |
---|
[3294] | 171 | |
---|
| 172 | SAVE |
---|
| 173 | |
---|
[4797] | 174 | PUBLIC :: bc_sa_t, ibc_sa_t, prho_reference, sa_surface, sa_vertical_gradient, & |
---|
| 175 | sa_vertical_gradient_level, sa_vertical_gradient_level_ind, stokes_force, & |
---|
| 176 | wave_breaking |
---|
[3294] | 177 | |
---|
| 178 | |
---|
| 179 | INTERFACE eqn_state_seawater |
---|
| 180 | MODULE PROCEDURE eqn_state_seawater |
---|
| 181 | MODULE PROCEDURE eqn_state_seawater_ij |
---|
| 182 | END INTERFACE eqn_state_seawater |
---|
| 183 | |
---|
| 184 | INTERFACE eqn_state_seawater_func |
---|
| 185 | MODULE PROCEDURE eqn_state_seawater_func |
---|
| 186 | END INTERFACE eqn_state_seawater_func |
---|
| 187 | |
---|
| 188 | INTERFACE ocean_check_parameters |
---|
| 189 | MODULE PROCEDURE ocean_check_parameters |
---|
| 190 | END INTERFACE ocean_check_parameters |
---|
| 191 | |
---|
| 192 | INTERFACE ocean_check_data_output |
---|
| 193 | MODULE PROCEDURE ocean_check_data_output |
---|
| 194 | END INTERFACE ocean_check_data_output |
---|
| 195 | |
---|
| 196 | INTERFACE ocean_check_data_output_pr |
---|
| 197 | MODULE PROCEDURE ocean_check_data_output_pr |
---|
| 198 | END INTERFACE ocean_check_data_output_pr |
---|
| 199 | |
---|
| 200 | INTERFACE ocean_define_netcdf_grid |
---|
| 201 | MODULE PROCEDURE ocean_define_netcdf_grid |
---|
| 202 | END INTERFACE ocean_define_netcdf_grid |
---|
| 203 | |
---|
| 204 | INTERFACE ocean_data_output_2d |
---|
| 205 | MODULE PROCEDURE ocean_data_output_2d |
---|
| 206 | END INTERFACE ocean_data_output_2d |
---|
| 207 | |
---|
| 208 | INTERFACE ocean_data_output_3d |
---|
| 209 | MODULE PROCEDURE ocean_data_output_3d |
---|
| 210 | END INTERFACE ocean_data_output_3d |
---|
| 211 | |
---|
[4731] | 212 | INTERFACE ocean_exchange_horiz |
---|
| 213 | MODULE PROCEDURE ocean_exchange_horiz |
---|
| 214 | END INTERFACE ocean_exchange_horiz |
---|
| 215 | |
---|
[3302] | 216 | INTERFACE ocean_header |
---|
| 217 | MODULE PROCEDURE ocean_header |
---|
| 218 | END INTERFACE ocean_header |
---|
| 219 | |
---|
[3294] | 220 | INTERFACE ocean_init |
---|
| 221 | MODULE PROCEDURE ocean_init |
---|
| 222 | END INTERFACE ocean_init |
---|
| 223 | |
---|
| 224 | INTERFACE ocean_init_arrays |
---|
| 225 | MODULE PROCEDURE ocean_init_arrays |
---|
| 226 | END INTERFACE ocean_init_arrays |
---|
| 227 | |
---|
| 228 | INTERFACE ocean_parin |
---|
| 229 | MODULE PROCEDURE ocean_parin |
---|
| 230 | END INTERFACE ocean_parin |
---|
| 231 | |
---|
[3873] | 232 | INTERFACE ocean_actions |
---|
| 233 | MODULE PROCEDURE ocean_actions |
---|
| 234 | MODULE PROCEDURE ocean_actions_ij |
---|
| 235 | END INTERFACE ocean_actions |
---|
| 236 | |
---|
[3294] | 237 | INTERFACE ocean_prognostic_equations |
---|
| 238 | MODULE PROCEDURE ocean_prognostic_equations |
---|
| 239 | MODULE PROCEDURE ocean_prognostic_equations_ij |
---|
| 240 | END INTERFACE ocean_prognostic_equations |
---|
| 241 | |
---|
[4272] | 242 | INTERFACE ocean_boundary_conditions |
---|
| 243 | MODULE PROCEDURE ocean_boundary_conditions |
---|
| 244 | END INTERFACE ocean_boundary_conditions |
---|
| 245 | |
---|
[3294] | 246 | INTERFACE ocean_swap_timelevel |
---|
| 247 | MODULE PROCEDURE ocean_swap_timelevel |
---|
| 248 | END INTERFACE ocean_swap_timelevel |
---|
| 249 | |
---|
| 250 | INTERFACE ocean_rrd_global |
---|
[4495] | 251 | MODULE PROCEDURE ocean_rrd_global_ftn |
---|
| 252 | MODULE PROCEDURE ocean_rrd_global_mpi |
---|
[3294] | 253 | END INTERFACE ocean_rrd_global |
---|
| 254 | |
---|
| 255 | INTERFACE ocean_rrd_local |
---|
[4517] | 256 | MODULE PROCEDURE ocean_rrd_local_ftn |
---|
| 257 | MODULE PROCEDURE ocean_rrd_local_mpi |
---|
[3294] | 258 | END INTERFACE ocean_rrd_local |
---|
| 259 | |
---|
| 260 | INTERFACE ocean_wrd_global |
---|
| 261 | MODULE PROCEDURE ocean_wrd_global |
---|
| 262 | END INTERFACE ocean_wrd_global |
---|
| 263 | |
---|
| 264 | INTERFACE ocean_wrd_local |
---|
| 265 | MODULE PROCEDURE ocean_wrd_local |
---|
| 266 | END INTERFACE ocean_wrd_local |
---|
| 267 | |
---|
| 268 | INTERFACE ocean_3d_data_averaging |
---|
| 269 | MODULE PROCEDURE ocean_3d_data_averaging |
---|
| 270 | END INTERFACE ocean_3d_data_averaging |
---|
| 271 | |
---|
[3302] | 272 | INTERFACE stokes_drift_terms |
---|
| 273 | MODULE PROCEDURE stokes_drift_terms |
---|
| 274 | MODULE PROCEDURE stokes_drift_terms_ij |
---|
| 275 | END INTERFACE stokes_drift_terms |
---|
[3294] | 276 | |
---|
[3302] | 277 | INTERFACE wave_breaking_term |
---|
| 278 | MODULE PROCEDURE wave_breaking_term |
---|
| 279 | MODULE PROCEDURE wave_breaking_term_ij |
---|
| 280 | END INTERFACE wave_breaking_term |
---|
| 281 | |
---|
[3294] | 282 | PRIVATE |
---|
| 283 | ! |
---|
| 284 | !-- Add INTERFACES that must be available to other modules (alphabetical order) |
---|
[4797] | 285 | PUBLIC eqn_state_seawater, ocean_actions, ocean_check_data_output, ocean_check_data_output_pr, & |
---|
| 286 | ocean_check_parameters, ocean_data_output_2d, ocean_data_output_3d, & |
---|
| 287 | ocean_exchange_horiz, ocean_define_netcdf_grid, ocean_header, ocean_init, & |
---|
| 288 | ocean_init_arrays, ocean_parin, ocean_prognostic_equations, ocean_swap_timelevel, & |
---|
| 289 | ocean_rrd_global, ocean_rrd_local, ocean_wrd_global, ocean_wrd_local, & |
---|
| 290 | ocean_3d_data_averaging, ocean_boundary_conditions, stokes_drift_terms, & |
---|
| 291 | wave_breaking_term |
---|
[3294] | 292 | |
---|
| 293 | |
---|
| 294 | CONTAINS |
---|
| 295 | |
---|
[4797] | 296 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 297 | ! Description: |
---|
| 298 | ! ------------ |
---|
[4797] | 299 | !> Equation of state for seawater as a function of potential temperature, salinity, and pressure. |
---|
[3294] | 300 | !> For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech. |
---|
| 301 | !> eqn_state_seawater calculates the potential density referred at hyp(0). |
---|
| 302 | !> eqn_state_seawater_func calculates density. |
---|
| 303 | !> TODO: so far, routine is not adjusted to use topography |
---|
[4797] | 304 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 305 | SUBROUTINE eqn_state_seawater |
---|
| 306 | |
---|
[4797] | 307 | USE arrays_3d, & |
---|
[3294] | 308 | ONLY: hyp, prho, pt_p, rho_ocean, sa_p |
---|
[4797] | 309 | USE indices, & |
---|
[3294] | 310 | ONLY: nxl, nxr, nyn, nys, nzb, nzt |
---|
| 311 | |
---|
| 312 | IMPLICIT NONE |
---|
| 313 | |
---|
| 314 | INTEGER(iwp) :: i !< running index x direction |
---|
| 315 | INTEGER(iwp) :: j !< running index y direction |
---|
| 316 | INTEGER(iwp) :: k !< running index z direction |
---|
| 317 | INTEGER(iwp) :: m !< running index surface elements |
---|
| 318 | |
---|
| 319 | REAL(wp) :: pden !< temporary scalar user for calculating density |
---|
| 320 | REAL(wp) :: pnom !< temporary scalar user for calculating density |
---|
| 321 | REAL(wp) :: p1 !< temporary scalar user for calculating density |
---|
| 322 | REAL(wp) :: p2 !< temporary scalar user for calculating density |
---|
| 323 | REAL(wp) :: p3 !< temporary scalar user for calculating density |
---|
| 324 | REAL(wp) :: pt1 !< temporary scalar user for calculating density |
---|
| 325 | REAL(wp) :: pt2 !< temporary scalar user for calculating density |
---|
| 326 | REAL(wp) :: pt3 !< temporary scalar user for calculating density |
---|
| 327 | REAL(wp) :: pt4 !< temporary scalar user for calculating density |
---|
| 328 | REAL(wp) :: sa1 !< temporary scalar user for calculating density |
---|
| 329 | REAL(wp) :: sa15 !< temporary scalar user for calculating density |
---|
| 330 | REAL(wp) :: sa2 !< temporary scalar user for calculating density |
---|
| 331 | |
---|
| 332 | |
---|
| 333 | DO i = nxl, nxr |
---|
| 334 | DO j = nys, nyn |
---|
| 335 | DO k = nzb+1, nzt |
---|
| 336 | ! |
---|
| 337 | !-- Pressure is needed in dbar |
---|
| 338 | p1 = hyp(k) * 1E-4_wp |
---|
| 339 | p2 = p1 * p1 |
---|
| 340 | p3 = p2 * p1 |
---|
| 341 | |
---|
| 342 | ! |
---|
| 343 | !-- Temperature needed in degree Celsius |
---|
| 344 | pt1 = pt_p(k,j,i) - 273.15_wp |
---|
| 345 | pt2 = pt1 * pt1 |
---|
| 346 | pt3 = pt1 * pt2 |
---|
| 347 | pt4 = pt2 * pt2 |
---|
| 348 | |
---|
| 349 | sa1 = sa_p(k,j,i) |
---|
| 350 | sa15 = sa1 * SQRT( sa1 ) |
---|
| 351 | sa2 = sa1 * sa1 |
---|
| 352 | |
---|
[4797] | 353 | pnom = nom(1) + nom(2)*pt1 + nom(3)*pt2 + & |
---|
| 354 | nom(4)*pt3 + nom(5)*sa1 + nom(6)*sa1*pt1 + & |
---|
[3294] | 355 | nom(7)*sa2 |
---|
| 356 | |
---|
[4797] | 357 | pden = den(1) + den(2)*pt1 + den(3)*pt2 + & |
---|
| 358 | den(4)*pt3 + den(5)*pt4 + den(6)*sa1 + & |
---|
| 359 | den(7)*sa1*pt1 + den(8)*sa1*pt3 + den(9)*sa15 + & |
---|
[3294] | 360 | den(10)*sa15*pt2 |
---|
| 361 | ! |
---|
| 362 | !-- Potential density (without pressure terms) |
---|
| 363 | prho(k,j,i) = pnom / pden |
---|
| 364 | |
---|
[4797] | 365 | pnom = pnom + nom(8)*p1 + nom(9)*p1*pt2 + & |
---|
[3294] | 366 | nom(10)*p1*sa1 + nom(11)*p2 + nom(12)*p2*pt2 |
---|
| 367 | |
---|
[4797] | 368 | pden = pden + den(11)*p1 + den(12)*p2*pt3 + & |
---|
[3294] | 369 | den(13)*p3*pt1 |
---|
| 370 | |
---|
| 371 | ! |
---|
| 372 | !-- In-situ density |
---|
| 373 | rho_ocean(k,j,i) = pnom / pden |
---|
| 374 | |
---|
| 375 | ENDDO |
---|
| 376 | |
---|
| 377 | ! |
---|
[4797] | 378 | !-- Neumann conditions are assumed at top boundary and currently also at bottom boundary |
---|
| 379 | !-- (see comment lines below) |
---|
[3294] | 380 | prho(nzt+1,j,i) = prho(nzt,j,i) |
---|
| 381 | rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i) |
---|
| 382 | |
---|
| 383 | ENDDO |
---|
| 384 | ENDDO |
---|
| 385 | ! |
---|
| 386 | !-- Neumann conditions at up/downward-facing surfaces |
---|
| 387 | !$OMP PARALLEL DO PRIVATE( i, j, k ) |
---|
| 388 | DO m = 1, bc_h(0)%ns |
---|
| 389 | i = bc_h(0)%i(m) |
---|
| 390 | j = bc_h(0)%j(m) |
---|
| 391 | k = bc_h(0)%k(m) |
---|
| 392 | prho(k-1,j,i) = prho(k,j,i) |
---|
| 393 | rho_ocean(k-1,j,i) = rho_ocean(k,j,i) |
---|
| 394 | ENDDO |
---|
| 395 | ! |
---|
| 396 | !-- Downward facing surfaces |
---|
| 397 | !$OMP PARALLEL DO PRIVATE( i, j, k ) |
---|
| 398 | DO m = 1, bc_h(1)%ns |
---|
| 399 | i = bc_h(1)%i(m) |
---|
| 400 | j = bc_h(1)%j(m) |
---|
| 401 | k = bc_h(1)%k(m) |
---|
| 402 | prho(k+1,j,i) = prho(k,j,i) |
---|
| 403 | rho_ocean(k+1,j,i) = rho_ocean(k,j,i) |
---|
| 404 | ENDDO |
---|
| 405 | |
---|
| 406 | END SUBROUTINE eqn_state_seawater |
---|
| 407 | |
---|
| 408 | |
---|
[4797] | 409 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 410 | ! Description: |
---|
| 411 | ! ------------ |
---|
| 412 | !> Same as above, but for grid point i,j |
---|
[4797] | 413 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 414 | SUBROUTINE eqn_state_seawater_ij( i, j ) |
---|
| 415 | |
---|
[4797] | 416 | USE arrays_3d, & |
---|
[3294] | 417 | ONLY: hyp, prho, pt_p, rho_ocean, sa_p |
---|
| 418 | |
---|
[4797] | 419 | USE indices, & |
---|
[3294] | 420 | ONLY: nzb, nzt |
---|
| 421 | |
---|
| 422 | IMPLICIT NONE |
---|
| 423 | |
---|
| 424 | INTEGER(iwp) :: i !< running index x direction |
---|
| 425 | INTEGER(iwp) :: j !< running index y direction |
---|
| 426 | INTEGER(iwp) :: k !< running index z direction |
---|
| 427 | INTEGER(iwp) :: m !< running index surface elements |
---|
| 428 | INTEGER(iwp) :: surf_e !< end index of surface elements at (j,i)-gridpoint |
---|
| 429 | INTEGER(iwp) :: surf_s !< start index of surface elements at (j,i)-gridpoint |
---|
| 430 | |
---|
| 431 | REAL(wp) :: pden !< temporary scalar user for calculating density |
---|
| 432 | REAL(wp) :: pnom !< temporary scalar user for calculating density |
---|
| 433 | REAL(wp) :: p1 !< temporary scalar user for calculating density |
---|
| 434 | REAL(wp) :: p2 !< temporary scalar user for calculating density |
---|
| 435 | REAL(wp) :: p3 !< temporary scalar user for calculating density |
---|
| 436 | REAL(wp) :: pt1 !< temporary scalar user for calculating density |
---|
| 437 | REAL(wp) :: pt2 !< temporary scalar user for calculating density |
---|
| 438 | REAL(wp) :: pt3 !< temporary scalar user for calculating density |
---|
| 439 | REAL(wp) :: pt4 !< temporary scalar user for calculating density |
---|
| 440 | REAL(wp) :: sa1 !< temporary scalar user for calculating density |
---|
| 441 | REAL(wp) :: sa15 !< temporary scalar user for calculating density |
---|
| 442 | REAL(wp) :: sa2 !< temporary scalar user for calculating density |
---|
| 443 | |
---|
[4797] | 444 | |
---|
[3294] | 445 | DO k = nzb+1, nzt |
---|
| 446 | ! |
---|
| 447 | !-- Pressure is needed in dbar |
---|
| 448 | p1 = hyp(k) * 1E-4_wp |
---|
| 449 | p2 = p1 * p1 |
---|
| 450 | p3 = p2 * p1 |
---|
| 451 | ! |
---|
| 452 | !-- Temperature needed in degree Celsius |
---|
| 453 | pt1 = pt_p(k,j,i) - 273.15_wp |
---|
| 454 | pt2 = pt1 * pt1 |
---|
| 455 | pt3 = pt1 * pt2 |
---|
| 456 | pt4 = pt2 * pt2 |
---|
| 457 | |
---|
| 458 | sa1 = sa_p(k,j,i) |
---|
| 459 | sa15 = sa1 * SQRT( sa1 ) |
---|
| 460 | sa2 = sa1 * sa1 |
---|
| 461 | |
---|
[4797] | 462 | pnom = nom(1) + nom(2)*pt1 + nom(3)*pt2 + & |
---|
[3294] | 463 | nom(4)*pt3 + nom(5)*sa1 + nom(6)*sa1*pt1 + nom(7)*sa2 |
---|
| 464 | |
---|
[4797] | 465 | pden = den(1) + den(2)*pt1 + den(3)*pt2 + & |
---|
| 466 | den(4)*pt3 + den(5)*pt4 + den(6)*sa1 + & |
---|
| 467 | den(7)*sa1*pt1 + den(8)*sa1*pt3 + den(9)*sa15 + & |
---|
[3294] | 468 | den(10)*sa15*pt2 |
---|
| 469 | ! |
---|
| 470 | !-- Potential density (without pressure terms) |
---|
| 471 | prho(k,j,i) = pnom / pden |
---|
| 472 | |
---|
[4797] | 473 | pnom = pnom + nom(8)*p1 + nom(9)*p1*pt2 + & |
---|
[3294] | 474 | nom(10)*p1*sa1 + nom(11)*p2 + nom(12)*p2*pt2 |
---|
[4797] | 475 | |
---|
| 476 | pden = pden + den(11)*p1 + den(12)*p2*pt3 + & |
---|
[3294] | 477 | den(13)*p3*pt1 |
---|
| 478 | |
---|
| 479 | ! |
---|
| 480 | !-- In-situ density |
---|
| 481 | rho_ocean(k,j,i) = pnom / pden |
---|
| 482 | |
---|
| 483 | ENDDO |
---|
| 484 | ! |
---|
| 485 | !-- Neumann conditions at up/downward-facing walls |
---|
| 486 | surf_s = bc_h(0)%start_index(j,i) |
---|
| 487 | surf_e = bc_h(0)%end_index(j,i) |
---|
| 488 | DO m = surf_s, surf_e |
---|
| 489 | k = bc_h(0)%k(m) |
---|
| 490 | prho(k-1,j,i) = prho(k,j,i) |
---|
| 491 | rho_ocean(k-1,j,i) = rho_ocean(k,j,i) |
---|
| 492 | ENDDO |
---|
| 493 | ! |
---|
| 494 | !-- Downward facing surfaces |
---|
| 495 | surf_s = bc_h(1)%start_index(j,i) |
---|
| 496 | surf_e = bc_h(1)%end_index(j,i) |
---|
| 497 | DO m = surf_s, surf_e |
---|
| 498 | k = bc_h(1)%k(m) |
---|
| 499 | prho(k+1,j,i) = prho(k,j,i) |
---|
| 500 | rho_ocean(k+1,j,i) = rho_ocean(k,j,i) |
---|
| 501 | ENDDO |
---|
| 502 | ! |
---|
| 503 | !-- Neumann condition are assumed at top boundary |
---|
| 504 | prho(nzt+1,j,i) = prho(nzt,j,i) |
---|
| 505 | rho_ocean(nzt+1,j,i) = rho_ocean(nzt,j,i) |
---|
| 506 | |
---|
| 507 | END SUBROUTINE eqn_state_seawater_ij |
---|
| 508 | |
---|
| 509 | |
---|
[4797] | 510 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 511 | ! Description: |
---|
| 512 | ! ------------ |
---|
| 513 | !> Equation of state for seawater as a function |
---|
[4797] | 514 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 515 | REAL(wp) FUNCTION eqn_state_seawater_func( p, pt, sa ) |
---|
| 516 | |
---|
| 517 | IMPLICIT NONE |
---|
| 518 | |
---|
| 519 | REAL(wp) :: p !< temporary scalar user for calculating density |
---|
| 520 | REAL(wp) :: p1 !< temporary scalar user for calculating density |
---|
| 521 | REAL(wp) :: p2 !< temporary scalar user for calculating density |
---|
| 522 | REAL(wp) :: p3 !< temporary scalar user for calculating density |
---|
| 523 | REAL(wp) :: pt !< temporary scalar user for calculating density |
---|
| 524 | REAL(wp) :: pt1 !< temporary scalar user for calculating density |
---|
| 525 | REAL(wp) :: pt2 !< temporary scalar user for calculating density |
---|
| 526 | REAL(wp) :: pt3 !< temporary scalar user for calculating density |
---|
| 527 | REAL(wp) :: pt4 !< temporary scalar user for calculating density |
---|
| 528 | REAL(wp) :: sa !< temporary scalar user for calculating density |
---|
| 529 | REAL(wp) :: sa15 !< temporary scalar user for calculating density |
---|
| 530 | REAL(wp) :: sa2 !< temporary scalar user for calculating density |
---|
| 531 | |
---|
[4797] | 532 | |
---|
[3294] | 533 | ! |
---|
| 534 | !-- Pressure is needed in dbar |
---|
| 535 | p1 = p * 1.0E-4_wp |
---|
| 536 | p2 = p1 * p1 |
---|
| 537 | p3 = p2 * p1 |
---|
| 538 | |
---|
| 539 | ! |
---|
| 540 | !-- Temperature needed in degree Celsius |
---|
| 541 | pt1 = pt - 273.15_wp |
---|
| 542 | pt2 = pt1 * pt1 |
---|
| 543 | pt3 = pt1 * pt2 |
---|
| 544 | pt4 = pt2 * pt2 |
---|
| 545 | |
---|
| 546 | sa15 = sa * SQRT( sa ) |
---|
| 547 | sa2 = sa * sa |
---|
| 548 | |
---|
| 549 | |
---|
[4797] | 550 | eqn_state_seawater_func = & |
---|
| 551 | ( nom(1) + nom(2)*pt1 + nom(3)*pt2 + nom(4)*pt3 + & |
---|
| 552 | nom(5)*sa + nom(6)*sa*pt1 + nom(7)*sa2 + nom(8)*p1 + & |
---|
| 553 | nom(9)*p1*pt2 + nom(10)*p1*sa + nom(11)*p2 + nom(12)*p2*pt2 & |
---|
| 554 | ) / & |
---|
| 555 | ( den(1) + den(2)*pt1 + den(3)*pt2 + den(4)*pt3 + & |
---|
| 556 | den(5)*pt4 + den(6)*sa + den(7)*sa*pt1 + den(8)*sa*pt3 + & |
---|
| 557 | den(9)*sa15 + den(10)*sa15*pt2 + den(11)*p1 + den(12)*p2*pt3 + & |
---|
| 558 | den(13)*p3*pt1 & |
---|
| 559 | ) |
---|
[3294] | 560 | |
---|
| 561 | |
---|
| 562 | END FUNCTION eqn_state_seawater_func |
---|
| 563 | |
---|
| 564 | |
---|
[4797] | 565 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 566 | ! Description: |
---|
| 567 | ! ------------ |
---|
| 568 | !> Reads the ocean parameters namelist |
---|
[4797] | 569 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 570 | SUBROUTINE ocean_parin |
---|
| 571 | |
---|
| 572 | IMPLICIT NONE |
---|
| 573 | |
---|
| 574 | CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file |
---|
| 575 | |
---|
| 576 | |
---|
[4797] | 577 | NAMELIST /ocean_parameters/ bc_sa_t, bottom_salinityflux, salinity, sa_surface, & |
---|
| 578 | sa_vertical_gradient, sa_vertical_gradient_level, & |
---|
| 579 | stokes_waveheight, stokes_wavelength, surface_cooling_spinup_time,& |
---|
| 580 | top_salinityflux, wall_salinityflux, wave_breaking |
---|
[3294] | 581 | |
---|
| 582 | ! |
---|
| 583 | !-- Try to find the namelist |
---|
| 584 | REWIND ( 11 ) |
---|
| 585 | line = ' ' |
---|
| 586 | DO WHILE ( INDEX( line, '&ocean_parameters' ) == 0 ) |
---|
| 587 | READ ( 11, '(A)', END=12 ) line |
---|
| 588 | ENDDO |
---|
| 589 | BACKSPACE ( 11 ) |
---|
| 590 | |
---|
| 591 | ! |
---|
| 592 | !-- Read namelist |
---|
| 593 | READ ( 11, ocean_parameters, ERR = 10 ) |
---|
| 594 | ! |
---|
| 595 | !-- Set switch that enables PALM's ocean mode |
---|
| 596 | ocean_mode = .TRUE. |
---|
| 597 | |
---|
| 598 | GOTO 12 |
---|
| 599 | |
---|
| 600 | 10 BACKSPACE( 11 ) |
---|
| 601 | READ( 11 , '(A)') line |
---|
| 602 | CALL parin_fail_message( 'ocean_parameters', line ) |
---|
| 603 | |
---|
| 604 | 12 CONTINUE |
---|
| 605 | |
---|
| 606 | END SUBROUTINE ocean_parin |
---|
| 607 | |
---|
[4797] | 608 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 609 | ! Description: |
---|
| 610 | ! ------------ |
---|
| 611 | !> Check parameters routine for the ocean mode |
---|
[4797] | 612 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 613 | SUBROUTINE ocean_check_parameters |
---|
| 614 | |
---|
[4797] | 615 | USE control_parameters, & |
---|
| 616 | ONLY: coupling_char, coupling_mode, initializing_actions, message_string, use_top_fluxes |
---|
[3294] | 617 | |
---|
[4797] | 618 | USE pmc_interface, & |
---|
[3311] | 619 | ONLY: nested_run |
---|
| 620 | |
---|
[3294] | 621 | IMPLICIT NONE |
---|
| 622 | |
---|
| 623 | |
---|
| 624 | ! |
---|
[3311] | 625 | !-- Check for invalid combinations |
---|
| 626 | IF ( nested_run ) THEN |
---|
| 627 | message_string = 'ocean mode not allowed for nesting' |
---|
| 628 | CALL message( 'ocean_check_parameters', 'PA0510', 1, 2, 0, 6, 0 ) |
---|
| 629 | ENDIF |
---|
| 630 | |
---|
| 631 | IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN |
---|
| 632 | message_string = 'ocean mode does not allow cyclic-fill initialization' |
---|
| 633 | CALL message( 'ocean_check_parameters', 'PA0511', 1, 2, 0, 6, 0 ) |
---|
| 634 | ENDIF |
---|
| 635 | |
---|
| 636 | ! |
---|
[3294] | 637 | !-- Check ocean setting |
---|
[4797] | 638 | IF ( TRIM( coupling_mode ) == 'uncoupled' .AND. TRIM( coupling_char ) == '_O' .AND. & |
---|
[3294] | 639 | .NOT. ocean_mode ) THEN |
---|
| 640 | |
---|
| 641 | ! |
---|
[4797] | 642 | !-- Check whether an (uncoupled) atmospheric run has been declared as an ocean run (this setting |
---|
| 643 | !-- is done via palmrun-option -y) |
---|
| 644 | message_string = 'ocean mode does not allow coupling_char = "' // TRIM( coupling_char ) // & |
---|
| 645 | '" set by palmrun-option "-y"' |
---|
[3311] | 646 | CALL message( 'ocean_check_parameters', 'PA0317', 1, 2, 0, 6, 0 ) |
---|
[3294] | 647 | |
---|
| 648 | ENDIF |
---|
| 649 | |
---|
| 650 | ! |
---|
| 651 | !-- Ocean version must use flux boundary conditions at the top |
---|
| 652 | IF ( .NOT. use_top_fluxes ) THEN |
---|
| 653 | message_string = 'use_top_fluxes must be .TRUE. in ocean mode' |
---|
| 654 | CALL message( 'ocean_check_parameters', 'PA0042', 1, 2, 0, 6, 0 ) |
---|
| 655 | ENDIF |
---|
| 656 | |
---|
| 657 | ! |
---|
| 658 | !-- Boundary conditions for salinity |
---|
| 659 | IF ( bc_sa_t == 'dirichlet' ) THEN |
---|
| 660 | ibc_sa_t = 0 |
---|
| 661 | ELSEIF ( bc_sa_t == 'neumann' ) THEN |
---|
| 662 | ibc_sa_t = 1 |
---|
| 663 | ELSE |
---|
[4797] | 664 | message_string = 'unknown boundary condition: bc_sa_t = "' // TRIM( bc_sa_t ) // '"' |
---|
[3294] | 665 | CALL message( 'ocean_check_parameters', 'PA0068', 1, 2, 0, 6, 0 ) |
---|
| 666 | ENDIF |
---|
| 667 | |
---|
| 668 | IF ( top_salinityflux == 9999999.9_wp ) constant_top_salinityflux = .FALSE. |
---|
| 669 | |
---|
[3303] | 670 | IF ( .NOT. salinity ) THEN |
---|
[4797] | 671 | IF ( ( bottom_salinityflux /= 0.0_wp .AND. bottom_salinityflux /= 9999999.9_wp ) .OR. & |
---|
| 672 | ( top_salinityflux /= 0.0_wp .AND. top_salinityflux /= 9999999.9_wp ) ) & |
---|
[3303] | 673 | THEN |
---|
[4797] | 674 | message_string = 'salinityflux must not be set for ocean run ' // 'without salinity' |
---|
[3311] | 675 | CALL message( 'ocean_check_parameters', 'PA0509', 1, 2, 0, 6, 0 ) |
---|
[3303] | 676 | ENDIF |
---|
| 677 | ENDIF |
---|
| 678 | |
---|
[3294] | 679 | IF ( ibc_sa_t == 1 .AND. top_salinityflux == 9999999.9_wp ) THEN |
---|
[4797] | 680 | message_string = 'boundary condition: bc_sa_t = "' // TRIM( bc_sa_t ) // & |
---|
| 681 | '" requires to set top_salinityflux' |
---|
[3294] | 682 | CALL message( 'ocean_check_parameters', 'PA0069', 1, 2, 0, 6, 0 ) |
---|
| 683 | ENDIF |
---|
| 684 | |
---|
| 685 | ! |
---|
[4797] | 686 | !-- A fixed salinity at the top implies Dirichlet boundary condition for salinity. In this case |
---|
| 687 | !-- specification of a constant salinity flux is forbidden. |
---|
| 688 | IF ( ibc_sa_t == 0 .AND. constant_top_salinityflux .AND. top_salinityflux /= 0.0_wp ) THEN |
---|
| 689 | message_string = 'boundary condition: bc_sa_t = "' // TRIM( bc_sa_t ) // & |
---|
| 690 | '" is not allowed with top_salinityflux /= 0.0' |
---|
[3294] | 691 | CALL message( 'ocean_check_parameters', 'PA0070', 1, 2, 0, 6, 0 ) |
---|
| 692 | ENDIF |
---|
| 693 | |
---|
[3302] | 694 | ! |
---|
| 695 | !-- Check if Stokes force is to be used |
---|
| 696 | IF ( stokes_waveheight /= 0.0_wp .AND. stokes_wavelength /= 0.0_wp ) THEN |
---|
| 697 | stokes_force = .TRUE. |
---|
| 698 | ELSE |
---|
[4797] | 699 | IF ( ( stokes_waveheight <= 0.0_wp .AND. stokes_wavelength > 0.0_wp ) .OR. & |
---|
| 700 | ( stokes_waveheight > 0.0_wp .AND. stokes_wavelength <= 0.0_wp ) .OR. & |
---|
| 701 | ( stokes_waveheight < 0.0_wp .AND. stokes_wavelength < 0.0_wp ) ) & |
---|
[3302] | 702 | THEN |
---|
[4797] | 703 | message_string = 'wrong settings for stokes_wavelength and/or stokes_waveheight' |
---|
[3302] | 704 | CALL message( 'ocean_check_parameters', 'PA0460', 1, 2, 0, 6, 0 ) |
---|
| 705 | ENDIF |
---|
| 706 | ENDIF |
---|
| 707 | |
---|
[3294] | 708 | END SUBROUTINE ocean_check_parameters |
---|
| 709 | |
---|
| 710 | |
---|
[4797] | 711 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 712 | ! Description: |
---|
| 713 | ! ------------ |
---|
| 714 | !> Check data output. |
---|
[4797] | 715 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 716 | SUBROUTINE ocean_check_data_output( var, unit ) |
---|
[4797] | 717 | |
---|
[3294] | 718 | IMPLICIT NONE |
---|
| 719 | |
---|
| 720 | CHARACTER (LEN=*) :: unit !< unit of output variable |
---|
| 721 | CHARACTER (LEN=*) :: var !< name of output variable |
---|
| 722 | |
---|
| 723 | |
---|
| 724 | SELECT CASE ( TRIM( var ) ) |
---|
| 725 | |
---|
[3421] | 726 | CASE ( 'rho_sea_water' ) |
---|
[3294] | 727 | unit = 'kg/m3' |
---|
| 728 | |
---|
| 729 | CASE ( 'sa' ) |
---|
| 730 | unit = 'psu' |
---|
| 731 | |
---|
| 732 | CASE DEFAULT |
---|
| 733 | unit = 'illegal' |
---|
| 734 | |
---|
| 735 | END SELECT |
---|
| 736 | |
---|
| 737 | END SUBROUTINE ocean_check_data_output |
---|
| 738 | |
---|
| 739 | |
---|
[4797] | 740 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 741 | ! Description: |
---|
| 742 | ! ------------ |
---|
| 743 | !> Check data output of profiles |
---|
[4797] | 744 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 745 | SUBROUTINE ocean_check_data_output_pr( variable, var_count, unit, dopr_unit ) |
---|
| 746 | |
---|
[4797] | 747 | USE arrays_3d, & |
---|
[3294] | 748 | ONLY: zu, zw |
---|
| 749 | |
---|
[4797] | 750 | USE control_parameters, & |
---|
[3614] | 751 | ONLY: data_output_pr |
---|
[3294] | 752 | |
---|
| 753 | USE indices |
---|
| 754 | |
---|
| 755 | USE profil_parameter |
---|
| 756 | |
---|
| 757 | USE statistics |
---|
| 758 | |
---|
| 759 | IMPLICIT NONE |
---|
| 760 | |
---|
[4797] | 761 | CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit |
---|
[3294] | 762 | CHARACTER (LEN=*) :: unit !< |
---|
| 763 | CHARACTER (LEN=*) :: variable !< |
---|
| 764 | |
---|
| 765 | INTEGER(iwp) :: var_count !< |
---|
| 766 | |
---|
[4797] | 767 | |
---|
[3294] | 768 | SELECT CASE ( TRIM( variable ) ) |
---|
| 769 | |
---|
| 770 | CASE ( 'prho' ) |
---|
| 771 | dopr_index(var_count) = 71 |
---|
| 772 | dopr_unit = 'kg/m3' |
---|
| 773 | hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 ) |
---|
| 774 | unit = dopr_unit |
---|
| 775 | |
---|
[3421] | 776 | CASE ( 'rho_sea_water' ) |
---|
[3294] | 777 | dopr_index(var_count) = 64 |
---|
| 778 | dopr_unit = 'kg/m3' |
---|
| 779 | hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 ) |
---|
| 780 | IF ( data_output_pr(var_count)(1:1) == '#' ) THEN |
---|
| 781 | dopr_initial_index(var_count) = 77 |
---|
| 782 | hom(:,2,77,:) = SPREAD( zu, 2, statistic_regions+1 ) |
---|
| 783 | hom(nzb,2,77,:) = 0.0_wp ! because zu(nzb) is negative |
---|
| 784 | data_output_pr(var_count) = data_output_pr(var_count)(2:) |
---|
| 785 | ENDIF |
---|
| 786 | unit = dopr_unit |
---|
| 787 | |
---|
| 788 | CASE ( 'sa', '#sa' ) |
---|
| 789 | dopr_index(var_count) = 23 |
---|
| 790 | dopr_unit = 'psu' |
---|
| 791 | hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 ) |
---|
| 792 | IF ( data_output_pr(var_count)(1:1) == '#' ) THEN |
---|
| 793 | dopr_initial_index(var_count) = 26 |
---|
| 794 | hom(:,2,26,:) = SPREAD( zu, 2, statistic_regions+1 ) |
---|
| 795 | hom(nzb,2,26,:) = 0.0_wp ! because zu(nzb) is negative |
---|
| 796 | data_output_pr(var_count) = data_output_pr(var_count)(2:) |
---|
| 797 | ENDIF |
---|
| 798 | unit = dopr_unit |
---|
| 799 | |
---|
| 800 | CASE ( 'w"sa"' ) |
---|
| 801 | dopr_index(var_count) = 65 |
---|
| 802 | dopr_unit = 'psu m/s' |
---|
| 803 | hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 ) |
---|
| 804 | unit = dopr_unit |
---|
| 805 | |
---|
| 806 | CASE ( 'w*sa*' ) |
---|
| 807 | dopr_index(var_count) = 66 |
---|
| 808 | dopr_unit = 'psu m/s' |
---|
| 809 | hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 ) |
---|
| 810 | unit = dopr_unit |
---|
| 811 | |
---|
| 812 | CASE ( 'wsa' ) |
---|
| 813 | dopr_index(var_count) = 67 |
---|
| 814 | dopr_unit = 'psu m/s' |
---|
| 815 | hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 ) |
---|
| 816 | unit = dopr_unit |
---|
| 817 | |
---|
| 818 | CASE DEFAULT |
---|
| 819 | unit = 'illegal' |
---|
| 820 | |
---|
| 821 | END SELECT |
---|
| 822 | |
---|
| 823 | |
---|
| 824 | END SUBROUTINE ocean_check_data_output_pr |
---|
| 825 | |
---|
| 826 | |
---|
[4797] | 827 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 828 | ! Description: |
---|
| 829 | ! ------------ |
---|
| 830 | !> Define appropriate grid for netcdf variables. |
---|
| 831 | !> It is called out from subroutine netcdf. |
---|
[4797] | 832 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 833 | SUBROUTINE ocean_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) |
---|
[4797] | 834 | |
---|
[3294] | 835 | IMPLICIT NONE |
---|
| 836 | |
---|
| 837 | CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< x grid of output variable |
---|
| 838 | CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< y grid of output variable |
---|
| 839 | CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< z grid of output variable |
---|
| 840 | CHARACTER (LEN=*), INTENT(IN) :: var !< name of output variable |
---|
| 841 | |
---|
| 842 | LOGICAL, INTENT(OUT) :: found !< flag if output variable is found |
---|
| 843 | |
---|
[4797] | 844 | |
---|
[3294] | 845 | found = .TRUE. |
---|
| 846 | |
---|
| 847 | ! |
---|
| 848 | !-- Check for the grid |
---|
| 849 | SELECT CASE ( TRIM( var ) ) |
---|
| 850 | |
---|
[4797] | 851 | CASE ( 'rho_sea_water', 'rho_sea_water_xy', 'rho_sea_water_xz', 'rho_sea_water_yz', 'sa', & |
---|
| 852 | 'sa_xy', 'sa_xz', 'sa_yz' ) |
---|
[3294] | 853 | grid_x = 'x' |
---|
| 854 | grid_y = 'y' |
---|
| 855 | grid_z = 'zu' |
---|
| 856 | |
---|
| 857 | CASE DEFAULT |
---|
| 858 | found = .FALSE. |
---|
| 859 | grid_x = 'none' |
---|
| 860 | grid_y = 'none' |
---|
| 861 | grid_z = 'none' |
---|
| 862 | |
---|
| 863 | END SELECT |
---|
| 864 | |
---|
| 865 | END SUBROUTINE ocean_define_netcdf_grid |
---|
| 866 | |
---|
| 867 | |
---|
[4797] | 868 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 869 | ! Description: |
---|
| 870 | ! ------------ |
---|
| 871 | !> Average 3D data. |
---|
[4797] | 872 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 873 | SUBROUTINE ocean_3d_data_averaging( mode, variable ) |
---|
| 874 | |
---|
[4797] | 875 | |
---|
| 876 | USE arrays_3d, & |
---|
[3294] | 877 | ONLY: rho_ocean, sa |
---|
| 878 | |
---|
[4797] | 879 | USE averaging, & |
---|
[3294] | 880 | ONLY: rho_ocean_av, sa_av |
---|
| 881 | |
---|
[4797] | 882 | USE control_parameters, & |
---|
[3294] | 883 | ONLY: average_count_3d |
---|
| 884 | |
---|
[4797] | 885 | USE indices, & |
---|
[3294] | 886 | ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt |
---|
| 887 | |
---|
| 888 | IMPLICIT NONE |
---|
| 889 | |
---|
| 890 | CHARACTER (LEN=*) :: mode !< flag defining mode 'allocate', 'sum' or 'average' |
---|
| 891 | CHARACTER (LEN=*) :: variable !< name of variable |
---|
| 892 | |
---|
| 893 | INTEGER(iwp) :: i !< loop index |
---|
| 894 | INTEGER(iwp) :: j !< loop index |
---|
| 895 | INTEGER(iwp) :: k !< loop index |
---|
| 896 | |
---|
[4797] | 897 | |
---|
[3294] | 898 | IF ( mode == 'allocate' ) THEN |
---|
| 899 | |
---|
| 900 | SELECT CASE ( TRIM( variable ) ) |
---|
| 901 | |
---|
[3421] | 902 | CASE ( 'rho_sea_water' ) |
---|
[3294] | 903 | IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN |
---|
| 904 | ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 905 | ENDIF |
---|
| 906 | rho_ocean_av = 0.0_wp |
---|
| 907 | |
---|
| 908 | CASE ( 'sa' ) |
---|
| 909 | IF ( .NOT. ALLOCATED( sa_av ) ) THEN |
---|
| 910 | ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 911 | ENDIF |
---|
| 912 | sa_av = 0.0_wp |
---|
| 913 | |
---|
| 914 | CASE DEFAULT |
---|
| 915 | CONTINUE |
---|
| 916 | |
---|
| 917 | END SELECT |
---|
| 918 | |
---|
| 919 | ELSEIF ( mode == 'sum' ) THEN |
---|
| 920 | |
---|
| 921 | SELECT CASE ( TRIM( variable ) ) |
---|
| 922 | |
---|
[3421] | 923 | CASE ( 'rho_sea_water' ) |
---|
[3294] | 924 | IF ( ALLOCATED( rho_ocean_av ) ) THEN |
---|
| 925 | DO i = nxlg, nxrg |
---|
| 926 | DO j = nysg, nyng |
---|
| 927 | DO k = nzb, nzt+1 |
---|
[4797] | 928 | rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) + rho_ocean(k,j,i) |
---|
[3294] | 929 | ENDDO |
---|
| 930 | ENDDO |
---|
| 931 | ENDDO |
---|
| 932 | ENDIF |
---|
| 933 | |
---|
| 934 | CASE ( 'sa' ) |
---|
| 935 | IF ( ALLOCATED( sa_av ) ) THEN |
---|
| 936 | DO i = nxlg, nxrg |
---|
| 937 | DO j = nysg, nyng |
---|
| 938 | DO k = nzb, nzt+1 |
---|
| 939 | sa_av(k,j,i) = sa_av(k,j,i) + sa(k,j,i) |
---|
| 940 | ENDDO |
---|
| 941 | ENDDO |
---|
| 942 | ENDDO |
---|
| 943 | ENDIF |
---|
| 944 | |
---|
| 945 | CASE DEFAULT |
---|
| 946 | CONTINUE |
---|
| 947 | |
---|
| 948 | END SELECT |
---|
| 949 | |
---|
| 950 | ELSEIF ( mode == 'average' ) THEN |
---|
| 951 | |
---|
| 952 | SELECT CASE ( TRIM( variable ) ) |
---|
| 953 | |
---|
[3421] | 954 | CASE ( 'rho_sea_water' ) |
---|
[4797] | 955 | IF ( ALLOCATED( rho_ocean_av ) ) THEN |
---|
[3294] | 956 | DO i = nxlg, nxrg |
---|
| 957 | DO j = nysg, nyng |
---|
| 958 | DO k = nzb, nzt+1 |
---|
[4797] | 959 | rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) / & |
---|
[3294] | 960 | REAL( average_count_3d, KIND=wp ) |
---|
| 961 | ENDDO |
---|
| 962 | ENDDO |
---|
| 963 | ENDDO |
---|
| 964 | ENDIF |
---|
| 965 | |
---|
| 966 | CASE ( 'sa' ) |
---|
[4797] | 967 | IF ( ALLOCATED( sa_av ) ) THEN |
---|
[3294] | 968 | DO i = nxlg, nxrg |
---|
| 969 | DO j = nysg, nyng |
---|
| 970 | DO k = nzb, nzt+1 |
---|
[4797] | 971 | sa_av(k,j,i) = sa_av(k,j,i) / REAL( average_count_3d, KIND=wp ) |
---|
[3294] | 972 | ENDDO |
---|
| 973 | ENDDO |
---|
| 974 | ENDDO |
---|
| 975 | ENDIF |
---|
| 976 | |
---|
| 977 | END SELECT |
---|
| 978 | |
---|
| 979 | ENDIF |
---|
| 980 | |
---|
| 981 | END SUBROUTINE ocean_3d_data_averaging |
---|
| 982 | |
---|
| 983 | |
---|
[4797] | 984 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 985 | ! Description: |
---|
| 986 | ! ------------ |
---|
| 987 | !> Define 2D output variables. |
---|
[4797] | 988 | !--------------------------------------------------------------------------------------------------! |
---|
| 989 | SUBROUTINE ocean_data_output_2d( av, variable, found, grid, mode, local_pf, nzb_do, nzt_do ) |
---|
| 990 | |
---|
| 991 | USE arrays_3d, & |
---|
[3294] | 992 | ONLY: rho_ocean, sa |
---|
| 993 | |
---|
[4797] | 994 | USE averaging, & |
---|
[3294] | 995 | ONLY: rho_ocean_av, sa_av |
---|
| 996 | |
---|
[4797] | 997 | USE indices, & |
---|
| 998 | ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_total_0 |
---|
[3294] | 999 | |
---|
| 1000 | IMPLICIT NONE |
---|
| 1001 | |
---|
| 1002 | CHARACTER (LEN=*) :: grid !< name of vertical grid |
---|
| 1003 | CHARACTER (LEN=*) :: mode !< either 'xy', 'xz' or 'yz' |
---|
| 1004 | CHARACTER (LEN=*) :: variable !< name of variable |
---|
| 1005 | |
---|
| 1006 | INTEGER(iwp) :: av !< flag for (non-)average output |
---|
| 1007 | INTEGER(iwp) :: flag_nr !< number of masking flag |
---|
| 1008 | INTEGER(iwp) :: i !< loop index |
---|
| 1009 | INTEGER(iwp) :: j !< loop index |
---|
| 1010 | INTEGER(iwp) :: k !< loop index |
---|
| 1011 | INTEGER(iwp) :: nzb_do !< vertical output index (bottom) |
---|
| 1012 | INTEGER(iwp) :: nzt_do !< vertical output index (top) |
---|
| 1013 | |
---|
| 1014 | LOGICAL :: found !< flag if output variable is found |
---|
| 1015 | LOGICAL :: resorted !< flag if output is already resorted |
---|
| 1016 | |
---|
| 1017 | REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute |
---|
| 1018 | |
---|
[4797] | 1019 | REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local array to which output data is resorted to |
---|
[3294] | 1020 | |
---|
| 1021 | REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable |
---|
[4797] | 1022 | |
---|
| 1023 | |
---|
[3294] | 1024 | found = .TRUE. |
---|
| 1025 | resorted = .FALSE. |
---|
| 1026 | ! |
---|
| 1027 | !-- Set masking flag for topography for not resorted arrays |
---|
| 1028 | flag_nr = 0 |
---|
| 1029 | |
---|
| 1030 | SELECT CASE ( TRIM( variable ) ) |
---|
| 1031 | |
---|
[3421] | 1032 | CASE ( 'rho_sea_water_xy', 'rho_sea_water_xz', 'rho_sea_water_yz' ) |
---|
[3294] | 1033 | IF ( av == 0 ) THEN |
---|
| 1034 | to_be_resorted => rho_ocean |
---|
| 1035 | ELSE |
---|
[4797] | 1036 | IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN |
---|
[3294] | 1037 | ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 1038 | rho_ocean_av = REAL( fill_value, KIND = wp ) |
---|
| 1039 | ENDIF |
---|
| 1040 | to_be_resorted => rho_ocean_av |
---|
| 1041 | ENDIF |
---|
| 1042 | |
---|
| 1043 | CASE ( 'sa_xy', 'sa_xz', 'sa_yz' ) |
---|
| 1044 | IF ( av == 0 ) THEN |
---|
| 1045 | to_be_resorted => sa |
---|
| 1046 | ELSE |
---|
[4797] | 1047 | IF ( .NOT. ALLOCATED( sa_av ) ) THEN |
---|
[3294] | 1048 | ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 1049 | sa_av = REAL( fill_value, KIND = wp ) |
---|
| 1050 | ENDIF |
---|
| 1051 | to_be_resorted => sa_av |
---|
| 1052 | ENDIF |
---|
| 1053 | IF ( mode == 'xy' ) grid = 'zu' |
---|
| 1054 | |
---|
| 1055 | CASE DEFAULT |
---|
| 1056 | found = .FALSE. |
---|
| 1057 | grid = 'none' |
---|
| 1058 | |
---|
| 1059 | END SELECT |
---|
| 1060 | |
---|
| 1061 | IF ( found .AND. .NOT. resorted ) THEN |
---|
[3303] | 1062 | DO i = nxl, nxr |
---|
| 1063 | DO j = nys, nyn |
---|
| 1064 | DO k = nzb_do, nzt_do |
---|
[4797] | 1065 | local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, KIND = wp ), & |
---|
| 1066 | BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) |
---|
[3303] | 1067 | ENDDO |
---|
| 1068 | ENDDO |
---|
| 1069 | ENDDO |
---|
| 1070 | resorted = .TRUE. |
---|
[3294] | 1071 | ENDIF |
---|
[4797] | 1072 | |
---|
[3294] | 1073 | END SUBROUTINE ocean_data_output_2d |
---|
| 1074 | |
---|
[4797] | 1075 | |
---|
| 1076 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1077 | ! Description: |
---|
| 1078 | ! ------------ |
---|
| 1079 | !> Define 3D output variables. |
---|
[4797] | 1080 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1081 | SUBROUTINE ocean_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) |
---|
| 1082 | |
---|
[4797] | 1083 | |
---|
| 1084 | USE arrays_3d, & |
---|
[3294] | 1085 | ONLY: rho_ocean, sa |
---|
| 1086 | |
---|
[4797] | 1087 | USE averaging, & |
---|
[3294] | 1088 | ONLY: rho_ocean_av, sa_av |
---|
| 1089 | |
---|
[4797] | 1090 | USE indices, & |
---|
| 1091 | ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_total_0 |
---|
[3294] | 1092 | |
---|
| 1093 | IMPLICIT NONE |
---|
| 1094 | |
---|
| 1095 | CHARACTER (LEN=*) :: variable !< name of variable |
---|
| 1096 | |
---|
| 1097 | INTEGER(iwp) :: av !< flag for (non-)average output |
---|
| 1098 | INTEGER(iwp) :: flag_nr !< number of masking flag |
---|
| 1099 | INTEGER(iwp) :: i !< loop index |
---|
| 1100 | INTEGER(iwp) :: j !< loop index |
---|
| 1101 | INTEGER(iwp) :: k !< loop index |
---|
| 1102 | INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0) |
---|
| 1103 | INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) |
---|
| 1104 | |
---|
| 1105 | LOGICAL :: found !< flag if output variable is found |
---|
| 1106 | LOGICAL :: resorted !< flag if output is already resorted |
---|
| 1107 | |
---|
| 1108 | REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute |
---|
| 1109 | |
---|
[4768] | 1110 | REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local |
---|
[3294] | 1111 | !< array to which output data is resorted to |
---|
| 1112 | |
---|
| 1113 | REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable |
---|
| 1114 | |
---|
[4797] | 1115 | |
---|
[3294] | 1116 | found = .TRUE. |
---|
| 1117 | resorted = .FALSE. |
---|
| 1118 | ! |
---|
| 1119 | !-- Set masking flag for topography for not resorted arrays |
---|
| 1120 | flag_nr = 0 |
---|
| 1121 | |
---|
| 1122 | SELECT CASE ( TRIM( variable ) ) |
---|
| 1123 | |
---|
[3421] | 1124 | CASE ( 'rho_sea_water' ) |
---|
[3294] | 1125 | IF ( av == 0 ) THEN |
---|
| 1126 | to_be_resorted => rho_ocean |
---|
| 1127 | ELSE |
---|
[4797] | 1128 | IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN |
---|
[3294] | 1129 | ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 1130 | rho_ocean_av = REAL( fill_value, KIND = wp ) |
---|
| 1131 | ENDIF |
---|
| 1132 | to_be_resorted => rho_ocean_av |
---|
| 1133 | ENDIF |
---|
| 1134 | |
---|
| 1135 | CASE ( 'sa' ) |
---|
| 1136 | IF ( av == 0 ) THEN |
---|
| 1137 | to_be_resorted => sa |
---|
| 1138 | ELSE |
---|
[4797] | 1139 | IF ( .NOT. ALLOCATED( sa_av ) ) THEN |
---|
[3294] | 1140 | ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 1141 | sa_av = REAL( fill_value, KIND = wp ) |
---|
| 1142 | ENDIF |
---|
| 1143 | to_be_resorted => sa_av |
---|
| 1144 | ENDIF |
---|
| 1145 | |
---|
| 1146 | CASE DEFAULT |
---|
| 1147 | found = .FALSE. |
---|
| 1148 | |
---|
| 1149 | END SELECT |
---|
| 1150 | |
---|
| 1151 | |
---|
| 1152 | IF ( found .AND. .NOT. resorted ) THEN |
---|
| 1153 | DO i = nxl, nxr |
---|
| 1154 | DO j = nys, nyn |
---|
| 1155 | DO k = nzb_do, nzt_do |
---|
[4797] | 1156 | local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, KIND = wp ), & |
---|
| 1157 | BTEST( wall_flags_total_0(k,j,i), flag_nr ) ) |
---|
[3294] | 1158 | ENDDO |
---|
| 1159 | ENDDO |
---|
| 1160 | ENDDO |
---|
| 1161 | resorted = .TRUE. |
---|
| 1162 | ENDIF |
---|
| 1163 | |
---|
| 1164 | END SUBROUTINE ocean_data_output_3d |
---|
| 1165 | |
---|
| 1166 | |
---|
[4731] | 1167 | !--------------------------------------------------------------------------------------------------! |
---|
| 1168 | ! Description: |
---|
| 1169 | ! ------------ |
---|
| 1170 | !> Exchange ghostpoints |
---|
| 1171 | !--------------------------------------------------------------------------------------------------! |
---|
| 1172 | SUBROUTINE ocean_exchange_horiz( location ) |
---|
| 1173 | |
---|
| 1174 | CHARACTER (LEN=*), INTENT(IN) :: location !< call location string |
---|
| 1175 | |
---|
| 1176 | SELECT CASE ( location ) |
---|
| 1177 | |
---|
| 1178 | CASE ( 'before_prognostic_equation' ) |
---|
| 1179 | |
---|
| 1180 | CASE ( 'after_prognostic_equation' ) |
---|
| 1181 | |
---|
| 1182 | CALL exchange_horiz( sa_p, nbgp ) |
---|
| 1183 | CALL exchange_horiz( rho_ocean, nbgp ) |
---|
| 1184 | CALL exchange_horiz( prho, nbgp ) |
---|
| 1185 | |
---|
| 1186 | END SELECT |
---|
| 1187 | |
---|
| 1188 | END SUBROUTINE ocean_exchange_horiz |
---|
| 1189 | |
---|
| 1190 | |
---|
| 1191 | |
---|
[4797] | 1192 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1193 | ! Description: |
---|
| 1194 | ! ------------ |
---|
[3302] | 1195 | !> Header output for ocean parameters |
---|
[4797] | 1196 | !--------------------------------------------------------------------------------------------------! |
---|
[3302] | 1197 | SUBROUTINE ocean_header( io ) |
---|
| 1198 | |
---|
| 1199 | IMPLICIT NONE |
---|
| 1200 | |
---|
| 1201 | INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file |
---|
| 1202 | |
---|
[4797] | 1203 | |
---|
[3302] | 1204 | ! |
---|
| 1205 | !-- Write ocean header |
---|
| 1206 | WRITE( io, 1 ) |
---|
| 1207 | IF ( stokes_force ) WRITE( io, 2 ) stokes_waveheight, stokes_wavelength |
---|
| 1208 | IF ( wave_breaking ) THEN |
---|
| 1209 | WRITE( io, 3 ) alpha_wave_breaking, timescale_wave_breaking |
---|
| 1210 | ENDIF |
---|
[3303] | 1211 | IF ( .NOT. salinity ) WRITE( io, 4 ) |
---|
[3381] | 1212 | IF ( surface_cooling_spinup_time /= 999999.9_wp ) THEN |
---|
| 1213 | WRITE( io, 5 ) surface_cooling_spinup_time |
---|
| 1214 | ENDIF |
---|
[3302] | 1215 | |
---|
[4797] | 1216 | 1 FORMAT (//' Ocean settings:'/' ------------------------------------------'/) |
---|
| 1217 | 2 FORMAT (' --> Craik-Leibovich vortex force and Stokes drift switched',' on'/ & |
---|
[3302] | 1218 | ' waveheight: ',F4.1,' m wavelength: ',F6.1,' m') |
---|
[4797] | 1219 | 3 FORMAT (' --> wave breaking generated turbulence switched on'/ & |
---|
| 1220 | ' alpha: ',F4.1/' timescale:',F5.1,' s') |
---|
[3303] | 1221 | 4 FORMAT (' --> prognostic salinity equation is switched off' ) |
---|
[3381] | 1222 | 5 FORMAT (' --> surface heat flux is switched off after ',F8.1,' s') |
---|
[3302] | 1223 | |
---|
| 1224 | END SUBROUTINE ocean_header |
---|
| 1225 | |
---|
| 1226 | |
---|
[4797] | 1227 | !--------------------------------------------------------------------------------------------------! |
---|
[3302] | 1228 | ! Description: |
---|
| 1229 | ! ------------ |
---|
[3294] | 1230 | !> Allocate arrays and assign pointers. |
---|
[4797] | 1231 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1232 | SUBROUTINE ocean_init_arrays |
---|
| 1233 | |
---|
[4797] | 1234 | USE indices, & |
---|
[3873] | 1235 | ONLY: nxlg, nxrg, nyn, nyng, nys, nysg, nzb, nzt |
---|
[3294] | 1236 | |
---|
| 1237 | IMPLICIT NONE |
---|
| 1238 | |
---|
[4797] | 1239 | ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & |
---|
| 1240 | rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & |
---|
| 1241 | sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & |
---|
[3294] | 1242 | sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 1243 | |
---|
[4797] | 1244 | |
---|
[3303] | 1245 | IF ( salinity ) THEN |
---|
| 1246 | ALLOCATE( sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 1247 | ENDIF |
---|
| 1248 | |
---|
[3873] | 1249 | IF ( ws_scheme_sca ) THEN |
---|
| 1250 | ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) ) |
---|
| 1251 | sums_wssas_ws_l = 0.0_wp |
---|
| 1252 | ENDIF |
---|
| 1253 | |
---|
| 1254 | IF ( loop_optimization /= 'vector' ) THEN |
---|
| 1255 | |
---|
| 1256 | IF ( ws_scheme_sca ) THEN |
---|
| 1257 | ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1) ) |
---|
| 1258 | ALLOCATE( diss_s_sa(nzb+1:nzt,0:threads_per_task-1) ) |
---|
| 1259 | ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) |
---|
| 1260 | ALLOCATE( diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) |
---|
| 1261 | ENDIF |
---|
| 1262 | |
---|
| 1263 | ENDIF |
---|
| 1264 | |
---|
[3294] | 1265 | prho => prho_1 |
---|
[4797] | 1266 | rho_ocean => rho_1 ! routines calc_mean_profile and diffusion_e require density to be a pointer |
---|
[3294] | 1267 | |
---|
| 1268 | ! |
---|
| 1269 | !-- Initial assignment of pointers |
---|
[3303] | 1270 | IF ( salinity ) THEN |
---|
| 1271 | sa => sa_1; sa_p => sa_2; tsa_m => sa_3 |
---|
| 1272 | ELSE |
---|
| 1273 | sa => sa_1; sa_p => sa_1; tsa_m => sa_3 |
---|
| 1274 | ENDIF |
---|
[3294] | 1275 | |
---|
| 1276 | END SUBROUTINE ocean_init_arrays |
---|
| 1277 | |
---|
| 1278 | |
---|
[4797] | 1279 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1280 | ! Description: |
---|
| 1281 | ! ------------ |
---|
| 1282 | !> Initialization of quantities needed for the ocean mode |
---|
[4797] | 1283 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1284 | SUBROUTINE ocean_init |
---|
| 1285 | |
---|
| 1286 | |
---|
[4797] | 1287 | USE arrays_3d, & |
---|
| 1288 | ONLY: dzu, dzw, hyp, pt_init, ref_state, u_stokes_zu, u_stokes_zw, v_stokes_zu, & |
---|
| 1289 | v_stokes_zw, zu, zw |
---|
[3294] | 1290 | |
---|
[4797] | 1291 | USE basic_constants_and_equations_mod, & |
---|
[3294] | 1292 | ONLY: g |
---|
| 1293 | |
---|
[4797] | 1294 | USE basic_constants_and_equations_mod, & |
---|
[3302] | 1295 | ONLY: pi |
---|
| 1296 | |
---|
[4797] | 1297 | USE control_parameters, & |
---|
| 1298 | ONLY: initializing_actions, molecular_viscosity, rho_surface, rho_reference, & |
---|
| 1299 | surface_pressure, top_momentumflux_u, top_momentumflux_v, use_single_reference_value |
---|
[3294] | 1300 | |
---|
[4797] | 1301 | USE indices, & |
---|
[3294] | 1302 | ONLY: nxl, nxlg, nxrg, nyng, nys, nysg, nzb, nzt |
---|
| 1303 | |
---|
| 1304 | USE kinds |
---|
| 1305 | |
---|
[4797] | 1306 | USE pegrid, & |
---|
[3302] | 1307 | ONLY: myid |
---|
[3294] | 1308 | |
---|
[4797] | 1309 | USE statistics, & |
---|
[3294] | 1310 | ONLY: hom, statistic_regions |
---|
| 1311 | |
---|
| 1312 | IMPLICIT NONE |
---|
| 1313 | |
---|
| 1314 | INTEGER(iwp) :: i !< loop index |
---|
| 1315 | INTEGER(iwp) :: j !< loop index |
---|
| 1316 | INTEGER(iwp) :: k !< loop index |
---|
| 1317 | INTEGER(iwp) :: n !< loop index |
---|
| 1318 | |
---|
[4797] | 1319 | REAL(wp) :: alpha !< angle of surface stress |
---|
| 1320 | REAL(wp) :: dum !< dummy argument |
---|
| 1321 | REAL(wp) :: pt_l !< local scalar for pt used in equation of state function |
---|
| 1322 | REAL(wp) :: sa_l !< local scalar for sa used in equation of state function |
---|
[3302] | 1323 | REAL(wp) :: velocity_amplitude !< local scalar for amplitude of Stokes drift velocity |
---|
[4797] | 1324 | REAL(wp) :: x !< temporary variable to store surface stress along x |
---|
| 1325 | REAL(wp) :: y !< temporary variable to store surface stress along y |
---|
[3294] | 1326 | |
---|
| 1327 | REAL(wp), DIMENSION(nzb:nzt+1) :: rho_ocean_init !< local array for initial density |
---|
| 1328 | |
---|
| 1329 | ALLOCATE( hyp(nzb:nzt+1) ) |
---|
| 1330 | |
---|
| 1331 | |
---|
| 1332 | ! |
---|
[4797] | 1333 | !-- In case of no restart run, calculate the inital salinity profile using the given salinity |
---|
| 1334 | !-- gradients. |
---|
[3294] | 1335 | IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN |
---|
| 1336 | |
---|
| 1337 | sa_init = sa_surface |
---|
| 1338 | ! |
---|
[4797] | 1339 | !-- Last arguments gives back the gradient at top level to be used as possible Neumann boundary |
---|
| 1340 | !-- condition. This is not realized for the ocean mode, therefore a dummy argument is used. |
---|
[3303] | 1341 | IF ( salinity ) THEN |
---|
[4797] | 1342 | CALL init_vertical_profiles( sa_vertical_gradient_level_ind, & |
---|
| 1343 | sa_vertical_gradient_level, & |
---|
| 1344 | sa_vertical_gradient, sa_init, & |
---|
[3303] | 1345 | sa_surface, dum ) |
---|
| 1346 | ENDIF |
---|
[3294] | 1347 | ENDIF |
---|
| 1348 | |
---|
| 1349 | ! |
---|
| 1350 | !-- Initialize required 3d-arrays |
---|
[4797] | 1351 | IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & |
---|
[3294] | 1352 | TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN |
---|
| 1353 | ! |
---|
| 1354 | !-- Initialization via computed 1D-model profiles |
---|
| 1355 | IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 ) THEN |
---|
| 1356 | |
---|
| 1357 | DO i = nxlg, nxrg |
---|
| 1358 | DO j = nysg, nyng |
---|
| 1359 | sa(:,j,i) = sa_init |
---|
| 1360 | ENDDO |
---|
| 1361 | ENDDO |
---|
| 1362 | |
---|
| 1363 | ENDIF |
---|
| 1364 | ! |
---|
| 1365 | !-- Store initial profiles for output purposes etc. |
---|
| 1366 | !-- Store initial salinity profile |
---|
| 1367 | hom(:,1,26,:) = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 ) |
---|
| 1368 | ! |
---|
| 1369 | !-- Initialize old and new time levels. |
---|
| 1370 | tsa_m = 0.0_wp |
---|
| 1371 | sa_p = sa |
---|
| 1372 | |
---|
| 1373 | ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN |
---|
| 1374 | |
---|
| 1375 | ! |
---|
[4797] | 1376 | !-- Initialize new time levels (only done in order to set boundary values including ghost points) |
---|
[3294] | 1377 | sa_p = sa |
---|
| 1378 | ! |
---|
[4797] | 1379 | !-- Allthough tendency arrays are set in prognostic_equations, they have to be predefined here |
---|
| 1380 | !-- because they are used (but multiplied with 0) there before they are set. |
---|
[3294] | 1381 | tsa_m = 0.0_wp |
---|
| 1382 | |
---|
| 1383 | ENDIF |
---|
| 1384 | |
---|
| 1385 | ! |
---|
| 1386 | !-- Set water density near the ocean surface |
---|
| 1387 | rho_surface = 1027.62_wp |
---|
| 1388 | |
---|
| 1389 | ! |
---|
| 1390 | !-- Set kinematic viscosity to sea water at 20C. |
---|
| 1391 | !-- This changes the default value that is given for air! |
---|
| 1392 | molecular_viscosity = 1.05E-6_wp |
---|
| 1393 | |
---|
| 1394 | ! |
---|
[4797] | 1395 | !-- Change sign of buoyancy/stability terms because density gradient is used instead of the |
---|
| 1396 | !-- potential temperature gradient to calculate the buoyancy. |
---|
[3294] | 1397 | atmos_ocean_sign = -1.0_wp |
---|
| 1398 | |
---|
| 1399 | ! |
---|
[4797] | 1400 | !-- Calculate initial vertical profile of hydrostatic pressure (in Pa) and the reference density |
---|
| 1401 | !-- (used later in buoyancy term). |
---|
| 1402 | !-- First step: Calculate pressure using reference density. |
---|
[3294] | 1403 | hyp(nzt+1) = surface_pressure * 100.0_wp |
---|
| 1404 | hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1) |
---|
| 1405 | rho_ocean_init(nzt) = rho_surface |
---|
| 1406 | rho_ocean_init(nzt+1) = rho_surface ! only required for output |
---|
| 1407 | |
---|
| 1408 | DO k = nzt-1, 1, -1 |
---|
| 1409 | hyp(k) = hyp(k+1) + rho_surface * g * dzu(k) |
---|
| 1410 | ENDDO |
---|
| 1411 | hyp(0) = hyp(1) + rho_surface * g * dzu(1) |
---|
| 1412 | |
---|
| 1413 | ! |
---|
[4797] | 1414 | !-- Second step: Iteratively calculate in situ density (based on presssure) and pressure |
---|
| 1415 | !-- (based on in situ density) |
---|
[3294] | 1416 | DO n = 1, 5 |
---|
| 1417 | |
---|
| 1418 | rho_reference = rho_surface * 0.5_wp * dzu(nzt+1) |
---|
| 1419 | |
---|
| 1420 | DO k = nzt, 0, -1 |
---|
| 1421 | |
---|
| 1422 | sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) ) |
---|
| 1423 | pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) ) |
---|
| 1424 | |
---|
| 1425 | rho_ocean_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l ) |
---|
| 1426 | |
---|
| 1427 | rho_reference = rho_reference + rho_ocean_init(k) * dzu(k+1) |
---|
| 1428 | |
---|
| 1429 | ENDDO |
---|
| 1430 | |
---|
| 1431 | rho_reference = rho_reference / ( zw(nzt) - zu(nzb) ) |
---|
| 1432 | |
---|
| 1433 | hyp(nzt) = hyp(nzt+1) + rho_surface * g * 0.5_wp * dzu(nzt+1) |
---|
| 1434 | DO k = nzt-1, 0, -1 |
---|
[4797] | 1435 | hyp(k) = hyp(k+1) + g * 0.5_wp * ( rho_ocean_init(k) + rho_ocean_init(k+1) ) * dzu(k+1) |
---|
[3294] | 1436 | ENDDO |
---|
| 1437 | |
---|
| 1438 | ENDDO |
---|
| 1439 | |
---|
| 1440 | ! |
---|
| 1441 | !-- Calculate the reference potential density |
---|
| 1442 | prho_reference = 0.0_wp |
---|
| 1443 | DO k = 0, nzt |
---|
| 1444 | |
---|
| 1445 | sa_l = 0.5_wp * ( sa_init(k) + sa_init(k+1) ) |
---|
| 1446 | pt_l = 0.5_wp * ( pt_init(k) + pt_init(k+1) ) |
---|
| 1447 | |
---|
[4797] | 1448 | prho_reference = prho_reference + dzu(k+1) * eqn_state_seawater_func( 0.0_wp, pt_l, sa_l ) |
---|
[3294] | 1449 | |
---|
| 1450 | ENDDO |
---|
| 1451 | |
---|
| 1452 | prho_reference = prho_reference / ( zu(nzt) - zu(nzb) ) |
---|
| 1453 | |
---|
| 1454 | ! |
---|
[4797] | 1455 | !-- Calculate the 3d array of initial in situ and potential density, based on the initial |
---|
| 1456 | !-- temperature and salinity profile. |
---|
[3294] | 1457 | CALL eqn_state_seawater |
---|
| 1458 | |
---|
| 1459 | ! |
---|
| 1460 | !-- Store initial density profile |
---|
| 1461 | hom(:,1,77,:) = SPREAD( rho_ocean_init(:), 2, statistic_regions+1 ) |
---|
| 1462 | |
---|
| 1463 | ! |
---|
| 1464 | !-- Set the reference state to be used in the buoyancy terms |
---|
| 1465 | IF ( use_single_reference_value ) THEN |
---|
| 1466 | ref_state(:) = rho_reference |
---|
| 1467 | ELSE |
---|
| 1468 | ref_state(:) = rho_ocean_init(:) |
---|
| 1469 | ENDIF |
---|
| 1470 | |
---|
[3302] | 1471 | ! |
---|
| 1472 | !-- Calculate the Stokes drift velocity profile |
---|
| 1473 | IF ( stokes_force ) THEN |
---|
[3294] | 1474 | |
---|
[3302] | 1475 | ! |
---|
| 1476 | !-- First, calculate angle of surface stress |
---|
| 1477 | x = -top_momentumflux_u |
---|
| 1478 | y = -top_momentumflux_v |
---|
| 1479 | IF ( x == 0.0_wp ) THEN |
---|
| 1480 | IF ( y > 0.0_wp ) THEN |
---|
| 1481 | alpha = pi / 2.0_wp |
---|
| 1482 | ELSEIF ( y < 0.0_wp ) THEN |
---|
| 1483 | alpha = 3.0_wp * pi / 2.0_wp |
---|
| 1484 | ENDIF |
---|
| 1485 | ELSE |
---|
| 1486 | IF ( x < 0.0_wp ) THEN |
---|
| 1487 | alpha = ATAN( y / x ) + pi |
---|
| 1488 | ELSE |
---|
| 1489 | IF ( y < 0.0_wp ) THEN |
---|
| 1490 | alpha = ATAN( y / x ) + 2.0_wp * pi |
---|
| 1491 | ELSE |
---|
| 1492 | alpha = ATAN( y / x ) |
---|
| 1493 | ENDIF |
---|
| 1494 | ENDIF |
---|
| 1495 | ENDIF |
---|
| 1496 | |
---|
[4797] | 1497 | velocity_amplitude = ( pi * stokes_waveheight / stokes_wavelength )**2 * & |
---|
[3302] | 1498 | SQRT( g * stokes_wavelength / ( 2.0_wp * pi ) ) |
---|
| 1499 | |
---|
| 1500 | DO k = nzb, nzt |
---|
[4797] | 1501 | u_stokes_zu(k) = velocity_amplitude * COS( alpha ) * & |
---|
[3302] | 1502 | EXP( 4.0_wp * pi * zu(k) / stokes_wavelength ) |
---|
[4797] | 1503 | u_stokes_zw(k) = velocity_amplitude * COS( alpha ) * & |
---|
[3302] | 1504 | EXP( 4.0_wp * pi * zw(k) / stokes_wavelength ) |
---|
[4797] | 1505 | v_stokes_zu(k) = velocity_amplitude * SIN( alpha ) * & |
---|
[3302] | 1506 | EXP( 4.0_wp * pi * zu(k) / stokes_wavelength ) |
---|
| 1507 | v_stokes_zw(k) = velocity_amplitude * SIN( alpha ) * & |
---|
| 1508 | EXP( 4.0_wp * pi * zw(k) / stokes_wavelength ) |
---|
| 1509 | ENDDO |
---|
| 1510 | u_stokes_zu(nzt+1) = u_stokes_zw(nzt) ! because zu(nzt+1) changes the sign |
---|
| 1511 | u_stokes_zw(nzt+1) = u_stokes_zw(nzt) ! because zw(nzt+1) changes the sign |
---|
| 1512 | v_stokes_zu(nzt+1) = v_stokes_zw(nzt) ! because zu(nzt+1) changes the sign |
---|
| 1513 | v_stokes_zw(nzt+1) = v_stokes_zw(nzt) ! because zw(nzt+1) changes the sign |
---|
| 1514 | |
---|
| 1515 | ENDIF |
---|
| 1516 | |
---|
| 1517 | ! |
---|
| 1518 | !-- Wave breaking effects |
---|
| 1519 | IF ( wave_breaking ) THEN |
---|
| 1520 | ! |
---|
| 1521 | !-- Calculate friction velocity at ocean surface |
---|
[4797] | 1522 | u_star_wave_breaking = SQRT( SQRT( top_momentumflux_u**2 + top_momentumflux_v**2 ) ) |
---|
[3302] | 1523 | ! |
---|
[4797] | 1524 | !-- Set the time scale of random forcing. The vertical grid spacing at the ocean surface is |
---|
| 1525 | !-- assumed as the length scale of turbulence. |
---|
[3302] | 1526 | !-- Formula follows Noh et al. (2004), JPO |
---|
[4797] | 1527 | timescale_wave_breaking = 0.1_wp * dzw(nzt) / alpha_wave_breaking / u_star_wave_breaking |
---|
[3302] | 1528 | ! |
---|
[4797] | 1529 | !-- Set random number seeds differently on the processor cores in order to create different |
---|
| 1530 | !-- random number sequences |
---|
[3302] | 1531 | iran_ocean = iran_ocean + myid |
---|
| 1532 | ENDIF |
---|
| 1533 | |
---|
[3294] | 1534 | END SUBROUTINE ocean_init |
---|
| 1535 | |
---|
| 1536 | |
---|
[4797] | 1537 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1538 | ! Description: |
---|
| 1539 | ! ------------ |
---|
[3873] | 1540 | !> Call for all grid points |
---|
[4797] | 1541 | !--------------------------------------------------------------------------------------------------! |
---|
[3873] | 1542 | SUBROUTINE ocean_actions( location ) |
---|
| 1543 | |
---|
| 1544 | |
---|
| 1545 | CHARACTER (LEN=*), INTENT(IN) :: location !< call location string |
---|
| 1546 | |
---|
| 1547 | SELECT CASE ( location ) |
---|
| 1548 | |
---|
| 1549 | CASE ( 'before_timestep' ) |
---|
| 1550 | |
---|
| 1551 | IF ( ws_scheme_sca ) THEN |
---|
| 1552 | sums_wssas_ws_l = 0.0_wp |
---|
| 1553 | ENDIF |
---|
| 1554 | |
---|
| 1555 | CASE DEFAULT |
---|
| 1556 | CONTINUE |
---|
| 1557 | |
---|
| 1558 | END SELECT |
---|
| 1559 | |
---|
| 1560 | END SUBROUTINE ocean_actions |
---|
| 1561 | |
---|
| 1562 | |
---|
[4797] | 1563 | !--------------------------------------------------------------------------------------------------! |
---|
[3873] | 1564 | ! Description: |
---|
| 1565 | ! ------------ |
---|
| 1566 | !> Call for grid points i,j |
---|
[4797] | 1567 | !--------------------------------------------------------------------------------------------------! |
---|
[3873] | 1568 | SUBROUTINE ocean_actions_ij( i, j, location ) |
---|
| 1569 | |
---|
| 1570 | CHARACTER (LEN=*), INTENT(IN) :: location !< call location string |
---|
| 1571 | |
---|
[4797] | 1572 | INTEGER(iwp) :: dummy !< call location string |
---|
| 1573 | INTEGER(iwp), INTENT(IN) :: i !< grid index in x-direction |
---|
| 1574 | INTEGER(iwp), INTENT(IN) :: j !< grid index in y-direction |
---|
| 1575 | |
---|
| 1576 | |
---|
[3873] | 1577 | IF ( ocean_mode ) dummy = i + j |
---|
| 1578 | |
---|
| 1579 | SELECT CASE ( location ) |
---|
| 1580 | |
---|
| 1581 | CASE ( 'before_timestep' ) |
---|
| 1582 | |
---|
| 1583 | IF ( ws_scheme_sca ) THEN |
---|
| 1584 | sums_wssas_ws_l = 0.0_wp |
---|
| 1585 | ENDIF |
---|
| 1586 | |
---|
| 1587 | CASE DEFAULT |
---|
| 1588 | CONTINUE |
---|
| 1589 | |
---|
| 1590 | END SELECT |
---|
| 1591 | |
---|
| 1592 | END SUBROUTINE ocean_actions_ij |
---|
| 1593 | |
---|
| 1594 | |
---|
[4797] | 1595 | !--------------------------------------------------------------------------------------------------! |
---|
[3873] | 1596 | ! Description: |
---|
| 1597 | ! ------------ |
---|
[3294] | 1598 | !> Prognostic equation for salinity. |
---|
| 1599 | !> Vector-optimized version |
---|
[4797] | 1600 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1601 | SUBROUTINE ocean_prognostic_equations |
---|
| 1602 | |
---|
[4797] | 1603 | USE advec_s_bc_mod, & |
---|
[3294] | 1604 | ONLY: advec_s_bc |
---|
| 1605 | |
---|
[4797] | 1606 | USE advec_s_pw_mod, & |
---|
[3294] | 1607 | ONLY: advec_s_pw |
---|
| 1608 | |
---|
[4797] | 1609 | USE advec_s_up_mod, & |
---|
[3294] | 1610 | ONLY: advec_s_up |
---|
| 1611 | |
---|
[4797] | 1612 | USE advec_ws, & |
---|
[3294] | 1613 | ONLY: advec_s_ws |
---|
| 1614 | |
---|
[4797] | 1615 | USE arrays_3d, & |
---|
[3294] | 1616 | ONLY: rdf_sc, tend, tsa_m |
---|
| 1617 | |
---|
[4797] | 1618 | USE control_parameters, & |
---|
| 1619 | ONLY: bc_dirichlet_l, & |
---|
| 1620 | bc_dirichlet_n, & |
---|
| 1621 | bc_dirichlet_r, & |
---|
| 1622 | bc_dirichlet_s, & |
---|
| 1623 | bc_radiation_l, & |
---|
| 1624 | bc_radiation_n, & |
---|
| 1625 | bc_radiation_r, & |
---|
| 1626 | bc_radiation_s, & |
---|
| 1627 | dt_3d, & |
---|
| 1628 | intermediate_timestep_count, & |
---|
| 1629 | intermediate_timestep_count_max, & |
---|
| 1630 | scalar_advec, & |
---|
| 1631 | simulated_time, & |
---|
| 1632 | timestep_scheme, & |
---|
| 1633 | tsc, & |
---|
| 1634 | ws_scheme_sca |
---|
[3294] | 1635 | |
---|
[4797] | 1636 | USE cpulog, & |
---|
[3719] | 1637 | ONLY: cpu_log, log_point_s |
---|
[3294] | 1638 | |
---|
[4797] | 1639 | USE diffusion_s_mod, & |
---|
[3294] | 1640 | ONLY: diffusion_s |
---|
| 1641 | |
---|
| 1642 | IMPLICIT NONE |
---|
| 1643 | |
---|
| 1644 | INTEGER(iwp) :: i !< loop index |
---|
| 1645 | INTEGER(iwp) :: j !< loop index |
---|
| 1646 | INTEGER(iwp) :: k !< loop index |
---|
| 1647 | |
---|
| 1648 | REAL(wp) :: sbt !< weighting factor for sub-time step |
---|
| 1649 | |
---|
[4797] | 1650 | |
---|
[3294] | 1651 | ! |
---|
[3381] | 1652 | !-- Switch of the surface heat flux, if requested |
---|
| 1653 | IF ( surface_cooling_spinup_time /= 999999.9_wp ) THEN |
---|
[4797] | 1654 | IF ( .NOT. surface_cooling_switched_off .AND. & |
---|
[3381] | 1655 | simulated_time >= surface_cooling_spinup_time ) THEN |
---|
| 1656 | |
---|
| 1657 | surf_def_h(2)%shf = 0.0_wp |
---|
| 1658 | surface_cooling_switched_off = .TRUE. |
---|
| 1659 | |
---|
| 1660 | ENDIF |
---|
| 1661 | ENDIF |
---|
| 1662 | |
---|
| 1663 | ! |
---|
[4797] | 1664 | !-- Compute prognostic equations for the ocean mode. |
---|
| 1665 | !-- First, start with salinity. |
---|
[3568] | 1666 | IF ( salinity ) THEN |
---|
[3303] | 1667 | |
---|
[3719] | 1668 | CALL cpu_log( log_point_s(20), 'sa-equation', 'start' ) |
---|
[3294] | 1669 | |
---|
| 1670 | ! |
---|
[3568] | 1671 | !-- sa-tendency terms with communication |
---|
| 1672 | sbt = tsc(2) |
---|
| 1673 | IF ( scalar_advec == 'bc-scheme' ) THEN |
---|
[3294] | 1674 | |
---|
[3568] | 1675 | IF ( timestep_scheme(1:5) /= 'runge' ) THEN |
---|
[3294] | 1676 | ! |
---|
[3568] | 1677 | !-- Bott-Chlond scheme always uses Euler time step. Thus: |
---|
| 1678 | sbt = 1.0_wp |
---|
| 1679 | ENDIF |
---|
| 1680 | tend = 0.0_wp |
---|
| 1681 | CALL advec_s_bc( sa, 'sa' ) |
---|
| 1682 | |
---|
[3294] | 1683 | ENDIF |
---|
| 1684 | |
---|
| 1685 | ! |
---|
[3568] | 1686 | !-- sa-tendency terms with no communication |
---|
| 1687 | IF ( scalar_advec /= 'bc-scheme' ) THEN |
---|
| 1688 | tend = 0.0_wp |
---|
| 1689 | IF ( timestep_scheme(1:5) == 'runge' ) THEN |
---|
| 1690 | IF ( ws_scheme_sca ) THEN |
---|
[4797] | 1691 | CALL advec_s_ws( advc_flags_s, sa, 'sa', & |
---|
| 1692 | bc_dirichlet_l .OR. bc_radiation_l, & |
---|
| 1693 | bc_dirichlet_n .OR. bc_radiation_n, & |
---|
| 1694 | bc_dirichlet_r .OR. bc_radiation_r, & |
---|
[4109] | 1695 | bc_dirichlet_s .OR. bc_radiation_s ) |
---|
[3568] | 1696 | ELSE |
---|
| 1697 | CALL advec_s_pw( sa ) |
---|
| 1698 | ENDIF |
---|
[3294] | 1699 | ELSE |
---|
[3568] | 1700 | CALL advec_s_up( sa ) |
---|
[3294] | 1701 | ENDIF |
---|
| 1702 | ENDIF |
---|
| 1703 | |
---|
[4797] | 1704 | CALL diffusion_s( sa, & |
---|
| 1705 | surf_def_h(0)%sasws, surf_def_h(1)%sasws, & |
---|
| 1706 | surf_def_h(2)%sasws, & |
---|
| 1707 | surf_lsm_h(0)%sasws, surf_lsm_h(1)%sasws, & |
---|
| 1708 | surf_usm_h(0)%sasws, surf_usm_h(1)%sasws, & |
---|
| 1709 | surf_def_v(0)%sasws, surf_def_v(1)%sasws, & |
---|
| 1710 | surf_def_v(2)%sasws, surf_def_v(3)%sasws, & |
---|
| 1711 | surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, & |
---|
| 1712 | surf_lsm_v(2)%sasws, surf_lsm_v(3)%sasws, & |
---|
| 1713 | surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, & |
---|
[3568] | 1714 | surf_usm_v(2)%sasws, surf_usm_v(3)%sasws ) |
---|
[3294] | 1715 | |
---|
[3684] | 1716 | ! CALL user_actions( 'sa-tendency' ) ToDo: find general solution for dependency between modules |
---|
[3294] | 1717 | |
---|
| 1718 | ! |
---|
[3568] | 1719 | !-- Prognostic equation for salinity |
---|
| 1720 | DO i = nxl, nxr |
---|
| 1721 | DO j = nys, nyn |
---|
[4370] | 1722 | !following directive is required to vectorize on Intel19 |
---|
| 1723 | !DIR$ IVDEP |
---|
[3568] | 1724 | DO k = nzb+1, nzt |
---|
[4797] | 1725 | sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) ) & |
---|
| 1726 | - tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) ) & |
---|
| 1727 | ) & |
---|
| 1728 | * MERGE( 1.0_wp, 0.0_wp, & |
---|
| 1729 | BTEST( wall_flags_total_0(k,j,i), 0 ) & |
---|
[3568] | 1730 | ) |
---|
| 1731 | IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i) |
---|
| 1732 | ENDDO |
---|
[3294] | 1733 | ENDDO |
---|
| 1734 | ENDDO |
---|
| 1735 | |
---|
| 1736 | ! |
---|
[3568] | 1737 | !-- Calculate tendencies for the next Runge-Kutta step |
---|
| 1738 | IF ( timestep_scheme(1:5) == 'runge' ) THEN |
---|
| 1739 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
| 1740 | DO i = nxl, nxr |
---|
| 1741 | DO j = nys, nyn |
---|
| 1742 | DO k = nzb+1, nzt |
---|
| 1743 | tsa_m(k,j,i) = tend(k,j,i) |
---|
| 1744 | ENDDO |
---|
[3294] | 1745 | ENDDO |
---|
| 1746 | ENDDO |
---|
[4797] | 1747 | ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN |
---|
[3568] | 1748 | DO i = nxl, nxr |
---|
| 1749 | DO j = nys, nyn |
---|
| 1750 | DO k = nzb+1, nzt |
---|
[4797] | 1751 | tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tsa_m(k,j,i) |
---|
[3568] | 1752 | ENDDO |
---|
[3294] | 1753 | ENDDO |
---|
| 1754 | ENDDO |
---|
[3568] | 1755 | ENDIF |
---|
[3294] | 1756 | ENDIF |
---|
[3568] | 1757 | |
---|
[3719] | 1758 | CALL cpu_log( log_point_s(20), 'sa-equation', 'stop' ) |
---|
[3568] | 1759 | |
---|
[3294] | 1760 | ENDIF |
---|
| 1761 | |
---|
| 1762 | ! |
---|
| 1763 | !-- Calculate density by the equation of state for seawater |
---|
[3719] | 1764 | CALL cpu_log( log_point_s(21), 'eqns-seawater', 'start' ) |
---|
[3294] | 1765 | CALL eqn_state_seawater |
---|
[3719] | 1766 | CALL cpu_log( log_point_s(21), 'eqns-seawater', 'stop' ) |
---|
[3294] | 1767 | |
---|
| 1768 | END SUBROUTINE ocean_prognostic_equations |
---|
| 1769 | |
---|
| 1770 | |
---|
[4797] | 1771 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1772 | ! Description: |
---|
| 1773 | ! ------------ |
---|
| 1774 | !> Prognostic equations for ocean mode (so far, salinity only) |
---|
| 1775 | !> Cache-optimized version |
---|
[4797] | 1776 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1777 | SUBROUTINE ocean_prognostic_equations_ij( i, j, i_omp_start, tn ) |
---|
| 1778 | |
---|
[4797] | 1779 | USE advec_s_pw_mod, & |
---|
[3294] | 1780 | ONLY: advec_s_pw |
---|
| 1781 | |
---|
[4797] | 1782 | USE advec_s_up_mod, & |
---|
[3294] | 1783 | ONLY: advec_s_up |
---|
| 1784 | |
---|
[4797] | 1785 | USE advec_ws, & |
---|
[3294] | 1786 | ONLY: advec_s_ws |
---|
| 1787 | |
---|
[4797] | 1788 | USE arrays_3d, & |
---|
[3294] | 1789 | ONLY: diss_l_sa, diss_s_sa, flux_l_sa, flux_s_sa, rdf_sc, tend, tsa_m |
---|
| 1790 | |
---|
[4797] | 1791 | USE control_parameters, & |
---|
| 1792 | ONLY: bc_dirichlet_l, & |
---|
| 1793 | bc_dirichlet_n, & |
---|
| 1794 | bc_dirichlet_r, & |
---|
| 1795 | bc_dirichlet_s, & |
---|
| 1796 | bc_radiation_l, & |
---|
| 1797 | bc_radiation_n, & |
---|
| 1798 | bc_radiation_r, & |
---|
| 1799 | bc_radiation_s, & |
---|
| 1800 | dt_3d, & |
---|
| 1801 | intermediate_timestep_count, & |
---|
| 1802 | intermediate_timestep_count_max, & |
---|
| 1803 | simulated_time, & |
---|
| 1804 | timestep_scheme, & |
---|
| 1805 | tsc, & |
---|
| 1806 | ws_scheme_sca |
---|
[3294] | 1807 | |
---|
[4797] | 1808 | USE diffusion_s_mod, & |
---|
[3294] | 1809 | ONLY: diffusion_s |
---|
| 1810 | |
---|
| 1811 | IMPLICIT NONE |
---|
| 1812 | |
---|
| 1813 | INTEGER(iwp) :: i !< loop index x direction |
---|
[4797] | 1814 | INTEGER(iwp) :: i_omp_start !< first loop index of i-loop in calling routine prognostic_equations |
---|
[3294] | 1815 | INTEGER(iwp) :: j !< loop index y direction |
---|
| 1816 | INTEGER(iwp) :: k !< loop index z direction |
---|
| 1817 | INTEGER(iwp) :: tn !< task number of openmp task |
---|
| 1818 | |
---|
[3381] | 1819 | |
---|
[3294] | 1820 | ! |
---|
[3381] | 1821 | !-- Switch of the surface heat flux, if requested |
---|
| 1822 | IF ( surface_cooling_spinup_time /= 999999.9_wp ) THEN |
---|
[4797] | 1823 | IF ( .NOT. surface_cooling_switched_off .AND. & |
---|
[3381] | 1824 | simulated_time >= surface_cooling_spinup_time ) THEN |
---|
| 1825 | |
---|
| 1826 | surf_def_h(2)%shf = 0.0_wp |
---|
| 1827 | surface_cooling_switched_off = .TRUE. |
---|
| 1828 | |
---|
| 1829 | ENDIF |
---|
| 1830 | ENDIF |
---|
| 1831 | |
---|
| 1832 | ! |
---|
[4797] | 1833 | !-- Compute prognostic equations for the ocean mode. |
---|
| 1834 | !-- First, start with tendency-terms for salinity. |
---|
[3568] | 1835 | IF ( salinity ) THEN |
---|
[3303] | 1836 | |
---|
[3568] | 1837 | tend(:,j,i) = 0.0_wp |
---|
[4797] | 1838 | IF ( timestep_scheme(1:5) == 'runge' ) THEN |
---|
[3568] | 1839 | IF ( ws_scheme_sca ) THEN |
---|
[4797] | 1840 | CALL advec_s_ws( advc_flags_s, & |
---|
| 1841 | i, j, sa, 'sa', flux_s_sa, diss_s_sa, flux_l_sa, diss_l_sa, & |
---|
| 1842 | i_omp_start, tn, & |
---|
| 1843 | bc_dirichlet_l .OR. bc_radiation_l, & |
---|
| 1844 | bc_dirichlet_n .OR. bc_radiation_n, & |
---|
| 1845 | bc_dirichlet_r .OR. bc_radiation_r, & |
---|
[4109] | 1846 | bc_dirichlet_s .OR. bc_radiation_s ) |
---|
[3568] | 1847 | ELSE |
---|
| 1848 | CALL advec_s_pw( i, j, sa ) |
---|
| 1849 | ENDIF |
---|
[3294] | 1850 | ELSE |
---|
[3568] | 1851 | CALL advec_s_up( i, j, sa ) |
---|
[3294] | 1852 | ENDIF |
---|
[4797] | 1853 | CALL diffusion_s( i, j, sa, & |
---|
| 1854 | surf_def_h(0)%sasws, surf_def_h(1)%sasws, surf_def_h(2)%sasws, & |
---|
| 1855 | surf_lsm_h(0)%sasws, surf_lsm_h(1)%sasws, & |
---|
| 1856 | surf_usm_h(0)%sasws, surf_usm_h(1)%sasws, & |
---|
| 1857 | surf_def_v(0)%sasws, surf_def_v(1)%sasws, surf_def_v(2)%sasws, & |
---|
| 1858 | surf_def_v(3)%sasws, & |
---|
| 1859 | surf_lsm_v(0)%sasws, surf_lsm_v(1)%sasws, surf_lsm_v(2)%sasws, & |
---|
| 1860 | surf_lsm_v(3)%sasws, & |
---|
| 1861 | surf_usm_v(0)%sasws, surf_usm_v(1)%sasws, surf_usm_v(2)%sasws, & |
---|
| 1862 | surf_usm_v(3)%sasws ) |
---|
[3294] | 1863 | |
---|
[3684] | 1864 | ! CALL user_actions( i, j, 'sa-tendency' ) ToDo: find general solution for dependency between modules |
---|
[3294] | 1865 | |
---|
| 1866 | ! |
---|
[3568] | 1867 | !-- Prognostic equation for salinity |
---|
| 1868 | DO k = nzb+1, nzt |
---|
[3294] | 1869 | |
---|
[4797] | 1870 | sa_p(k,j,i) = sa(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) ) & |
---|
| 1871 | - tsc(5) * rdf_sc(k) * ( sa(k,j,i) - sa_init(k) ) & |
---|
| 1872 | ) * MERGE( 1.0_wp, 0.0_wp, & |
---|
| 1873 | BTEST( wall_flags_total_0(k,j,i), 0 ) & |
---|
| 1874 | ) |
---|
[3294] | 1875 | |
---|
[3568] | 1876 | IF ( sa_p(k,j,i) < 0.0_wp ) sa_p(k,j,i) = 0.1_wp * sa(k,j,i) |
---|
[3294] | 1877 | |
---|
[3568] | 1878 | ENDDO |
---|
[3294] | 1879 | |
---|
| 1880 | ! |
---|
[3568] | 1881 | !-- Calculate tendencies for the next Runge-Kutta step |
---|
| 1882 | IF ( timestep_scheme(1:5) == 'runge' ) THEN |
---|
| 1883 | IF ( intermediate_timestep_count == 1 ) THEN |
---|
| 1884 | DO k = nzb+1, nzt |
---|
| 1885 | tsa_m(k,j,i) = tend(k,j,i) |
---|
| 1886 | ENDDO |
---|
[4797] | 1887 | ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max ) THEN |
---|
[3568] | 1888 | DO k = nzb+1, nzt |
---|
[4797] | 1889 | tsa_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tsa_m(k,j,i) |
---|
[3568] | 1890 | ENDDO |
---|
| 1891 | ENDIF |
---|
[3294] | 1892 | ENDIF |
---|
[3568] | 1893 | |
---|
[3294] | 1894 | ENDIF |
---|
| 1895 | |
---|
| 1896 | ! |
---|
| 1897 | !-- Calculate density by the equation of state for seawater |
---|
| 1898 | CALL eqn_state_seawater( i, j ) |
---|
| 1899 | |
---|
| 1900 | END SUBROUTINE ocean_prognostic_equations_ij |
---|
| 1901 | |
---|
[4797] | 1902 | !--------------------------------------------------------------------------------------------------! |
---|
[4272] | 1903 | ! Description: |
---|
| 1904 | ! ------------ |
---|
| 1905 | !> Boundary conditions for ocean model |
---|
[4797] | 1906 | !--------------------------------------------------------------------------------------------------! |
---|
[4272] | 1907 | SUBROUTINE ocean_boundary_conditions |
---|
[3294] | 1908 | |
---|
[4272] | 1909 | IMPLICIT NONE |
---|
| 1910 | |
---|
| 1911 | INTEGER(iwp) :: i !< grid index x direction. |
---|
| 1912 | INTEGER(iwp) :: j !< grid index y direction. |
---|
| 1913 | INTEGER(iwp) :: k !< grid index z direction. |
---|
| 1914 | INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls. |
---|
| 1915 | INTEGER(iwp) :: m !< running index surface elements. |
---|
| 1916 | |
---|
| 1917 | ! |
---|
[4797] | 1918 | !-- Boundary conditions for salinity. |
---|
| 1919 | !-- Bottom boundary: Neumann condition because salinity flux is always given. |
---|
[4272] | 1920 | DO l = 0, 1 |
---|
| 1921 | !$OMP PARALLEL DO PRIVATE( i, j, k ) |
---|
| 1922 | DO m = 1, bc_h(l)%ns |
---|
| 1923 | i = bc_h(l)%i(m) |
---|
| 1924 | j = bc_h(l)%j(m) |
---|
| 1925 | k = bc_h(l)%k(m) |
---|
| 1926 | sa_p(k+bc_h(l)%koff,j,i) = sa_p(k,j,i) |
---|
| 1927 | ENDDO |
---|
| 1928 | ENDDO |
---|
| 1929 | ! |
---|
| 1930 | !-- Top boundary: Dirichlet or Neumann |
---|
| 1931 | IF ( ibc_sa_t == 0 ) THEN |
---|
| 1932 | sa_p(nzt+1,:,:) = sa(nzt+1,:,:) |
---|
| 1933 | ELSEIF ( ibc_sa_t == 1 ) THEN |
---|
| 1934 | sa_p(nzt+1,:,:) = sa_p(nzt,:,:) |
---|
| 1935 | ENDIF |
---|
| 1936 | |
---|
| 1937 | END SUBROUTINE ocean_boundary_conditions |
---|
| 1938 | |
---|
[4797] | 1939 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1940 | ! Description: |
---|
| 1941 | ! ------------ |
---|
| 1942 | !> Swapping of timelevels. |
---|
[4797] | 1943 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1944 | SUBROUTINE ocean_swap_timelevel( mod_count ) |
---|
| 1945 | |
---|
| 1946 | IMPLICIT NONE |
---|
| 1947 | |
---|
| 1948 | INTEGER, INTENT(IN) :: mod_count !< flag defining where pointers point to |
---|
| 1949 | |
---|
| 1950 | |
---|
| 1951 | SELECT CASE ( mod_count ) |
---|
| 1952 | |
---|
| 1953 | CASE ( 0 ) |
---|
[3303] | 1954 | IF ( salinity ) THEN |
---|
| 1955 | sa => sa_1; sa_p => sa_2 |
---|
| 1956 | ENDIF |
---|
[3294] | 1957 | |
---|
| 1958 | CASE ( 1 ) |
---|
[3303] | 1959 | IF ( salinity ) THEN |
---|
| 1960 | sa => sa_2; sa_p => sa_1 |
---|
| 1961 | ENDIF |
---|
[3294] | 1962 | |
---|
| 1963 | END SELECT |
---|
| 1964 | |
---|
| 1965 | END SUBROUTINE ocean_swap_timelevel |
---|
| 1966 | |
---|
| 1967 | |
---|
[4797] | 1968 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 1969 | ! Description: |
---|
| 1970 | ! ------------ |
---|
[4495] | 1971 | !> Read module-specific global restart data (Fortran binary format). |
---|
[4797] | 1972 | !--------------------------------------------------------------------------------------------------! |
---|
[4495] | 1973 | SUBROUTINE ocean_rrd_global_ftn( found ) |
---|
[3294] | 1974 | |
---|
| 1975 | |
---|
[4797] | 1976 | USE control_parameters, & |
---|
[3294] | 1977 | ONLY: length, restart_string |
---|
| 1978 | |
---|
| 1979 | |
---|
| 1980 | IMPLICIT NONE |
---|
| 1981 | |
---|
| 1982 | LOGICAL, INTENT(OUT) :: found |
---|
| 1983 | |
---|
| 1984 | |
---|
| 1985 | found = .TRUE. |
---|
| 1986 | |
---|
| 1987 | SELECT CASE ( restart_string(1:length) ) |
---|
| 1988 | |
---|
| 1989 | CASE ( 'bc_sa_t' ) |
---|
| 1990 | READ ( 13 ) bc_sa_t |
---|
| 1991 | |
---|
| 1992 | CASE ( 'bottom_salinityflux' ) |
---|
| 1993 | READ ( 13 ) bottom_salinityflux |
---|
| 1994 | |
---|
[3303] | 1995 | CASE ( 'salinity' ) |
---|
| 1996 | READ ( 13 ) salinity |
---|
| 1997 | |
---|
[3294] | 1998 | CASE ( 'sa_init' ) |
---|
| 1999 | READ ( 13 ) sa_init |
---|
| 2000 | |
---|
| 2001 | CASE ( 'sa_surface' ) |
---|
| 2002 | READ ( 13 ) sa_surface |
---|
| 2003 | |
---|
| 2004 | CASE ( 'sa_vertical_gradient' ) |
---|
| 2005 | READ ( 13 ) sa_vertical_gradient |
---|
| 2006 | |
---|
| 2007 | CASE ( 'sa_vertical_gradient_level' ) |
---|
| 2008 | READ ( 13 ) sa_vertical_gradient_level |
---|
| 2009 | |
---|
| 2010 | CASE ( 'stokes_waveheight' ) |
---|
| 2011 | READ ( 13 ) stokes_waveheight |
---|
| 2012 | |
---|
| 2013 | CASE ( 'stokes_wavelength' ) |
---|
| 2014 | READ ( 13 ) stokes_wavelength |
---|
| 2015 | |
---|
[3381] | 2016 | CASE ( 'surface_cooling_spinup_time' ) |
---|
| 2017 | READ ( 13 ) surface_cooling_spinup_time |
---|
| 2018 | |
---|
[3294] | 2019 | CASE ( 'top_salinityflux' ) |
---|
| 2020 | READ ( 13 ) top_salinityflux |
---|
| 2021 | |
---|
| 2022 | CASE ( 'wall_salinityflux' ) |
---|
| 2023 | READ ( 13 ) wall_salinityflux |
---|
| 2024 | |
---|
[3302] | 2025 | CASE ( 'wave_breaking' ) |
---|
| 2026 | READ ( 13 ) wave_breaking |
---|
| 2027 | |
---|
[3294] | 2028 | CASE DEFAULT |
---|
| 2029 | |
---|
| 2030 | found = .FALSE. |
---|
| 2031 | |
---|
| 2032 | END SELECT |
---|
| 2033 | |
---|
[4495] | 2034 | END SUBROUTINE ocean_rrd_global_ftn |
---|
[3294] | 2035 | |
---|
| 2036 | |
---|
[4797] | 2037 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 2038 | ! Description: |
---|
| 2039 | ! ------------ |
---|
[4495] | 2040 | !> Read module-specific global restart data (MPI-IO). |
---|
[4797] | 2041 | !--------------------------------------------------------------------------------------------------! |
---|
[4495] | 2042 | SUBROUTINE ocean_rrd_global_mpi |
---|
| 2043 | |
---|
| 2044 | CALL rrd_mpi_io( 'bc_sa_t', bc_sa_t ) |
---|
| 2045 | CALL rrd_mpi_io( 'bottom_salinityflux', bottom_salinityflux ) |
---|
| 2046 | CALL rrd_mpi_io( 'salinity', salinity ) |
---|
| 2047 | CALL rrd_mpi_io_global_array( 'sa_init', sa_init ) |
---|
| 2048 | CALL rrd_mpi_io( 'sa_surface', sa_surface ) |
---|
| 2049 | CALL rrd_mpi_io_global_array( 'sa_vertical_gradient', sa_vertical_gradient ) |
---|
| 2050 | CALL rrd_mpi_io_global_array( 'sa_vertical_gradient_level', sa_vertical_gradient_level ) |
---|
| 2051 | CALL rrd_mpi_io_global_array( 'sa_vertical_gradient_level_ind', sa_vertical_gradient_level_ind ) |
---|
| 2052 | CALL rrd_mpi_io( 'stokes_waveheight', stokes_waveheight ) |
---|
| 2053 | CALL rrd_mpi_io( 'stokes_wavelength', stokes_wavelength ) |
---|
| 2054 | CALL rrd_mpi_io( 'surface_cooling_spinup_time', surface_cooling_spinup_time ) |
---|
| 2055 | CALL rrd_mpi_io( 'top_salinityflux', top_salinityflux ) |
---|
| 2056 | CALL rrd_mpi_io_global_array( 'wall_salinityflux', wall_salinityflux ) |
---|
| 2057 | CALL rrd_mpi_io( 'wave_breaking', wave_breaking ) |
---|
| 2058 | |
---|
| 2059 | END SUBROUTINE ocean_rrd_global_mpi |
---|
| 2060 | |
---|
| 2061 | |
---|
[4797] | 2062 | !--------------------------------------------------------------------------------------------------! |
---|
[4495] | 2063 | ! Description: |
---|
| 2064 | ! ------------ |
---|
[4517] | 2065 | !> Read module-specific local restart data arrays (Fortran binary format). |
---|
[4797] | 2066 | !--------------------------------------------------------------------------------------------------! |
---|
| 2067 | SUBROUTINE ocean_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync, & |
---|
| 2068 | nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found ) |
---|
[3294] | 2069 | |
---|
[4797] | 2070 | USE averaging, & |
---|
[3294] | 2071 | ONLY: rho_ocean_av, sa_av |
---|
| 2072 | |
---|
[4797] | 2073 | USE control_parameters, & |
---|
[3294] | 2074 | ONLY: length, restart_string |
---|
| 2075 | |
---|
[4797] | 2076 | USE indices, & |
---|
[3294] | 2077 | ONLY: nbgp, nxlg, nxrg, nyng, nysg, nzb, nzt |
---|
| 2078 | |
---|
| 2079 | USE pegrid |
---|
| 2080 | |
---|
| 2081 | |
---|
| 2082 | IMPLICIT NONE |
---|
| 2083 | |
---|
| 2084 | INTEGER(iwp) :: k !< |
---|
| 2085 | INTEGER(iwp) :: nxlc !< |
---|
| 2086 | INTEGER(iwp) :: nxlf !< |
---|
| 2087 | INTEGER(iwp) :: nxl_on_file !< |
---|
| 2088 | INTEGER(iwp) :: nxrc !< |
---|
| 2089 | INTEGER(iwp) :: nxrf !< |
---|
| 2090 | INTEGER(iwp) :: nxr_on_file !< |
---|
| 2091 | INTEGER(iwp) :: nync !< |
---|
| 2092 | INTEGER(iwp) :: nynf !< |
---|
| 2093 | INTEGER(iwp) :: nyn_on_file !< |
---|
| 2094 | INTEGER(iwp) :: nysc !< |
---|
| 2095 | INTEGER(iwp) :: nysf !< |
---|
| 2096 | INTEGER(iwp) :: nys_on_file !< |
---|
| 2097 | |
---|
| 2098 | LOGICAL, INTENT(OUT) :: found |
---|
| 2099 | |
---|
[4797] | 2100 | REAL(wp), & |
---|
| 2101 | DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) & |
---|
| 2102 | :: tmp_3d !< |
---|
[3294] | 2103 | |
---|
| 2104 | |
---|
| 2105 | found = .TRUE. |
---|
| 2106 | |
---|
| 2107 | SELECT CASE ( restart_string(1:length) ) |
---|
| 2108 | |
---|
| 2109 | CASE ( 'rho_ocean_av' ) |
---|
| 2110 | IF ( .NOT. ALLOCATED( rho_ocean_av ) ) THEN |
---|
| 2111 | ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 2112 | ENDIF |
---|
| 2113 | IF ( k == 1 ) READ ( 13 ) tmp_3d |
---|
[4797] | 2114 | rho_ocean_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & |
---|
| 2115 | tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) |
---|
[3294] | 2116 | |
---|
| 2117 | CASE ( 'sa' ) |
---|
| 2118 | IF ( k == 1 ) READ ( 13 ) tmp_3d |
---|
[4797] | 2119 | sa(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & |
---|
| 2120 | tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) |
---|
[3294] | 2121 | |
---|
| 2122 | CASE ( 'sa_av' ) |
---|
| 2123 | IF ( .NOT. ALLOCATED( sa_av ) ) THEN |
---|
| 2124 | ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 2125 | ENDIF |
---|
| 2126 | IF ( k == 1 ) READ ( 13 ) tmp_3d |
---|
[4797] | 2127 | sa_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & |
---|
| 2128 | tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) |
---|
[3294] | 2129 | |
---|
| 2130 | CASE DEFAULT |
---|
| 2131 | found = .FALSE. |
---|
| 2132 | |
---|
| 2133 | END SELECT |
---|
| 2134 | |
---|
[4517] | 2135 | END SUBROUTINE ocean_rrd_local_ftn |
---|
[3294] | 2136 | |
---|
| 2137 | |
---|
[4797] | 2138 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 2139 | ! Description: |
---|
| 2140 | ! ------------ |
---|
[4517] | 2141 | !> Read module-specific local restart data arrays (MPI-IO). |
---|
[4797] | 2142 | !--------------------------------------------------------------------------------------------------! |
---|
[4517] | 2143 | SUBROUTINE ocean_rrd_local_mpi |
---|
| 2144 | |
---|
[4797] | 2145 | USE averaging, & |
---|
[4517] | 2146 | ONLY: rho_ocean_av, sa_av |
---|
| 2147 | |
---|
[4797] | 2148 | USE indices, & |
---|
[4517] | 2149 | ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt |
---|
| 2150 | |
---|
| 2151 | IMPLICIT NONE |
---|
| 2152 | |
---|
| 2153 | LOGICAL :: array_found !< |
---|
| 2154 | |
---|
| 2155 | |
---|
| 2156 | CALL rd_mpi_io_check_array( 'rho_ocean_av' , found = array_found ) |
---|
| 2157 | IF ( array_found ) THEN |
---|
[4797] | 2158 | IF ( .NOT. ALLOCATED( rho_ocean_av ) ) & |
---|
| 2159 | ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
[4517] | 2160 | CALL rrd_mpi_io( 'rho_ocean_av', rho_ocean_av ) |
---|
| 2161 | ENDIF |
---|
| 2162 | |
---|
| 2163 | CALL rrd_mpi_io( 'sa', sa ) |
---|
| 2164 | |
---|
| 2165 | CALL rd_mpi_io_check_array( 'sa_av' , found = array_found ) |
---|
| 2166 | IF ( array_found ) THEN |
---|
| 2167 | IF ( .NOT. ALLOCATED( sa_av ) ) ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) |
---|
| 2168 | CALL rrd_mpi_io( 'sa_av', sa_av ) |
---|
| 2169 | ENDIF |
---|
| 2170 | |
---|
| 2171 | END SUBROUTINE ocean_rrd_local_mpi |
---|
| 2172 | |
---|
| 2173 | |
---|
[4797] | 2174 | !--------------------------------------------------------------------------------------------------! |
---|
[4517] | 2175 | ! Description: |
---|
| 2176 | ! ------------ |
---|
[3294] | 2177 | !> This routine writes the respective restart data for the ocean module. |
---|
[4797] | 2178 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 2179 | SUBROUTINE ocean_wrd_global |
---|
| 2180 | |
---|
| 2181 | IMPLICIT NONE |
---|
| 2182 | |
---|
| 2183 | |
---|
[4495] | 2184 | IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN |
---|
[3294] | 2185 | |
---|
[4495] | 2186 | CALL wrd_write_string( 'bc_sa_t' ) |
---|
| 2187 | WRITE ( 14 ) bc_sa_t |
---|
[3303] | 2188 | |
---|
[4495] | 2189 | CALL wrd_write_string( 'bottom_salinityflux' ) |
---|
| 2190 | WRITE ( 14 ) bottom_salinityflux |
---|
[3294] | 2191 | |
---|
[4495] | 2192 | CALL wrd_write_string( 'salinity' ) |
---|
| 2193 | WRITE ( 14 ) salinity |
---|
[3294] | 2194 | |
---|
[4495] | 2195 | CALL wrd_write_string( 'sa_init' ) |
---|
| 2196 | WRITE ( 14 ) sa_init |
---|
[3294] | 2197 | |
---|
[4495] | 2198 | CALL wrd_write_string( 'sa_surface' ) |
---|
| 2199 | WRITE ( 14 ) sa_surface |
---|
[3294] | 2200 | |
---|
[4495] | 2201 | CALL wrd_write_string( 'sa_vertical_gradient' ) |
---|
| 2202 | WRITE ( 14 ) sa_vertical_gradient |
---|
[3294] | 2203 | |
---|
[4495] | 2204 | CALL wrd_write_string( 'sa_vertical_gradient_level' ) |
---|
| 2205 | WRITE ( 14 ) sa_vertical_gradient_level |
---|
[3294] | 2206 | |
---|
[4495] | 2207 | CALL wrd_write_string( 'stokes_waveheight' ) |
---|
| 2208 | WRITE ( 14 ) stokes_waveheight |
---|
[3381] | 2209 | |
---|
[4495] | 2210 | CALL wrd_write_string( 'stokes_wavelength' ) |
---|
| 2211 | WRITE ( 14 ) stokes_wavelength |
---|
[3294] | 2212 | |
---|
[4495] | 2213 | CALL wrd_write_string( 'surface_cooling_spinup_time' ) |
---|
| 2214 | WRITE ( 14 ) surface_cooling_spinup_time |
---|
[3294] | 2215 | |
---|
[4495] | 2216 | CALL wrd_write_string( 'top_salinityflux' ) |
---|
| 2217 | WRITE ( 14 ) top_salinityflux |
---|
[3302] | 2218 | |
---|
[4495] | 2219 | CALL wrd_write_string( 'wall_salinityflux' ) |
---|
| 2220 | WRITE ( 14 ) wall_salinityflux |
---|
| 2221 | |
---|
| 2222 | CALL wrd_write_string( 'wave_breaking' ) |
---|
| 2223 | WRITE ( 14 ) wave_breaking |
---|
| 2224 | |
---|
[4535] | 2225 | ELSEIF ( restart_data_format_output(1:3) == 'mpi' ) THEN |
---|
[4495] | 2226 | |
---|
| 2227 | CALL wrd_mpi_io( 'bc_sa_t', bc_sa_t ) |
---|
| 2228 | CALL wrd_mpi_io( 'bottom_salinityflux', bottom_salinityflux ) |
---|
| 2229 | CALL wrd_mpi_io( 'salinity', salinity ) |
---|
| 2230 | CALL wrd_mpi_io_global_array( 'sa_init', sa_init ) |
---|
| 2231 | CALL wrd_mpi_io( 'sa_surface', sa_surface ) |
---|
| 2232 | CALL wrd_mpi_io_global_array( 'sa_vertical_gradient', sa_vertical_gradient ) |
---|
| 2233 | CALL wrd_mpi_io_global_array( 'sa_vertical_gradient_level', sa_vertical_gradient_level ) |
---|
[4797] | 2234 | CALL wrd_mpi_io_global_array( 'sa_vertical_gradient_level_ind', & |
---|
| 2235 | sa_vertical_gradient_level_ind ) |
---|
[4495] | 2236 | CALL wrd_mpi_io( 'stokes_waveheight', stokes_waveheight ) |
---|
| 2237 | CALL wrd_mpi_io( 'stokes_wavelength', stokes_wavelength ) |
---|
| 2238 | CALL wrd_mpi_io( 'surface_cooling_spinup_time', surface_cooling_spinup_time ) |
---|
| 2239 | CALL wrd_mpi_io( 'top_salinityflux', top_salinityflux ) |
---|
| 2240 | CALL wrd_mpi_io_global_array( 'wall_salinityflux', wall_salinityflux ) |
---|
| 2241 | CALL wrd_mpi_io( 'wave_breaking', wave_breaking ) |
---|
| 2242 | |
---|
| 2243 | ENDIF |
---|
| 2244 | |
---|
[3294] | 2245 | END SUBROUTINE ocean_wrd_global |
---|
| 2246 | |
---|
| 2247 | |
---|
[4797] | 2248 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 2249 | ! Description: |
---|
| 2250 | ! ------------ |
---|
| 2251 | !> This routine writes the respective restart data for the ocean module. |
---|
[4797] | 2252 | !--------------------------------------------------------------------------------------------------! |
---|
[3294] | 2253 | SUBROUTINE ocean_wrd_local |
---|
| 2254 | |
---|
[4797] | 2255 | USE averaging, & |
---|
[3294] | 2256 | ONLY: rho_ocean_av, sa_av |
---|
| 2257 | |
---|
| 2258 | IMPLICIT NONE |
---|
| 2259 | |
---|
| 2260 | |
---|
[4495] | 2261 | IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN |
---|
[3294] | 2262 | |
---|
[4495] | 2263 | IF ( ALLOCATED( rho_ocean_av ) ) THEN |
---|
| 2264 | CALL wrd_write_string( 'rho_ocean_av' ) |
---|
| 2265 | WRITE ( 14 ) rho_ocean_av |
---|
| 2266 | ENDIF |
---|
| 2267 | |
---|
| 2268 | CALL wrd_write_string( 'sa' ) |
---|
| 2269 | WRITE ( 14 ) sa |
---|
| 2270 | |
---|
| 2271 | IF ( ALLOCATED( sa_av ) ) THEN |
---|
| 2272 | CALL wrd_write_string( 'sa_av' ) |
---|
| 2273 | WRITE ( 14 ) sa_av |
---|
| 2274 | ENDIF |
---|
| 2275 | |
---|
[4535] | 2276 | ELSEIF ( restart_data_format_output(1:3) == 'mpi' ) THEN |
---|
[4495] | 2277 | |
---|
| 2278 | IF ( ALLOCATED( rho_ocean_av ) ) CALL wrd_mpi_io( 'rho_ocean_av', rho_ocean_av ) |
---|
| 2279 | CALL wrd_mpi_io( 'sa', sa ) |
---|
| 2280 | IF ( ALLOCATED( sa_av ) ) CALL wrd_mpi_io( 'sa_av', sa_av ) |
---|
| 2281 | |
---|
[3294] | 2282 | ENDIF |
---|
| 2283 | |
---|
| 2284 | END SUBROUTINE ocean_wrd_local |
---|
| 2285 | |
---|
| 2286 | |
---|
[4797] | 2287 | !--------------------------------------------------------------------------------------------------! |
---|
[3302] | 2288 | ! Description: |
---|
| 2289 | ! ------------ |
---|
[4797] | 2290 | !> This routine calculates the Craik Leibovich vortex force and the additional effect of the Stokes |
---|
| 2291 | !> drift on the Coriolis force. |
---|
[3302] | 2292 | !> Call for all gridpoints. |
---|
[4797] | 2293 | !--------------------------------------------------------------------------------------------------! |
---|
[3302] | 2294 | SUBROUTINE stokes_drift_terms( component ) |
---|
| 2295 | |
---|
[4797] | 2296 | USE arrays_3d, & |
---|
| 2297 | ONLY: ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu, v_stokes_zw, w, tend |
---|
[3302] | 2298 | |
---|
[4797] | 2299 | USE basic_constants_and_equations_mod, & |
---|
[4196] | 2300 | ONLY: pi |
---|
| 2301 | |
---|
[4797] | 2302 | USE control_parameters, & |
---|
[4196] | 2303 | ONLY: f, fs, message_string, rotation_angle |
---|
[3302] | 2304 | |
---|
[4797] | 2305 | USE grid_variables, & |
---|
[3302] | 2306 | ONLY: ddx, ddy |
---|
| 2307 | |
---|
[4797] | 2308 | USE indices, & |
---|
[3302] | 2309 | ONLY: nxl, nxr, nys, nysv, nyn, nzb, nzt |
---|
| 2310 | |
---|
| 2311 | IMPLICIT NONE |
---|
| 2312 | |
---|
[4196] | 2313 | INTEGER(iwp) :: component !< component of momentum equation |
---|
| 2314 | INTEGER(iwp) :: i !< loop index along x |
---|
| 2315 | INTEGER(iwp) :: j !< loop index along y |
---|
| 2316 | INTEGER(iwp) :: k !< loop index along z |
---|
[3302] | 2317 | |
---|
[4196] | 2318 | REAL(wp) :: cos_rot_angle !< cosine of model rotation angle |
---|
| 2319 | REAL(wp) :: sin_rot_angle !< sine of model rotation angle |
---|
[3302] | 2320 | |
---|
[4797] | 2321 | |
---|
[3302] | 2322 | ! |
---|
| 2323 | !-- Compute Stokes terms for the respective velocity components |
---|
| 2324 | SELECT CASE ( component ) |
---|
| 2325 | |
---|
| 2326 | ! |
---|
| 2327 | !-- u-component |
---|
| 2328 | CASE ( 1 ) |
---|
| 2329 | DO i = nxl, nxr |
---|
| 2330 | DO j = nysv, nyn |
---|
| 2331 | DO k = nzb+1, nzt |
---|
[4797] | 2332 | tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * ( & |
---|
| 2333 | 0.5 * ( v(k,j+1,i) - v(k,j+1,i-1) & |
---|
| 2334 | + v(k,j,i) - v(k,j,i-1) ) * ddx & |
---|
| 2335 | - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) ) * ddy & |
---|
| 2336 | ) & |
---|
[3302] | 2337 | + f * v_stokes_zu(k) |
---|
| 2338 | ENDDO |
---|
| 2339 | ENDDO |
---|
| 2340 | ENDDO |
---|
| 2341 | |
---|
| 2342 | ! |
---|
| 2343 | !-- v-component |
---|
| 2344 | CASE ( 2 ) |
---|
| 2345 | DO i = nxl, nxr |
---|
| 2346 | DO j = nysv, nyn |
---|
| 2347 | DO k = nzb+1, nzt |
---|
[4797] | 2348 | tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * ( & |
---|
| 2349 | 0.5 * ( v(k,j,i+1) - v(k,j,i-1) ) * ddx & |
---|
| 2350 | - 0.5 * ( u(k,j,i) - u(k,j-1,i) & |
---|
| 2351 | + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & |
---|
[3302] | 2352 | ) & |
---|
| 2353 | - f * u_stokes_zu(k) |
---|
| 2354 | ENDDO |
---|
| 2355 | ENDDO |
---|
| 2356 | ENDDO |
---|
| 2357 | |
---|
| 2358 | ! |
---|
| 2359 | !-- w-component |
---|
| 2360 | CASE ( 3 ) |
---|
[4196] | 2361 | |
---|
| 2362 | ! |
---|
| 2363 | !-- Precalculate cosine and sine of rotation angle |
---|
| 2364 | cos_rot_angle = COS( rotation_angle * pi / 180.0_wp ) |
---|
| 2365 | sin_rot_angle = SIN( rotation_angle * pi / 180.0_wp ) |
---|
| 2366 | |
---|
[3302] | 2367 | DO i = nxl, nxr |
---|
| 2368 | DO j = nys, nyn |
---|
| 2369 | DO k = nzb+1, nzt |
---|
[4797] | 2370 | tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * ( & |
---|
| 2371 | 0.5 * ( u(k+1,j,i) - u(k,j,i) & |
---|
| 2372 | + u(k+1,j,i+1) - u(k,j,i+1) & |
---|
| 2373 | ) * ddzu(k+1) & |
---|
| 2374 | - 0.5 * ( w(k,j,i+1) - w(k,j,i-1) & |
---|
| 2375 | ) * ddx ) & |
---|
| 2376 | - v_stokes_zw(k) * ( & |
---|
| 2377 | 0.5 * ( w(k,j+1,i) - w(k,j-1,i) & |
---|
| 2378 | ) * ddy & |
---|
| 2379 | - 0.5 * ( v(k+1,j,i) - v(k,j,i) & |
---|
| 2380 | + v(k+1,j+1,i) - v(k,j+1,i) & |
---|
| 2381 | ) * ddzu(k) ) & |
---|
| 2382 | + fs * ( sin_rot_angle * v_stokes_zw(k) & |
---|
| 2383 | + cos_rot_angle * u_stokes_zw(k) & |
---|
| 2384 | ) |
---|
[3302] | 2385 | ENDDO |
---|
| 2386 | ENDDO |
---|
| 2387 | ENDDO |
---|
| 2388 | |
---|
| 2389 | CASE DEFAULT |
---|
[4797] | 2390 | WRITE( message_string, * ) 'wrong component of Stokes force: ', component |
---|
[3302] | 2391 | CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 ) |
---|
| 2392 | |
---|
| 2393 | END SELECT |
---|
| 2394 | |
---|
| 2395 | END SUBROUTINE stokes_drift_terms |
---|
| 2396 | |
---|
| 2397 | |
---|
[4797] | 2398 | !--------------------------------------------------------------------------------------------------! |
---|
[3302] | 2399 | ! Description: |
---|
| 2400 | ! ------------ |
---|
[4797] | 2401 | !> This routine calculates the Craik Leibovich vortex force and the additional effect of the Stokes |
---|
| 2402 | !> drift on the Coriolis force. |
---|
[3302] | 2403 | !> Call for gridpoints i,j. |
---|
[4797] | 2404 | !--------------------------------------------------------------------------------------------------! |
---|
[3302] | 2405 | |
---|
| 2406 | SUBROUTINE stokes_drift_terms_ij( i, j, component ) |
---|
| 2407 | |
---|
[4797] | 2408 | USE arrays_3d, & |
---|
| 2409 | ONLY: ddzu, u, u_stokes_zu, u_stokes_zw, v, v_stokes_zu, v_stokes_zw, w, tend |
---|
[3302] | 2410 | |
---|
[4797] | 2411 | USE basic_constants_and_equations_mod, & |
---|
[4196] | 2412 | ONLY: pi |
---|
| 2413 | |
---|
[4797] | 2414 | USE control_parameters, & |
---|
[4196] | 2415 | ONLY: f, fs, message_string, rotation_angle |
---|
[3302] | 2416 | |
---|
[4797] | 2417 | USE grid_variables, & |
---|
[3302] | 2418 | ONLY: ddx, ddy |
---|
| 2419 | |
---|
[4797] | 2420 | USE indices, & |
---|
[3614] | 2421 | ONLY: nzb, nzt |
---|
[3302] | 2422 | |
---|
| 2423 | IMPLICIT NONE |
---|
| 2424 | |
---|
[4196] | 2425 | INTEGER(iwp) :: component !< component of momentum equation |
---|
| 2426 | INTEGER(iwp) :: i !< loop index along x |
---|
| 2427 | INTEGER(iwp) :: j !< loop index along y |
---|
| 2428 | INTEGER(iwp) :: k !< loop incex along z |
---|
[3302] | 2429 | |
---|
[4196] | 2430 | REAL(wp) :: cos_rot_angle !< cosine of model rotation angle |
---|
| 2431 | REAL(wp) :: sin_rot_angle !< sine of model rotation angle |
---|
[3302] | 2432 | |
---|
[4797] | 2433 | |
---|
[3302] | 2434 | ! |
---|
| 2435 | !-- Compute Stokes terms for the respective velocity components |
---|
| 2436 | SELECT CASE ( component ) |
---|
| 2437 | |
---|
| 2438 | ! |
---|
| 2439 | !-- u-component |
---|
| 2440 | CASE ( 1 ) |
---|
| 2441 | DO k = nzb+1, nzt |
---|
[4797] | 2442 | tend(k,j,i) = tend(k,j,i) + v_stokes_zu(k) * ( & |
---|
| 2443 | 0.5 * ( v(k,j+1,i) - v(k,j+1,i-1) & |
---|
| 2444 | + v(k,j,i) - v(k,j,i-1) ) * ddx & |
---|
| 2445 | - 0.5 * ( u(k,j+1,i) - u(k,j-1,i) ) * ddy & |
---|
| 2446 | ) & |
---|
[3302] | 2447 | + f * v_stokes_zu(k) |
---|
| 2448 | ENDDO |
---|
| 2449 | ! |
---|
| 2450 | !-- v-component |
---|
| 2451 | CASE ( 2 ) |
---|
| 2452 | DO k = nzb+1, nzt |
---|
[4797] | 2453 | tend(k,j,i) = tend(k,j,i) - u_stokes_zu(k) * ( & |
---|
| 2454 | 0.5 * ( v(k,j,i+1) - v(k,j,i-1) ) * ddx & |
---|
| 2455 | - 0.5 * ( u(k,j,i) - u(k,j-1,i) & |
---|
| 2456 | + u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & |
---|
| 2457 | ) & |
---|
[3302] | 2458 | - f * u_stokes_zu(k) |
---|
| 2459 | ENDDO |
---|
| 2460 | |
---|
| 2461 | ! |
---|
| 2462 | !-- w-component |
---|
| 2463 | CASE ( 3 ) |
---|
[4196] | 2464 | |
---|
| 2465 | ! |
---|
| 2466 | !-- Precalculate cosine and sine of rotation angle |
---|
| 2467 | cos_rot_angle = COS( rotation_angle * pi / 180.0_wp ) |
---|
| 2468 | sin_rot_angle = SIN( rotation_angle * pi / 180.0_wp ) |
---|
| 2469 | |
---|
[3302] | 2470 | DO k = nzb+1, nzt |
---|
[4797] | 2471 | tend(k,j,i) = tend(k,j,i) + u_stokes_zw(k) * ( 0.5 * ( u(k+1,j,i) - u(k,j,i) & |
---|
| 2472 | + u(k+1,j,i+1) - u(k,j,i+1) & |
---|
| 2473 | ) * ddzu(k+1) & |
---|
| 2474 | - 0.5 * ( w(k,j,i+1) - w(k,j,i-1) & |
---|
| 2475 | ) * ddx ) & |
---|
| 2476 | - v_stokes_zw(k) * ( 0.5 * ( w(k,j+1,i) - w(k,j-1,i) & |
---|
| 2477 | ) * ddy & |
---|
| 2478 | - 0.5 * ( v(k+1,j,i) - v(k,j,i) & |
---|
| 2479 | + v(k+1,j+1,i) - v(k,j+1,i) & |
---|
| 2480 | ) * ddzu(k) ) & |
---|
| 2481 | + fs * ( sin_rot_angle * v_stokes_zw(k) & |
---|
| 2482 | + cos_rot_angle * u_stokes_zw(k) & |
---|
[4196] | 2483 | ) |
---|
[3302] | 2484 | ENDDO |
---|
| 2485 | |
---|
| 2486 | CASE DEFAULT |
---|
| 2487 | WRITE( message_string, * ) ' wrong component: ', component |
---|
| 2488 | CALL message( 'stokes_drift_terms', 'PA0091', 1, 2, 0, 6, 0 ) |
---|
| 2489 | |
---|
| 2490 | END SELECT |
---|
| 2491 | |
---|
| 2492 | END SUBROUTINE stokes_drift_terms_ij |
---|
| 2493 | |
---|
| 2494 | |
---|
[4797] | 2495 | !--------------------------------------------------------------------------------------------------! |
---|
[3302] | 2496 | ! Description: |
---|
| 2497 | ! ------------ |
---|
[4797] | 2498 | !> This routine calculates turbulence generated by wave breaking near the ocean surface, following |
---|
| 2499 | !> a parameterization given in Noh et al. (2004), JPO |
---|
[3302] | 2500 | !> Call for all gridpoints. |
---|
[4797] | 2501 | !> TODO: so far, this routine only works if the model time step has about the same value as the time |
---|
| 2502 | !> scale of wave breaking! |
---|
| 2503 | !--------------------------------------------------------------------------------------------------! |
---|
[3302] | 2504 | SUBROUTINE wave_breaking_term( component ) |
---|
| 2505 | |
---|
[4797] | 2506 | USE arrays_3d, & |
---|
[3302] | 2507 | ONLY: u_p, v_p |
---|
| 2508 | |
---|
[4797] | 2509 | USE control_parameters, & |
---|
[3302] | 2510 | ONLY: dt_3d, message_string |
---|
| 2511 | |
---|
[4797] | 2512 | USE indices, & |
---|
[3302] | 2513 | ONLY: nxl, nxlu, nxr, nys, nysv, nyn, nzt |
---|
| 2514 | |
---|
| 2515 | IMPLICIT NONE |
---|
| 2516 | |
---|
| 2517 | INTEGER(iwp) :: component !< component of momentum equation |
---|
| 2518 | INTEGER(iwp) :: i !< loop index along x |
---|
| 2519 | INTEGER(iwp) :: j !< loop index along y |
---|
| 2520 | |
---|
[4797] | 2521 | REAL(wp) :: random_gauss !< function that creates a random number with a Gaussian distribution |
---|
[3302] | 2522 | |
---|
| 2523 | |
---|
| 2524 | ! |
---|
| 2525 | !-- Compute wave breaking terms for the respective velocity components. |
---|
| 2526 | !-- Velocities are directly manipulated, since this is not a real force |
---|
| 2527 | SELECT CASE ( component ) |
---|
| 2528 | |
---|
| 2529 | ! |
---|
| 2530 | !-- u-component |
---|
| 2531 | CASE ( 1 ) |
---|
| 2532 | DO i = nxlu, nxr |
---|
| 2533 | DO j = nys, nyn |
---|
[4797] | 2534 | u_p(nzt,j,i) = u_p(nzt,j,i) + & |
---|
| 2535 | ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) & |
---|
| 2536 | * alpha_wave_breaking * u_star_wave_breaking & |
---|
[3302] | 2537 | / timescale_wave_breaking * dt_3d |
---|
| 2538 | ENDDO |
---|
| 2539 | ENDDO |
---|
| 2540 | ! |
---|
| 2541 | !-- v-component |
---|
| 2542 | CASE ( 2 ) |
---|
| 2543 | DO i = nxl, nxr |
---|
| 2544 | DO j = nysv, nyn |
---|
[4797] | 2545 | v_p(nzt,j,i) = v_p(nzt,j,i) + & |
---|
| 2546 | ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) & |
---|
| 2547 | * alpha_wave_breaking * u_star_wave_breaking & |
---|
[3302] | 2548 | / timescale_wave_breaking * dt_3d |
---|
| 2549 | ENDDO |
---|
| 2550 | ENDDO |
---|
| 2551 | |
---|
| 2552 | CASE DEFAULT |
---|
[4797] | 2553 | WRITE( message_string, * ) 'wrong component of wave breaking: ', component |
---|
[3302] | 2554 | CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 ) |
---|
| 2555 | |
---|
| 2556 | END SELECT |
---|
| 2557 | |
---|
| 2558 | END SUBROUTINE wave_breaking_term |
---|
| 2559 | |
---|
| 2560 | |
---|
[4797] | 2561 | !--------------------------------------------------------------------------------------------------! |
---|
[3302] | 2562 | ! Description: |
---|
| 2563 | ! ------------ |
---|
[4797] | 2564 | !> This routine calculates turbulence generated by wave breaking near the ocean surface, following a |
---|
| 2565 | !> parameterization given in Noh et al. (2004), JPO |
---|
[3302] | 2566 | !> Call for gridpoint i,j. |
---|
[4797] | 2567 | !> TODO: so far, this routine only works if the model time step has about the same value as the time |
---|
| 2568 | !> scale of wave breaking! |
---|
| 2569 | !--------------------------------------------------------------------------------------------------! |
---|
[3302] | 2570 | SUBROUTINE wave_breaking_term_ij( i, j, component ) |
---|
| 2571 | |
---|
[4797] | 2572 | USE arrays_3d, & |
---|
[3302] | 2573 | ONLY: u_p, v_p |
---|
| 2574 | |
---|
[4797] | 2575 | USE control_parameters, & |
---|
[3302] | 2576 | ONLY: dt_3d, message_string |
---|
| 2577 | |
---|
[4797] | 2578 | USE indices, & |
---|
[3302] | 2579 | ONLY: nzt |
---|
| 2580 | |
---|
| 2581 | IMPLICIT NONE |
---|
| 2582 | |
---|
| 2583 | INTEGER(iwp) :: component !< component of momentum equation |
---|
| 2584 | INTEGER(iwp) :: i !< loop index along x |
---|
| 2585 | INTEGER(iwp) :: j !< loop index along y |
---|
| 2586 | |
---|
[4797] | 2587 | REAL(wp) :: random_gauss !< function that creates a random number with a Gaussian distribution |
---|
[3302] | 2588 | |
---|
| 2589 | ! |
---|
| 2590 | !-- Compute wave breaking terms for the respective velocity components |
---|
| 2591 | SELECT CASE ( component ) |
---|
| 2592 | |
---|
| 2593 | ! |
---|
| 2594 | !-- u-/v-component |
---|
| 2595 | CASE ( 1 ) |
---|
[4797] | 2596 | u_p(nzt,j,i) = u_p(nzt,j,i) + & |
---|
| 2597 | ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) & |
---|
| 2598 | * alpha_wave_breaking * u_star_wave_breaking & |
---|
[3302] | 2599 | / timescale_wave_breaking * dt_3d |
---|
| 2600 | |
---|
| 2601 | CASE ( 2 ) |
---|
[4797] | 2602 | v_p(nzt,j,i) = v_p(nzt,j,i) + & |
---|
| 2603 | ( random_gauss( iran_ocean, 1.0_wp ) - 1.0_wp ) & |
---|
| 2604 | * alpha_wave_breaking * u_star_wave_breaking & |
---|
[3302] | 2605 | / timescale_wave_breaking * dt_3d |
---|
| 2606 | |
---|
| 2607 | CASE DEFAULT |
---|
[4797] | 2608 | WRITE( message_string, * ) 'wrong component of wave breaking: ', component |
---|
[3302] | 2609 | CALL message( 'stokes_drift_terms', 'PA0466', 1, 2, 0, 6, 0 ) |
---|
| 2610 | |
---|
| 2611 | END SELECT |
---|
| 2612 | |
---|
| 2613 | END SUBROUTINE wave_breaking_term_ij |
---|
| 2614 | |
---|
| 2615 | |
---|
[3294] | 2616 | END MODULE ocean_mod |
---|