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