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