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