Changeset 4562 for palm/trunk
- Timestamp:
- Jun 12, 2020 8:38:47 AM (4 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
check_parameters.f90 (modified) (2 diffs)
-
surface_layer_fluxes_mod.f90 (modified) (9 diffs)
-
synthetic_turbulence_generator_mod.f90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/check_parameters.f90
r4559 r4562 24 24 ! ----------------- 25 25 ! $Id$ 26 ! bugfix: revised error message for exceeding allow number of time series 27 ! 28 ! 4559 2020-06-11 08:51:48Z raasch 26 29 ! file re-formatted to follow the PALM coding standard 27 30 ! … … 671 674 WRITE( message_string, * ) 'number of time series quantities exceeds', & 672 675 ' its maximum of dots_max = ', dots_max, & 673 '&Please increase dots_max in modules.f90.'674 CALL message( ' init_3d_model', 'PA0194', 1, 2, 0, 6, 0 )676 '&Please increase dots_max in netcdf_interface_mod.f90.' 677 CALL message( 'check_parameters', 'PA0194', 1, 2, 0, 6, 0 ) 675 678 ENDIF 676 679 -
palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
r4519 r4562 1 1 !> @file surface_layer_fluxes_mod.f90 2 !------------------------------------------------------------------------------ !2 !--------------------------------------------------------------------------------------------------! 3 3 ! This file is part of the PALM model system. 4 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/>. 5 ! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General 6 ! Public License as published by the Free Software Foundation, either version 3 of the License, or 7 ! (at your option) any later version. 8 ! 9 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the 10 ! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 11 ! Public License for more details. 12 ! 13 ! You should have received a copy of the GNU General Public License along with PALM. If not, see 14 ! <http://www.gnu.org/licenses/>. 16 15 ! 17 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 18 ! 19 ! ------------------------------------------------------------------------------!17 !--------------------------------------------------------------------------------------------------! 18 ! 20 19 ! 21 20 ! Current revisions: 22 ! ----------------- -21 ! ----------------- 23 22 ! 24 23 ! … … 26 25 ! ----------------- 27 26 ! $Id$ 27 ! File re-formatted to follow the PALM coding standard 28 ! 29 ! 4519 2020-05-05 17:33:30Z suehring 28 30 ! Add missing computation of passive scalar scaling parameter 29 ! 31 ! 30 32 ! 4370 2020-01-10 14:00:44Z raasch 31 ! bugfix: openacc porting for vector version of OL calculation added32 ! 33 ! Bugfix: openacc porting for vector version of OL calculation added 34 ! 33 35 ! 4366 2020-01-09 08:12:43Z raasch 34 ! vector version for calculation of Obukhov length via Newton iteration added35 ! 36 ! Vector version for calculation of Obukhov length via Newton iteration added 37 ! 36 38 ! 4360 2020-01-07 11:25:50Z suehring 37 ! Calculation of diagnostic-only 2-m potential temperature moved to 38 ! diagnostic_output_quantities. 39 ! 39 ! Calculation of diagnostic-only 2-m potential temperature moved to diagnostic_output_quantities. 40 ! 40 41 ! 4298 2019-11-21 15:59:16Z suehring 41 ! Calculation of 2-m temperature adjusted to the case the 2-m level is above 42 ! the first grid point. 43 ! 42 ! Calculation of 2-m temperature adjusted to the case the 2-m level is above the first grid point. 43 ! 44 44 ! 4258 2019-10-07 13:29:08Z suehring 45 45 ! Initialization of Obukhov lenght also at vertical surfaces (if allocated). 46 ! 46 ! 47 47 ! 4237 2019-09-25 11:33:42Z knoop 48 48 ! Added missing OpenMP directives 49 ! 49 ! 50 50 ! 4186 2019-08-23 16:06:14Z suehring 51 ! - To enable limitation of Obukhov length, move it before exit-cycle construct. 52 ! Further, change the limit to 10E-5 in order to get rid-off unrealistic 53 ! peaks in the heat fluxesduring nighttime51 ! - To enable limitation of Obukhov length, move it before exit-cycle construct. 52 ! Further, change the limit to 10E-5 in order to get rid-off unrealistic peaks in the heat fluxes 53 ! during nighttime 54 54 ! - Unused variable removed 55 ! 55 ! 56 56 ! 4182 2019-08-22 15:20:23Z scharf 57 57 ! Corrected "Former revisions" section 58 ! 58 ! 59 59 ! 3987 2019-05-22 09:52:13Z kanani 60 60 ! Introduce alternative switch for debug output during timestepping 61 ! 61 ! 62 62 ! 3885 2019-04-11 11:29:34Z kanani 63 ! Changes related to global restructuring of location messages and introduction 64 ! of additional debugmessages65 ! 63 ! Changes related to global restructuring of location messages and introduction of additional debug 64 ! messages 65 ! 66 66 ! 3881 2019-04-10 09:31:22Z suehring 67 67 ! Assure that Obukhov length does not become zero 68 ! 68 ! 69 69 ! 3834 2019-03-28 15:40:15Z forkel 70 ! added USE chem_gasphase_mod71 ! 70 ! Added USE chem_gasphase_mod 71 ! 72 72 ! 3787 2019-03-07 08:43:54Z raasch 73 ! unused variables removed74 ! 73 ! Unused variables removed 74 ! 75 75 ! 3745 2019-02-15 18:57:56Z suehring 76 ! Bugfix, missing calculation of 10cm temperature at vertical building walls, 77 ! required for indoormodel78 ! 76 ! Bugfix, missing calculation of 10cm temperature at vertical building walls, required for indoor 77 ! model 78 ! 79 79 ! 3744 2019-02-15 18:38:58Z suehring 80 80 ! Some interface calls moved to module_interface + cleanup 81 ! 81 ! 82 82 ! 3668 2019-01-14 12:49:24Z maronga 83 83 ! Removed methods "circular" and "lookup" 84 ! 84 ! 85 85 ! 3655 2019-01-07 16:51:22Z knoop 86 86 ! OpenACC port for SPEC … … 92 92 ! Description: 93 93 ! ------------ 94 !> Diagnostic computation of vertical fluxes in the constant flux layer from the 95 !> va lues of the variables at grid point k=1 based on Newton iteration94 !> Diagnostic computation of vertical fluxes in the constant flux layer from the values of the 95 !> variables at grid point k=1 based on Newton iteration. 96 96 !> 97 !> @todo ( re)move large_scale_forcing actions98 !> @todo check/optimize OpenMP directives99 !> @todo simplify if conditions (which flux need to be computed in which case)100 !------------------------------------------------------------------------------ !97 !> @todo (Re)move large_scale_forcing actions 98 !> @todo Check/optimize OpenMP directives 99 !> @todo Simplify if conditions (which flux need to be computed in which case) 100 !--------------------------------------------------------------------------------------------------! 101 101 MODULE surface_layer_fluxes_mod 102 102 103 USE arrays_3d, & 104 ONLY: e, kh, nc, nr, pt, q, ql, qc, qr, s, u, v, vpt, w, zu, zw, & 105 drho_air_zw, rho_air_zw, d_exner 106 107 USE basic_constants_and_equations_mod, & 108 ONLY: g, kappa, lv_d_cp, pi, rd_d_rv 109 110 USE chem_gasphase_mod, & 103 USE arrays_3d, & 104 ONLY: d_exner, & 105 drho_air_zw, & 106 e, & 107 kh, & 108 nc, & 109 nr, & 110 pt, & 111 q, & 112 ql, & 113 qc, & 114 qr, & 115 s, & 116 u, & 117 v, & 118 vpt, & 119 w, & 120 zu, & 121 zw, & 122 rho_air_zw 123 124 125 USE basic_constants_and_equations_mod, & 126 ONLY: g, & 127 kappa, & 128 lv_d_cp, & 129 pi, & 130 rd_d_rv 131 132 USE chem_gasphase_mod, & 111 133 ONLY: nvar 112 134 113 USE chem_modules, &135 USE chem_modules, & 114 136 ONLY: constant_csflux 115 137 116 138 USE cpulog 117 139 118 USE control_parameters, & 119 ONLY: air_chemistry, cloud_droplets, & 120 constant_heatflux, constant_scalarflux, & 121 constant_waterflux, coupling_mode, & 122 debug_output_timestep, & 123 humidity, loop_optimization, & 124 ibc_e_b, ibc_pt_b, indoor_model, & 125 land_surface, large_scale_forcing, lsf_surf, message_string, & 126 neutral, passive_scalar, pt_surface, q_surface, & 127 run_coupled, surface_pressure, simulated_time, & 128 time_since_reference_point, urban_surface, & 129 use_free_convection_scaling, zeta_max, zeta_min 130 131 USE grid_variables, & 132 ONLY: dx, dy 133 134 USE indices, & 140 USE control_parameters, & 141 ONLY: air_chemistry, & 142 cloud_droplets, & 143 constant_heatflux, & 144 constant_scalarflux, & 145 constant_waterflux, & 146 coupling_mode, & 147 debug_output_timestep, & 148 humidity, & 149 ibc_e_b, & 150 ibc_pt_b, & 151 indoor_model, & 152 land_surface, & 153 large_scale_forcing, & 154 loop_optimization, & 155 lsf_surf, & 156 message_string, & 157 neutral, & 158 passive_scalar, & 159 pt_surface, & 160 q_surface, & 161 run_coupled, & 162 surface_pressure, & 163 simulated_time, & 164 time_since_reference_point, & 165 urban_surface, & 166 use_free_convection_scaling, & 167 zeta_max, & 168 zeta_min 169 170 USE grid_variables, & 171 ONLY: dx, & 172 dy 173 174 USE indices, & 135 175 ONLY: nzt 136 176 137 177 USE kinds 138 178 139 USE bulk_cloud_model_mod, & 140 ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert 179 USE bulk_cloud_model_mod, & 180 ONLY: bulk_cloud_model, & 181 microphysics_morrison, & 182 microphysics_seifert 141 183 142 184 USE pegrid 143 185 144 USE land_surface_model_mod, & 145 ONLY: aero_resist_kray, skip_time_do_lsm 146 147 USE surface_mod, & 148 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_type, & 149 surf_usm_h, surf_usm_v 150 186 USE land_surface_model_mod, & 187 ONLY: aero_resist_kray, & 188 skip_time_do_lsm 189 190 USE surface_mod, & 191 ONLY : surf_def_h, & 192 surf_def_v, & 193 surf_lsm_h, & 194 surf_lsm_v, & 195 surf_type, & 196 surf_usm_h, & 197 surf_usm_v 198 151 199 152 200 IMPLICIT NONE 153 201 154 INTEGER(iwp) :: i 155 INTEGER(iwp) :: j 156 INTEGER(iwp) :: k 157 INTEGER(iwp) :: l 158 159 LOGICAL :: coupled_run!< Flag for coupled atmosphere-ocean runs160 LOGICAL :: downward = .FALSE.!< Flag indicating downward-facing horizontal surface161 LOGICAL :: mom_uv = .FALSE.!< Flag indicating calculation of usvs and vsus at vertical surfaces162 LOGICAL :: mom_w = .FALSE.!< Flag indicating calculation of wsus and wsvs at vertical surfaces163 LOGICAL :: mom_tke = .FALSE. !< Flag indicating calculation of momentum fluxes at vertical surfaces used for TKE production164 LOGICAL :: surf_vertical!< Flag indicating vertical surfaces165 166 REAL(wp) :: e_s, &!< Saturation water vapor pressure167 ol_max = 1.0E6_wp, &!< Maximum Obukhov length168 z_mo!< Height of the constant flux layer where MOST is assumed169 170 TYPE(surf_type), POINTER :: surf 202 INTEGER(iwp) :: i !< loop index x direction 203 INTEGER(iwp) :: j !< loop index y direction 204 INTEGER(iwp) :: k !< loop index z direction 205 INTEGER(iwp) :: l !< loop index for surf type 206 207 LOGICAL :: coupled_run !< Flag for coupled atmosphere-ocean runs 208 LOGICAL :: downward = .FALSE. !< Flag indicating downward-facing horizontal surface 209 LOGICAL :: mom_uv = .FALSE. !< Flag indicating calculation of usvs and vsus at vertical surfaces 210 LOGICAL :: mom_w = .FALSE. !< Flag indicating calculation of wsus and wsvs at vertical surfaces 211 LOGICAL :: mom_tke = .FALSE. !< Flag indicating calculation of momentum fluxes at vertical surfaces used for TKE production 212 LOGICAL :: surf_vertical !< Flag indicating vertical surfaces 213 214 REAL(wp) :: e_s !< Saturation water vapor pressure 215 REAL(wp) :: ol_max = 1.0E6_wp !< Maximum Obukhov length 216 REAL(wp) :: z_mo !< Height of the constant flux layer where MOST is assumed 217 218 TYPE(surf_type), POINTER :: surf !< surf-type array, used to generalize subroutines 171 219 172 220 … … 175 223 PRIVATE 176 224 177 PUBLIC init_surface_layer_fluxes, &178 phi_m, &179 psi_h, &180 psi_m, &225 PUBLIC init_surface_layer_fluxes, & 226 phi_m, & 227 psi_h, & 228 psi_m, & 181 229 surface_layer_fluxes 182 230 … … 205 253 206 254 207 !------------------------------------------------------------------------------ !255 !--------------------------------------------------------------------------------------------------! 208 256 ! Description: 209 257 ! ------------ 210 !> Main routine to compute the surface fluxes 211 !------------------------------------------------------------------------------! 212 SUBROUTINE surface_layer_fluxes 213 214 IMPLICIT NONE 215 216 217 IF ( debug_output_timestep ) CALL debug_message( 'surface_layer_fluxes', 'start' ) 218 219 surf_vertical = .FALSE. !< flag indicating vertically orientated surface elements 220 downward = .FALSE. !< flag indicating downward-facing surface elements 221 ! 222 !-- Derive potential temperature and specific humidity at first grid level 223 !-- from the fields pt and q 224 ! 225 !-- First call for horizontal default-type surfaces (l=0 - upward facing, 226 !-- l=1 - downward facing) 227 DO l = 0, 1 228 IF ( surf_def_h(l)%ns >= 1 ) THEN 229 surf => surf_def_h(l) 230 CALL calc_pt_q 231 IF ( .NOT. neutral ) THEN 232 CALL calc_pt_surface 233 IF ( humidity ) THEN 234 CALL calc_q_surface 235 CALL calc_vpt_surface 236 ENDIF 258 !> Main routine to compute the surface fluxes. 259 !--------------------------------------------------------------------------------------------------! 260 SUBROUTINE surface_layer_fluxes 261 262 IMPLICIT NONE 263 264 265 IF ( debug_output_timestep ) CALL debug_message( 'surface_layer_fluxes', 'start' ) 266 267 surf_vertical = .FALSE. !< flag indicating vertically orientated surface elements 268 downward = .FALSE. !< flag indicating downward-facing surface elements 269 ! 270 !-- Derive potential temperature and specific humidity at first grid level from the fields pt and q 271 ! 272 !-- First call for horizontal default-type surfaces (l=0 - upward facing, l=1 - downward facing) 273 DO l = 0, 1 274 IF ( surf_def_h(l)%ns >= 1 ) THEN 275 surf => surf_def_h(l) 276 CALL calc_pt_q 277 IF ( .NOT. neutral ) THEN 278 CALL calc_pt_surface 279 IF ( humidity ) THEN 280 CALL calc_q_surface 281 CALL calc_vpt_surface 237 282 ENDIF 238 283 ENDIF 239 ENDDO 240 ! 241 !-- Call for natural-type horizontal surfaces 242 IF ( surf_lsm_h%ns >= 1 ) THEN 243 surf => surf_lsm_h 284 ENDIF 285 ENDDO 286 ! 287 !-- Call for natural-type horizontal surfaces 288 IF ( surf_lsm_h%ns >= 1 ) THEN 289 surf => surf_lsm_h 290 CALL calc_pt_q 291 ENDIF 292 293 ! 294 !-- Call for urban-type horizontal surfaces 295 IF ( surf_usm_h%ns >= 1 ) THEN 296 surf => surf_usm_h 297 CALL calc_pt_q 298 ENDIF 299 300 ! 301 !-- Call for natural-type vertical surfaces 302 DO l = 0, 3 303 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 304 surf => surf_lsm_v(l) 244 305 CALL calc_pt_q 245 306 ENDIF 246 307 247 ! 248 !-- Call for urban-type horizontal surfaces 249 IF ( surf_usm_h%ns >= 1 ) THEN 250 surf => surf_usm_h 308 !-- Call for urban-type vertical surfaces 309 IF ( surf_usm_v(l)%ns >= 1 ) THEN 310 surf => surf_usm_v(l) 251 311 CALL calc_pt_q 252 312 ENDIF 253 254 ! 255 !-- Call for natural-type vertical surfaces 256 DO l = 0, 3 257 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 313 ENDDO 314 315 ! 316 !-- First, calculate the new Obukhov length, then new friction velocity, followed by the new scaling 317 !-- parameters (th*, q*, etc.), and the new surface fluxes, if required. Note, each routine is called 318 !-- for different surface types. First call for default-type horizontal surfaces, for natural- and 319 !-- urban-type surfaces. Note, here only upward-facing horizontal surfaces are treated. 320 321 ! 322 !-- Default-type upward-facing horizontal surfaces 323 IF ( surf_def_h(0)%ns >= 1 ) THEN 324 surf => surf_def_h(0) 325 CALL calc_uvw_abs 326 IF ( .NOT. neutral ) CALL calc_ol 327 CALL calc_us 328 CALL calc_scaling_parameters 329 CALL calc_surface_fluxes 330 ENDIF 331 ! 332 !-- Natural-type horizontal surfaces 333 IF ( surf_lsm_h%ns >= 1 ) THEN 334 surf => surf_lsm_h 335 CALL calc_uvw_abs 336 IF ( .NOT. neutral ) CALL calc_ol 337 CALL calc_us 338 CALL calc_scaling_parameters 339 CALL calc_surface_fluxes 340 ENDIF 341 ! 342 !-- Urban-type horizontal surfaces 343 IF ( surf_usm_h%ns >= 1 ) THEN 344 surf => surf_usm_h 345 CALL calc_uvw_abs 346 IF ( .NOT. neutral ) CALL calc_ol 347 CALL calc_us 348 CALL calc_scaling_parameters 349 CALL calc_surface_fluxes 350 ! 351 !-- Calculate 10cm temperature, required in indoor model 352 IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' ) 353 ENDIF 354 355 ! 356 !-- Treat downward-facing horizontal surfaces. Note, so far, these are always default type. 357 !-- Stratification is not considered in this case, hence, no further distinction between different 358 !-- most_method is required. 359 IF ( surf_def_h(1)%ns >= 1 ) THEN 360 downward = .TRUE. 361 surf => surf_def_h(1) 362 CALL calc_uvw_abs 363 CALL calc_us 364 CALL calc_surface_fluxes 365 downward = .FALSE. 366 ENDIF 367 ! 368 !-- Calculate surfaces fluxes at vertical surfaces for momentum and subgrid-scale TKE. No stability 369 !-- is considered. Therefore, scaling parameters and Obukhov length do not need to be calculated and 370 !-- no distinction in 'circular', 'Newton' or 'lookup' is necessary so far. Note, this will change 371 !-- if stability is once considered. 372 surf_vertical = .TRUE. 373 ! 374 !-- Calculate horizontal momentum fluxes at north- and south-facing surfaces(usvs). 375 !-- For default-type surfaces 376 mom_uv = .TRUE. 377 DO l = 0, 1 378 IF ( surf_def_v(l)%ns >= 1 ) THEN 379 surf => surf_def_v(l) 380 ! 381 !-- Compute surface-parallel velocity 382 CALL calc_uvw_abs_v_ugrid 383 ! 384 !-- Compute respective friction velocity on staggered grid 385 CALL calc_us 386 ! 387 !-- Compute respective surface fluxes for momentum and TKE 388 CALL calc_surface_fluxes 389 ENDIF 390 ENDDO 391 ! 392 !-- For natural-type surfaces. Please note, even though stability is not considered for the 393 !-- calculation of momentum fluxes at vertical surfaces, scaling parameters and Obukhov length are 394 !-- calculated nevertheless in this case. This is due to the requirement of ts in parameterization 395 !-- of heat flux in land-surface model in case that aero_resist_kray is not true. 396 IF ( .NOT. aero_resist_kray ) THEN 397 DO l = 0, 1 398 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 258 399 surf => surf_lsm_v(l) 259 CALL calc_pt_q260 ENDIF261 262 !-- Call for urban-type vertical surfaces263 IF ( surf_usm_v(l)%ns >= 1 ) THEN264 surf => surf_usm_v(l)265 CALL calc_pt_q266 ENDIF267 ENDDO268 269 !270 !-- First, calculate the new Obukhov length, then new friction velocity,271 !-- followed by the new scaling parameters (th*, q*, etc.), and the new272 !-- surface fluxes if required. Note, each routine is called for different surface types.273 !-- First call for default-type horizontal surfaces, for natural- and274 !-- urban-type surfaces. Note, at this place only upward-facing horizontal275 !-- surfaces are treated.276 277 !278 !-- Default-type upward-facing horizontal surfaces279 IF ( surf_def_h(0)%ns >= 1 ) THEN280 surf => surf_def_h(0)281 CALL calc_uvw_abs282 IF ( .NOT. neutral ) CALL calc_ol283 CALL calc_us284 CALL calc_scaling_parameters285 CALL calc_surface_fluxes286 ENDIF287 !288 !-- Natural-type horizontal surfaces289 IF ( surf_lsm_h%ns >= 1 ) THEN290 surf => surf_lsm_h291 CALL calc_uvw_abs292 IF ( .NOT. neutral ) CALL calc_ol293 CALL calc_us294 CALL calc_scaling_parameters295 CALL calc_surface_fluxes296 ENDIF297 !298 !-- Urban-type horizontal surfaces299 IF ( surf_usm_h%ns >= 1 ) THEN300 surf => surf_usm_h301 CALL calc_uvw_abs302 IF ( .NOT. neutral ) CALL calc_ol303 CALL calc_us304 CALL calc_scaling_parameters305 CALL calc_surface_fluxes306 !307 !-- Calculate 10cm temperature, required in indoor model308 IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' )309 ENDIF310 311 !312 !-- Treat downward-facing horizontal surfaces. Note, so far, these are313 !-- always default type. Stratification is not considered314 !-- in this case, hence, no further distinction between different315 !-- most_method is required.316 IF ( surf_def_h(1)%ns >= 1 ) THEN317 downward = .TRUE.318 surf => surf_def_h(1)319 CALL calc_uvw_abs320 CALL calc_us321 CALL calc_surface_fluxes322 downward = .FALSE.323 ENDIF324 !325 !-- Calculate surfaces fluxes at vertical surfaces for momentum326 !-- and subgrid-scale TKE.327 !-- No stability is considered. Therefore, scaling parameters and Obukhov-328 !-- length do not need to be calculated and no distinction in 'circular',329 !-- 'Newton' or 'lookup' is necessary so far.330 !-- Note, this will change if stability is once considered.331 surf_vertical = .TRUE.332 !333 !-- Calculate horizontal momentum fluxes at north- and south-facing334 !-- surfaces(usvs).335 !-- For default-type surfaces336 mom_uv = .TRUE.337 DO l = 0, 1338 IF ( surf_def_v(l)%ns >= 1 ) THEN339 surf => surf_def_v(l)340 400 ! 341 401 !-- Compute surface-parallel velocity 342 402 CALL calc_uvw_abs_v_ugrid 343 403 ! 404 !-- Compute Obukhov length 405 IF ( .NOT. neutral ) CALL calc_ol 406 ! 344 407 !-- Compute respective friction velocity on staggered grid 345 408 CALL calc_us 409 ! 410 !-- Compute scaling parameters 411 CALL calc_scaling_parameters 346 412 ! 347 413 !-- Compute respective surface fluxes for momentum and TKE … … 350 416 ENDDO 351 417 ! 352 !-- For natural-type surfaces. Please note, even though stability is not 353 !-- considered for the calculation of momentum fluxes at vertical surfaces, 354 !-- scaling parameters and Obukov length are calculated nevertheless in this 355 !-- case. This is due to the requirement of ts in parameterization of heat 356 !-- flux in land-surface model in case of aero_resist_kray is not true. 357 IF ( .NOT. aero_resist_kray ) THEN 358 DO l = 0, 1 359 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 360 surf => surf_lsm_v(l) 361 ! 362 !-- Compute surface-parallel velocity 363 CALL calc_uvw_abs_v_ugrid 364 ! 365 !-- Compute Obukhov length 366 IF ( .NOT. neutral ) CALL calc_ol 367 ! 368 !-- Compute respective friction velocity on staggered grid 369 CALL calc_us 370 ! 371 !-- Compute scaling parameters 372 CALL calc_scaling_parameters 373 ! 374 !-- Compute respective surface fluxes for momentum and TKE 375 CALL calc_surface_fluxes 376 ENDIF 377 ENDDO 378 ! 379 !-- No ts is required, so scaling parameters and Obukhov length do not need 380 !-- to be computed. 381 ELSE 382 DO l = 0, 1 383 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 384 surf => surf_lsm_v(l) 385 ! 386 !-- Compute surface-parallel velocity 387 CALL calc_uvw_abs_v_ugrid 388 ! 389 !-- Compute respective friction velocity on staggered grid 390 CALL calc_us 391 ! 392 !-- Compute respective surface fluxes for momentum and TKE 393 CALL calc_surface_fluxes 394 ENDIF 395 ENDDO 396 ENDIF 397 ! 398 !-- For urban-type surfaces 418 !-- No ts is required, so scaling parameters and Obukhov length do not need to be computed. 419 ELSE 399 420 DO l = 0, 1 400 IF ( surf_ usm_v(l)%ns >= 1) THEN401 surf => surf_ usm_v(l)421 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 422 surf => surf_lsm_v(l) 402 423 ! 403 424 !-- Compute surface-parallel velocity … … 409 430 !-- Compute respective surface fluxes for momentum and TKE 410 431 CALL calc_surface_fluxes 411 !412 !-- Calculate 10cm temperature, required in indoor model413 IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' )414 432 ENDIF 415 433 ENDDO 416 ! 417 !-- Calculate horizontal momentum fluxes at east- and west-facing 418 !-- surfaces (vsus). 419 !-- For default-type surfaces 434 ENDIF 435 ! 436 !-- For urban-type surfaces 437 DO l = 0, 1 438 IF ( surf_usm_v(l)%ns >= 1 ) THEN 439 surf => surf_usm_v(l) 440 ! 441 !-- Compute surface-parallel velocity 442 CALL calc_uvw_abs_v_ugrid 443 ! 444 !-- Compute respective friction velocity on staggered grid 445 CALL calc_us 446 ! 447 !-- Compute respective surface fluxes for momentum and TKE 448 CALL calc_surface_fluxes 449 ! 450 !-- Calculate 10cm temperature, required in indoor model 451 IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' ) 452 ENDIF 453 ENDDO 454 ! 455 !-- Calculate horizontal momentum fluxes at east- and west-facing surfaces (vsus). 456 !-- For default-type surfaces 457 DO l = 2, 3 458 IF ( surf_def_v(l)%ns >= 1 ) THEN 459 surf => surf_def_v(l) 460 ! 461 !-- Compute surface-parallel velocity 462 CALL calc_uvw_abs_v_vgrid 463 ! 464 !-- Compute respective friction velocity on staggered grid 465 CALL calc_us 466 ! 467 !-- Compute respective surface fluxes for momentum and TKE 468 CALL calc_surface_fluxes 469 ENDIF 470 ENDDO 471 ! 472 !-- For natural-type surfaces. Please note, even though stability is not considered for the 473 !-- calculation of momentum fluxes at vertical surfaces, scaling parameters and Obukov length are 474 !-- calculated nevertheless in this case. This is due to the requirement of ts in parameterization 475 !-- of heat flux in land-surface model in case that aero_resist_kray is not true. 476 IF ( .NOT. aero_resist_kray ) THEN 420 477 DO l = 2, 3 421 IF ( surf_ def_v(l)%ns >= 1) THEN422 surf => surf_ def_v(l)478 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 479 surf => surf_lsm_v(l) 423 480 ! 424 481 !-- Compute surface-parallel velocity 425 482 CALL calc_uvw_abs_v_vgrid 426 483 ! 484 !-- Compute Obukhov length 485 IF ( .NOT. neutral ) CALL calc_ol 486 ! 427 487 !-- Compute respective friction velocity on staggered grid 428 488 CALL calc_us 429 489 ! 490 !-- Compute scaling parameters 491 CALL calc_scaling_parameters 492 ! 430 493 !-- Compute respective surface fluxes for momentum and TKE 431 494 CALL calc_surface_fluxes 432 433 495 ENDIF 434 496 ENDDO 435 ! 436 !-- For natural-type surfaces. Please note, even though stability is not 437 !-- considered for the calculation of momentum fluxes at vertical surfaces, 438 !-- scaling parameters and Obukov length are calculated nevertheless in this 439 !-- case. This is due to the requirement of ts in parameterization of heat 440 !-- flux in land-surface model in case of aero_resist_kray is not true. 441 IF ( .NOT. aero_resist_kray ) THEN 442 DO l = 2, 3 443 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 444 surf => surf_lsm_v(l) 445 ! 446 !-- Compute surface-parallel velocity 447 CALL calc_uvw_abs_v_vgrid 448 ! 449 !-- Compute Obukhov length 450 IF ( .NOT. neutral ) CALL calc_ol 451 ! 452 !-- Compute respective friction velocity on staggered grid 453 CALL calc_us 454 ! 455 !-- Compute scaling parameters 456 CALL calc_scaling_parameters 457 ! 458 !-- Compute respective surface fluxes for momentum and TKE 459 CALL calc_surface_fluxes 460 ENDIF 461 ENDDO 462 ELSE 463 DO l = 2, 3 464 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 465 surf => surf_lsm_v(l) 466 ! 467 !-- Compute surface-parallel velocity 468 CALL calc_uvw_abs_v_vgrid 469 ! 470 !-- Compute respective friction velocity on staggered grid 471 CALL calc_us 472 ! 473 !-- Compute respective surface fluxes for momentum and TKE 474 CALL calc_surface_fluxes 475 ENDIF 476 ENDDO 477 ENDIF 478 ! 479 !-- For urban-type surfaces 497 ELSE 480 498 DO l = 2, 3 481 IF ( surf_ usm_v(l)%ns >= 1) THEN482 surf => surf_ usm_v(l)499 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 500 surf => surf_lsm_v(l) 483 501 ! 484 502 !-- Compute surface-parallel velocity … … 490 508 !-- Compute respective surface fluxes for momentum and TKE 491 509 CALL calc_surface_fluxes 492 !493 !-- Calculate 10cm temperature, required in indoor model494 IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' )495 510 ENDIF 496 511 ENDDO 497 mom_uv = .FALSE. 498 ! 499 !-- Calculate horizontal momentum fluxes of w (wsus and wsvs) at vertial 500 !-- surfaces. 501 mom_w = .TRUE. 502 ! 503 !-- Default-type surfaces 504 DO l = 0, 3 505 IF ( surf_def_v(l)%ns >= 1 ) THEN 506 surf => surf_def_v(l) 507 CALL calc_uvw_abs_v_wgrid 508 CALL calc_us 509 CALL calc_surface_fluxes 510 ENDIF 511 ENDDO 512 ! 513 !-- Natural-type surfaces 514 DO l = 0, 3 515 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 516 surf => surf_lsm_v(l) 517 CALL calc_uvw_abs_v_wgrid 518 CALL calc_us 519 CALL calc_surface_fluxes 520 ENDIF 521 ENDDO 522 ! 523 !-- Urban-type surfaces 524 DO l = 0, 3 525 IF ( surf_usm_v(l)%ns >= 1 ) THEN 526 surf => surf_usm_v(l) 527 CALL calc_uvw_abs_v_wgrid 528 CALL calc_us 529 CALL calc_surface_fluxes 530 ENDIF 531 ENDDO 532 mom_w = .FALSE. 533 ! 534 !-- Calculate momentum fluxes usvs, vsus, wsus and wsvs at vertical 535 !-- surfaces for TKE production. Note, here, momentum fluxes are defined 536 !-- at grid center and are not staggered as before. 537 mom_tke = .TRUE. 538 ! 539 !-- Default-type surfaces 540 DO l = 0, 3 541 IF ( surf_def_v(l)%ns >= 1 ) THEN 542 surf => surf_def_v(l) 543 CALL calc_uvw_abs_v_sgrid 544 CALL calc_us 545 CALL calc_surface_fluxes 546 ENDIF 547 ENDDO 548 ! 549 !-- Natural-type surfaces 550 DO l = 0, 3 551 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 552 surf => surf_lsm_v(l) 553 CALL calc_uvw_abs_v_sgrid 554 CALL calc_us 555 CALL calc_surface_fluxes 556 ENDIF 557 ENDDO 558 ! 559 !-- Urban-type surfaces 560 DO l = 0, 3 561 IF ( surf_usm_v(l)%ns >= 1 ) THEN 562 surf => surf_usm_v(l) 563 CALL calc_uvw_abs_v_sgrid 564 CALL calc_us 565 CALL calc_surface_fluxes 566 ENDIF 567 ENDDO 568 mom_tke = .FALSE. 569 570 IF ( debug_output_timestep ) CALL debug_message( 'surface_layer_fluxes', 'end' ) 571 572 END SUBROUTINE surface_layer_fluxes 573 574 575 !------------------------------------------------------------------------------! 512 ENDIF 513 ! 514 !-- For urban-type surfaces 515 DO l = 2, 3 516 IF ( surf_usm_v(l)%ns >= 1 ) THEN 517 surf => surf_usm_v(l) 518 ! 519 !-- Compute surface-parallel velocity 520 CALL calc_uvw_abs_v_vgrid 521 ! 522 !-- Compute respective friction velocity on staggered grid 523 CALL calc_us 524 ! 525 !-- Compute respective surface fluxes for momentum and TKE 526 CALL calc_surface_fluxes 527 ! 528 !-- Calculate 10cm temperature, required in indoor model 529 IF ( indoor_model ) CALL calc_pt_near_surface ( '10cm' ) 530 ENDIF 531 ENDDO 532 mom_uv = .FALSE. 533 ! 534 !-- Calculate horizontal momentum fluxes of w (wsus and wsvs) at vertial surfaces. 535 mom_w = .TRUE. 536 ! 537 !-- Default-type surfaces 538 DO l = 0, 3 539 IF ( surf_def_v(l)%ns >= 1 ) THEN 540 surf => surf_def_v(l) 541 CALL calc_uvw_abs_v_wgrid 542 CALL calc_us 543 CALL calc_surface_fluxes 544 ENDIF 545 ENDDO 546 ! 547 !-- Natural-type surfaces 548 DO l = 0, 3 549 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 550 surf => surf_lsm_v(l) 551 CALL calc_uvw_abs_v_wgrid 552 CALL calc_us 553 CALL calc_surface_fluxes 554 ENDIF 555 ENDDO 556 ! 557 !-- Urban-type surfaces 558 DO l = 0, 3 559 IF ( surf_usm_v(l)%ns >= 1 ) THEN 560 surf => surf_usm_v(l) 561 CALL calc_uvw_abs_v_wgrid 562 CALL calc_us 563 CALL calc_surface_fluxes 564 ENDIF 565 ENDDO 566 mom_w = .FALSE. 567 ! 568 !-- Calculate momentum fluxes usvs, vsus, wsus and wsvs at vertical surfaces for TKE production. 569 !-- Note, here, momentum fluxes are defined at grid center and are not staggered as before. 570 mom_tke = .TRUE. 571 ! 572 !-- Default-type surfaces 573 DO l = 0, 3 574 IF ( surf_def_v(l)%ns >= 1 ) THEN 575 surf => surf_def_v(l) 576 CALL calc_uvw_abs_v_sgrid 577 CALL calc_us 578 CALL calc_surface_fluxes 579 ENDIF 580 ENDDO 581 ! 582 !-- Natural-type surfaces 583 DO l = 0, 3 584 IF ( surf_lsm_v(l)%ns >= 1 ) THEN 585 surf => surf_lsm_v(l) 586 CALL calc_uvw_abs_v_sgrid 587 CALL calc_us 588 CALL calc_surface_fluxes 589 ENDIF 590 ENDDO 591 ! 592 !-- Urban-type surfaces 593 DO l = 0, 3 594 IF ( surf_usm_v(l)%ns >= 1 ) THEN 595 surf => surf_usm_v(l) 596 CALL calc_uvw_abs_v_sgrid 597 CALL calc_us 598 CALL calc_surface_fluxes 599 ENDIF 600 ENDDO 601 mom_tke = .FALSE. 602 603 IF ( debug_output_timestep ) CALL debug_message( 'surface_layer_fluxes', 'end' ) 604 605 END SUBROUTINE surface_layer_fluxes 606 607 608 !--------------------------------------------------------------------------------------------------! 576 609 ! Description: 577 610 ! ------------ 578 611 !> Initializing actions for the surface layer routine. 579 !------------------------------------------------------------------------------! 580 SUBROUTINE init_surface_layer_fluxes 581 582 IMPLICIT NONE 583 584 INTEGER(iwp) :: l !< running index for vertical surface orientation 585 586 CALL location_message( 'initializing surface layer', 'start' ) 587 588 ! 589 !-- In case of runs with neutral statification, set Obukhov length to a 590 !-- large value 591 IF ( neutral ) THEN 592 IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%ol = 1.0E10_wp 593 IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%ol = 1.0E10_wp 594 IF ( surf_usm_h%ns >= 1 ) surf_usm_h%ol = 1.0E10_wp 595 596 DO l = 0, 3 597 IF ( surf_def_v(l)%ns >= 1 .AND. & 598 ALLOCATED( surf_def_v(l)%ol ) ) surf_def_v(l)%ol = 1.0E10_wp 599 IF ( surf_lsm_v(l)%ns >= 1 .AND. & 600 ALLOCATED( surf_lsm_v(l)%ol ) ) surf_lsm_v(l)%ol = 1.0E10_wp 601 IF ( surf_usm_v(l)%ns >= 1 .AND. & 602 ALLOCATED( surf_usm_v(l)%ol ) ) surf_usm_v(l)%ol = 1.0E10_wp 603 ENDDO 604 605 ENDIF 606 607 CALL location_message( 'initializing surface layer', 'finished' ) 608 609 END SUBROUTINE init_surface_layer_fluxes 610 611 612 !------------------------------------------------------------------------------! 612 !--------------------------------------------------------------------------------------------------! 613 SUBROUTINE init_surface_layer_fluxes 614 615 IMPLICIT NONE 616 617 INTEGER(iwp) :: l !< running index for vertical surface orientation 618 619 CALL location_message( 'initializing surface layer', 'start' ) 620 621 ! 622 !-- In case of runs with neutral statification, set Obukhov length to a large value 623 IF ( neutral ) THEN 624 IF ( surf_def_h(0)%ns >= 1 ) surf_def_h(0)%ol = 1.0E10_wp 625 IF ( surf_lsm_h%ns >= 1 ) surf_lsm_h%ol = 1.0E10_wp 626 IF ( surf_usm_h%ns >= 1 ) surf_usm_h%ol = 1.0E10_wp 627 628 DO l = 0, 3 629 IF ( surf_def_v(l)%ns >= 1 .AND. & 630 ALLOCATED( surf_def_v(l)%ol ) ) surf_def_v(l)%ol = 1.0E10_wp 631 IF ( surf_lsm_v(l)%ns >= 1 .AND. & 632 ALLOCATED( surf_lsm_v(l)%ol ) ) surf_lsm_v(l)%ol = 1.0E10_wp 633 IF ( surf_usm_v(l)%ns >= 1 .AND. & 634 ALLOCATED( surf_usm_v(l)%ol ) ) surf_usm_v(l)%ol = 1.0E10_wp 635 ENDDO 636 637 ENDIF 638 639 CALL location_message( 'initializing surface layer', 'finished' ) 640 641 END SUBROUTINE init_surface_layer_fluxes 642 643 644 !--------------------------------------------------------------------------------------------------! 613 645 ! Description: 614 646 ! ------------ 615 !> Compute the absolute value of the horizontal velocity (relative to the 616 !> surface) for horizontal surface elements. This is required by all methods. 617 !------------------------------------------------------------------------------! 618 SUBROUTINE calc_uvw_abs 619 620 IMPLICIT NONE 621 622 INTEGER(iwp) :: i !< running index x direction 623 INTEGER(iwp) :: ibit !< flag to mask computation of relative velocity in case of downward-facing surfaces 624 INTEGER(iwp) :: j !< running index y direction 625 INTEGER(iwp) :: k !< running index z direction 626 INTEGER(iwp) :: m !< running index surface elements 627 628 REAL(wp) :: w_lfc !< local free convection velocity scale 629 ! 630 !-- ibit is 1 for upward-facing surfaces, zero for downward-facing surfaces. 631 ibit = MERGE( 1, 0, .NOT. downward ) 632 633 !$OMP PARALLEL DO PRIVATE(i, j, k, w_lfc) 634 !$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) & 635 !$ACC PRESENT(surf, u, v) 647 !> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal 648 !> surface elements. This is required by all methods. 649 !--------------------------------------------------------------------------------------------------! 650 SUBROUTINE calc_uvw_abs 651 652 IMPLICIT NONE 653 654 INTEGER(iwp) :: i !< running index x direction 655 INTEGER(iwp) :: ibit !< flag to mask computation of relative velocity in case of downward-facing surfaces 656 INTEGER(iwp) :: j !< running index y direction 657 INTEGER(iwp) :: k !< running index z direction 658 INTEGER(iwp) :: m !< running index surface elements 659 660 REAL(wp) :: w_lfc !< local free convection velocity scale 661 ! 662 !-- ibit is 1 for upward-facing surfaces, zero for downward-facing surfaces. 663 ibit = MERGE( 1, 0, .NOT. downward ) 664 665 !$OMP PARALLEL DO PRIVATE(i, j, k, w_lfc) 666 !$ACC PARALLEL LOOP PRIVATE(i, j, k, w_lfc) & 667 !$ACC PRESENT(surf, u, v) 668 DO m = 1, surf%ns 669 i = surf%i(m) 670 j = surf%j(m) 671 k = surf%k(m) 672 673 ! 674 !-- Calculate free convection velocity scale w_lfc is use_free_convection_scaling = .T.. This 675 !-- will maintain a horizontal velocity even for very weak wind convective conditions. SIGN is 676 !-- used to set w_lfc to zero under stable conditions. 677 IF ( use_free_convection_scaling ) THEN 678 w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m)) 679 w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp) 680 ELSE 681 w_lfc = 0.0_wp 682 ENDIF 683 684 ! 685 !-- Compute the absolute value of the horizontal velocity. (relative to the surface in case the 686 !-- lower surface is the ocean). Please note, in new surface modelling concept the index values 687 !-- changed, i.e. the reference grid point is not the surface-grid point itself but the first 688 !-- grid point outside of the topography. Note, in case of coupled ocean-atmosphere simulations 689 !-- relative velocity with respect to the ocean surface is used, hence, (k-1,j,i) values are used 690 !-- to calculate the absolute velocity. However, this does not apply for downward-facing walls. 691 !-- To mask this, use ibit, which checks for upward/downward-facing surfaces. 692 surf%uvw_abs(m) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) - ( u(k-1,j,i) + u(k-1,j,i+1) ) & 693 * ibit ) )**2 & 694 + ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) - ( v(k-1,j,i) + v(k-1,j+1,i) & 695 ) * ibit ) & 696 )**2 + w_lfc**2 ) 697 ENDDO 698 699 END SUBROUTINE calc_uvw_abs 700 701 702 !--------------------------------------------------------------------------------------------------! 703 ! Description: 704 ! ------------ 705 !> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal 706 !> surface elements. This is required by all methods. 707 !--------------------------------------------------------------------------------------------------! 708 SUBROUTINE calc_uvw_abs_v_ugrid 709 710 IMPLICIT NONE 711 712 INTEGER(iwp) :: i !< running index x direction 713 INTEGER(iwp) :: j !< running index y direction 714 INTEGER(iwp) :: k !< running index z direction 715 INTEGER(iwp) :: m !< running index surface elements 716 717 REAL(wp) :: u_i !< u-component on xu-grid 718 REAL(wp) :: w_i !< w-component on xu-grid 719 720 721 DO m = 1, surf%ns 722 i = surf%i(m) 723 j = surf%j(m) 724 k = surf%k(m) 725 ! 726 !-- Compute the absolute value of the surface parallel velocity on u-grid. 727 u_i = u(k,j,i) 728 w_i = 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) 729 730 surf%uvw_abs(m) = SQRT( u_i**2 + w_i**2 ) 731 ENDDO 732 733 END SUBROUTINE calc_uvw_abs_v_ugrid 734 735 !--------------------------------------------------------------------------------------------------! 736 ! Description: 737 ! ------------ 738 !> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal 739 !> surface elements. This is required by all methods. 740 !--------------------------------------------------------------------------------------------------! 741 SUBROUTINE calc_uvw_abs_v_vgrid 742 743 IMPLICIT NONE 744 745 INTEGER(iwp) :: i !< running index x direction 746 INTEGER(iwp) :: j !< running index y direction 747 INTEGER(iwp) :: k !< running index z direction 748 INTEGER(iwp) :: m !< running index surface elements 749 750 REAL(wp) :: v_i !< v-component on yv-grid 751 REAL(wp) :: w_i !< w-component on yv-grid 752 753 754 DO m = 1, surf%ns 755 i = surf%i(m) 756 j = surf%j(m) 757 k = surf%k(m) 758 759 v_i = u(k,j,i) 760 w_i = 0.25_wp * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) 761 762 surf%uvw_abs(m) = SQRT( v_i**2 + w_i**2 ) 763 ENDDO 764 765 END SUBROUTINE calc_uvw_abs_v_vgrid 766 767 !--------------------------------------------------------------------------------------------------! 768 ! Description: 769 ! ------------ 770 !> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal 771 !> surface elements. This is required by all methods. 772 !--------------------------------------------------------------------------------------------------! 773 SUBROUTINE calc_uvw_abs_v_wgrid 774 775 IMPLICIT NONE 776 777 INTEGER(iwp) :: i !< running index x direction 778 INTEGER(iwp) :: j !< running index y direction 779 INTEGER(iwp) :: k !< running index z direction 780 INTEGER(iwp) :: m !< running index surface elements 781 782 REAL(wp) :: u_i !< u-component on x-zw-grid 783 REAL(wp) :: v_i !< v-component on y-zw-grid 784 REAL(wp) :: w_i !< w-component on zw-grid 785 ! 786 !-- North- (l=0) and south-facing (l=1) surfaces 787 IF ( l == 0 .OR. l == 1 ) THEN 636 788 DO m = 1, surf%ns 637 638 i = surf%i(m) 789 i = surf%i(m) 639 790 j = surf%j(m) 640 791 k = surf%k(m) 641 642 ! 643 !-- Calculate free convection velocity scale w_lfc is 644 !-- use_free_convection_scaling = .T.. This will maintain a horizontal 645 !-- velocity even for very weak wind convective conditions. SIGN is used 646 !-- to set w_lfc to zero under stable conditions. 647 IF ( use_free_convection_scaling ) THEN 648 w_lfc = ABS(g / surf%pt1(m) * surf%z_mo(m) * surf%shf(m)) 649 w_lfc = ( 0.5_wp * ( w_lfc + SIGN(w_lfc,surf%shf(m)) ) )**(0.33333_wp) 650 ELSE 651 w_lfc = 0.0_wp 652 ENDIF 653 654 ! 655 !-- Compute the absolute value of the horizontal velocity. 656 !-- (relative to the surface in case the lower surface is the ocean). 657 !-- Please note, in new surface modelling concept the index values changed, 658 !-- i.e. the reference grid point is not the surface-grid point itself but 659 !-- the first grid point outside of the topography. 660 !-- Note, in case of coupled ocean-atmosphere simulations relative velocity 661 !-- with respect to the ocean surface is used, hence, (k-1,j,i) values 662 !-- are used to calculate the absolute velocity. 663 !-- However, this do not apply for downward-facing walls. To mask this, 664 !-- use ibit, which checks for upward/downward-facing surfaces. 665 surf%uvw_abs(m) = SQRT( & 666 ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) & 667 - ( u(k-1,j,i) + u(k-1,j,i+1) & 668 ) * ibit & 669 ) & 670 )**2 + & 671 ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) & 672 - ( v(k-1,j,i) + v(k-1,j+1,i) & 673 ) * ibit & 674 ) & 675 )**2 + w_lfc**2 & 676 ) 677 678 679 792 793 u_i = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) 794 v_i = 0.0_wp 795 w_i = w(k,j,i) 796 797 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 680 798 ENDDO 681 682 END SUBROUTINE calc_uvw_abs 683 684 685 !------------------------------------------------------------------------------! 799 ! 800 !-- East- (l=2) and west-facing (l=3) surfaces 801 ELSE 802 DO m = 1, surf%ns 803 i = surf%i(m) 804 j = surf%j(m) 805 k = surf%k(m) 806 807 u_i = 0.0_wp 808 v_i = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) 809 w_i = w(k,j,i) 810 811 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 812 ENDDO 813 ENDIF 814 815 END SUBROUTINE calc_uvw_abs_v_wgrid 816 817 !--------------------------------------------------------------------------------------------------! 686 818 ! Description: 687 819 ! ------------ 688 !> Compute the absolute value of the horizontal velocity (relative to the 689 !> surface) for horizontal surface elements. This is required by all methods. 690 !------------------------------------------------------------------------------! 691 SUBROUTINE calc_uvw_abs_v_ugrid 692 693 IMPLICIT NONE 694 695 INTEGER(iwp) :: i !< running index x direction 696 INTEGER(iwp) :: j !< running index y direction 697 INTEGER(iwp) :: k !< running index z direction 698 INTEGER(iwp) :: m !< running index surface elements 699 700 REAL(wp) :: u_i !< u-component on xu-grid 701 REAL(wp) :: w_i !< w-component on xu-grid 702 703 820 !> Compute the absolute value of the horizontal velocity (relative to the surface) for horizontal 821 !> surface elements. This is required by all methods. 822 !--------------------------------------------------------------------------------------------------! 823 SUBROUTINE calc_uvw_abs_v_sgrid 824 825 IMPLICIT NONE 826 827 INTEGER(iwp) :: i !< running index x direction 828 INTEGER(iwp) :: j !< running index y direction 829 INTEGER(iwp) :: k !< running index z direction 830 INTEGER(iwp) :: m !< running index surface elements 831 832 REAL(wp) :: u_i !< u-component on scalar grid 833 REAL(wp) :: v_i !< v-component on scalar grid 834 REAL(wp) :: w_i !< w-component on scalar grid 835 836 ! 837 !-- North- (l=0) and south-facing (l=1) walls 838 IF ( l == 0 .OR. l == 1 ) THEN 704 839 DO m = 1, surf%ns 705 i = surf%i(m) 840 i = surf%i(m) 706 841 j = surf%j(m) 707 842 k = surf%k(m) 708 ! 709 !-- Compute the absolute value of the surface parallel velocity on u-grid. 710 u_i = u(k,j,i) 711 w_i = 0.25_wp * ( w(k-1,j,i-1) + w(k-1,j,i) + & 712 w(k,j,i-1) + w(k,j,i) ) 713 714 surf%uvw_abs(m) = SQRT( u_i**2 + w_i**2 ) 715 843 844 u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) 845 v_i = 0.0_wp 846 w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) 847 848 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 716 849 ENDDO 717 718 END SUBROUTINE calc_uvw_abs_v_ugrid 719 720 !------------------------------------------------------------------------------! 721 ! Description: 722 ! ------------ 723 !> Compute the absolute value of the horizontal velocity (relative to the 724 !> surface) for horizontal surface elements. This is required by all methods. 725 !------------------------------------------------------------------------------! 726 SUBROUTINE calc_uvw_abs_v_vgrid 727 728 IMPLICIT NONE 729 730 INTEGER(iwp) :: i !< running index x direction 731 INTEGER(iwp) :: j !< running index y direction 732 INTEGER(iwp) :: k !< running index z direction 733 INTEGER(iwp) :: m !< running index surface elements 734 735 REAL(wp) :: v_i !< v-component on yv-grid 736 REAL(wp) :: w_i !< w-component on yv-grid 737 738 850 ! 851 !-- East- (l=2) and west-facing (l=3) walls 852 ELSE 739 853 DO m = 1, surf%ns 740 i = surf%i(m) 854 i = surf%i(m) 741 855 j = surf%j(m) 742 856 k = surf%k(m) 743 857 744 v_i = u(k,j,i) 745 w_i = 0.25_wp * ( w(k-1,j-1,i) + w(k-1,j,i) + & 746 w(k,j-1,i) + w(k,j,i) ) 747 748 surf%uvw_abs(m) = SQRT( v_i**2 + w_i**2 ) 749 858 u_i = 0.0_wp 859 v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) 860 w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) 861 862 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 750 863 ENDDO 751 752 END SUBROUTINE calc_uvw_abs_v_vgrid 753 754 !------------------------------------------------------------------------------! 755 ! Description: 756 ! ------------ 757 !> Compute the absolute value of the horizontal velocity (relative to the 758 !> surface) for horizontal surface elements. This is required by all methods. 759 !------------------------------------------------------------------------------! 760 SUBROUTINE calc_uvw_abs_v_wgrid 761 762 IMPLICIT NONE 763 764 INTEGER(iwp) :: i !< running index x direction 765 INTEGER(iwp) :: j !< running index y direction 766 INTEGER(iwp) :: k !< running index z direction 767 INTEGER(iwp) :: m !< running index surface elements 768 769 REAL(wp) :: u_i !< u-component on x-zw-grid 770 REAL(wp) :: v_i !< v-component on y-zw-grid 771 REAL(wp) :: w_i !< w-component on zw-grid 772 ! 773 !-- North- (l=0) and south-facing (l=1) surfaces 774 IF ( l == 0 .OR. l == 1 ) THEN 775 DO m = 1, surf%ns 776 i = surf%i(m) 777 j = surf%j(m) 778 k = surf%k(m) 779 780 u_i = 0.25_wp * ( u(k+1,j,i+1) + u(k+1,j,i) + & 781 u(k,j,i+1) + u(k,j,i) ) 782 v_i = 0.0_wp 783 w_i = w(k,j,i) 784 785 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 786 ENDDO 787 ! 788 !-- East- (l=2) and west-facing (l=3) surfaces 789 ELSE 790 DO m = 1, surf%ns 791 i = surf%i(m) 792 j = surf%j(m) 793 k = surf%k(m) 794 795 u_i = 0.0_wp 796 v_i = 0.25_wp * ( v(k+1,j+1,i) + v(k+1,j,i) + & 797 v(k,j+1,i) + v(k,j,i) ) 798 w_i = w(k,j,i) 799 800 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 801 ENDDO 802 ENDIF 803 804 END SUBROUTINE calc_uvw_abs_v_wgrid 805 806 !------------------------------------------------------------------------------! 807 ! Description: 808 ! ------------ 809 !> Compute the absolute value of the horizontal velocity (relative to the 810 !> surface) for horizontal surface elements. This is required by all methods. 811 !------------------------------------------------------------------------------! 812 SUBROUTINE calc_uvw_abs_v_sgrid 813 814 IMPLICIT NONE 815 816 INTEGER(iwp) :: i !< running index x direction 817 INTEGER(iwp) :: j !< running index y direction 818 INTEGER(iwp) :: k !< running index z direction 819 INTEGER(iwp) :: m !< running index surface elements 820 821 REAL(wp) :: u_i !< u-component on scalar grid 822 REAL(wp) :: v_i !< v-component on scalar grid 823 REAL(wp) :: w_i !< w-component on scalar grid 824 825 ! 826 !-- North- (l=0) and south-facing (l=1) walls 827 IF ( l == 0 .OR. l == 1 ) THEN 828 DO m = 1, surf%ns 829 i = surf%i(m) 830 j = surf%j(m) 831 k = surf%k(m) 832 833 u_i = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) 834 v_i = 0.0_wp 835 w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) 836 837 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 838 ENDDO 839 ! 840 !-- East- (l=2) and west-facing (l=3) walls 841 ELSE 842 DO m = 1, surf%ns 843 i = surf%i(m) 844 j = surf%j(m) 845 k = surf%k(m) 846 847 u_i = 0.0_wp 848 v_i = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) 849 w_i = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) 850 851 surf%uvw_abs(m) = SQRT( u_i**2 + v_i**2 + w_i**2 ) 852 ENDDO 853 ENDIF 854 855 END SUBROUTINE calc_uvw_abs_v_sgrid 856 857 858 !------------------------------------------------------------------------------! 864 ENDIF 865 866 END SUBROUTINE calc_uvw_abs_v_sgrid 867 868 869 !--------------------------------------------------------------------------------------------------! 859 870 ! Description: 860 871 ! ------------ 861 872 !> Calculate the Obukhov length (L) and Richardson flux number (z/L) 862 !------------------------------------------------------------------------------! 863 SUBROUTINE calc_ol 864 865 IMPLICIT NONE 866 867 INTEGER(iwp) :: iter !< Newton iteration step 868 INTEGER(iwp) :: m !< loop variable over all horizontal wall elements 869 870 LOGICAL, DIMENSION(surf%ns) :: convergence_reached !< convergence switch for vectorization 871 !$ACC DECLARE CREATE( convergence_reached ) 872 873 REAL(wp) :: f, & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0 874 f_d_ol, & !< Derivative of f 875 ol_l, & !< Lower bound of L for Newton iteration 876 ol_m, & !< Previous value of L for Newton iteration 877 ol_old, & !< Previous time step value of L 878 ol_u !< Upper bound of L for Newton iteration 879 880 REAL(wp), DIMENSION(surf%ns) :: ol_old_vec !< temporary array required for vectorization 881 REAL(wp), DIMENSION(surf%ns) :: z_mo_vec !< temporary array required for vectorization 882 !$ACC DECLARE CREATE( ol_old_vec, z_mo_vec ) 883 884 ! 885 !-- Evaluate bulk Richardson number (calculation depends on 886 !-- definition based on setting of boundary conditions 887 IF ( ibc_pt_b /= 1 ) THEN 888 IF ( humidity ) THEN 889 !$OMP PARALLEL DO PRIVATE( z_mo ) 890 DO m = 1, surf%ns 891 892 z_mo = surf%z_mo(m) 893 894 surf%rib(m) = g * z_mo & 895 * ( surf%vpt1(m) - surf%vpt_surface(m) ) & 896 / ( surf%uvw_abs(m)**2 * surf%vpt1(m) & 897 + 1.0E-20_wp ) 898 ENDDO 899 ELSE 900 !$OMP PARALLEL DO PRIVATE( z_mo ) 901 DO m = 1, surf%ns 902 903 z_mo = surf%z_mo(m) 904 905 surf%rib(m) = g * z_mo & 906 * ( surf%pt1(m) - surf%pt_surface(m) ) & 907 / ( surf%uvw_abs(m)**2 * surf%pt1(m) + 1.0E-20_wp ) 908 ENDDO 873 !--------------------------------------------------------------------------------------------------! 874 SUBROUTINE calc_ol 875 876 IMPLICIT NONE 877 878 INTEGER(iwp) :: iter !< Newton iteration step 879 INTEGER(iwp) :: m !< loop variable over all horizontal wall elements 880 881 LOGICAL, DIMENSION(surf%ns) :: convergence_reached !< convergence switch for vectorization 882 !$ACC DECLARE CREATE( convergence_reached ) 883 884 REAL(wp) :: f !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0 885 REAL(wp) :: f_d_ol !< Derivative of f 886 REAL(wp) :: ol_l !< Lower bound of L for Newton iteration 887 REAL(wp) :: ol_m !< Previous value of L for Newton iteration 888 REAL(wp) :: ol_old !< Previous time step value of L 889 REAL(wp) :: ol_u !< Upper bound of L for Newton iteration 890 891 REAL(wp), DIMENSION(surf%ns) :: ol_old_vec !< temporary array required for vectorization 892 REAL(wp), DIMENSION(surf%ns) :: z_mo_vec !< temporary array required for vectorization 893 !$ACC DECLARE CREATE( ol_old_vec, z_mo_vec ) 894 895 ! 896 !-- Evaluate bulk Richardson number (calculation depends on definition based on setting of boundary 897 !-- conditions) 898 IF ( ibc_pt_b /= 1 ) THEN 899 IF ( humidity ) THEN 900 !$OMP PARALLEL DO PRIVATE( z_mo ) 901 DO m = 1, surf%ns 902 z_mo = surf%z_mo(m) 903 surf%rib(m) = g * z_mo * ( surf%vpt1(m) - surf%vpt_surface(m) ) / & 904 ( surf%uvw_abs(m)**2 * surf%vpt1(m) + 1.0E-20_wp ) 905 ENDDO 906 ELSE 907 !$OMP PARALLEL DO PRIVATE( z_mo ) 908 DO m = 1, surf%ns 909 z_mo = surf%z_mo(m) 910 surf%rib(m) = g * z_mo * ( surf%pt1(m) - surf%pt_surface(m) ) / & 911 ( surf%uvw_abs(m)**2 * surf%pt1(m) + 1.0E-20_wp ) 912 ENDDO 913 ENDIF 914 ELSE 915 IF ( humidity ) THEN 916 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 917 DO m = 1, surf%ns 918 k = surf%k(m) 919 z_mo = surf%z_mo(m) 920 surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp * surf%qv1(m) ) * & 921 surf%shf(m) + 0.61_wp * surf%pt1(m) * surf%qsws(m) ) * & 922 drho_air_zw(k-1) / ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2 & 923 + 1.0E-20_wp ) 924 ENDDO 925 ELSE 926 !$OMP PARALLEL DO PRIVATE( k, z_mo ) 927 !$ACC PARALLEL LOOP PRIVATE(k, z_mo) & 928 !$ACC PRESENT(surf, drho_air_zw) 929 DO m = 1, surf%ns 930 k = surf%k(m) 931 z_mo = surf%z_mo(m) 932 surf%rib(m) = - g * z_mo * surf%shf(m) * drho_air_zw(k-1) / & 933 ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2 + 1.0E-20_wp ) 934 ENDDO 935 ENDIF 936 ENDIF 937 938 IF ( loop_optimization == 'cache' ) THEN 939 ! 940 !-- Calculate the Obukhov length using Newton iteration 941 !$OMP PARALLEL DO PRIVATE(i, j, z_mo) & 942 !$OMP PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) 943 !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) & 944 !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) & 945 !$ACC PRESENT(surf) 946 DO m = 1, surf%ns 947 i = surf%i(m) 948 j = surf%j(m) 949 950 z_mo = surf%z_mo(m) 951 952 ! 953 !-- Store current value in case the Newton iteration fails 954 ol_old = surf%ol(m) 955 956 ! 957 !-- Ensure that the bulk Richardson number and the Obukhov length have the same sign 958 IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. ABS( surf%ol(m) ) == ol_max ) THEN 959 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp 960 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp 909 961 ENDIF 910 ELSE911 IF ( humidity ) THEN912 !$OMP PARALLEL DO PRIVATE( k, z_mo )913 DO m = 1, surf%ns914 915 k = surf%k(m)916 917 z_mo = surf%z_mo(m)918 919 surf%rib(m) = - g * z_mo * ( ( 1.0_wp + 0.61_wp &920 * surf%qv1(m) ) * surf%shf(m) + 0.61_wp &921 * surf%pt1(m) * surf%qsws(m) ) * &922 drho_air_zw(k-1) / &923 ( surf%uvw_abs(m)**3 * surf%vpt1(m) * kappa**2 &924 + 1.0E-20_wp )925 ENDDO926 ELSE927 !$OMP PARALLEL DO PRIVATE( k, z_mo )928 !$ACC PARALLEL LOOP PRIVATE(k, z_mo) &929 !$ACC PRESENT(surf, drho_air_zw)930 DO m = 1, surf%ns931 932 k = surf%k(m)933 934 z_mo = surf%z_mo(m)935 936 surf%rib(m) = - g * z_mo * surf%shf(m) * &937 drho_air_zw(k-1) / &938 ( surf%uvw_abs(m)**3 * surf%pt1(m) * kappa**2 &939 + 1.0E-20_wp )940 ENDDO941 ENDIF942 ENDIF943 944 IF ( loop_optimization == 'cache' ) THEN945 !946 !-- Calculate the Obukhov length using Newton iteration947 !$OMP PARALLEL DO PRIVATE(i, j, z_mo) &948 !$OMP PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol)949 !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &950 !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &951 !$ACC PRESENT(surf)952 DO m = 1, surf%ns953 954 i = surf%i(m)955 j = surf%j(m)956 957 z_mo = surf%z_mo(m)958 959 !960 !-- Store current value in case the Newton iteration fails961 ol_old = surf%ol(m)962 963 !964 !-- Ensure that the bulk Richardson number and the Obukhov965 !-- length have the same sign966 IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. ABS( surf%ol(m) ) == ol_max ) THEN967 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp968 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp969 ENDIF970 !971 !-- Iteration to find Obukhov length972 iter = 0973 DO974 iter = iter + 1975 !976 !-- In case of divergence, use the value of the previous time step977 IF ( iter > 1000 ) THEN978 surf%ol(m) = ol_old979 EXIT980 ENDIF981 982 ol_m = surf%ol(m)983 ol_l = ol_m - 0.001_wp * ol_m984 ol_u = ol_m + 0.001_wp * ol_m985 986 987 IF ( ibc_pt_b /= 1 ) THEN988 !989 !-- Calculate f = Ri - [...]/[...]^2 = 0990 f = surf%rib(m) - ( z_mo / ol_m ) * ( LOG( z_mo / surf%z0h(m) ) &991 - psi_h( z_mo / ol_m ) &992 + psi_h( surf%z0h(m) / ol_m ) &993 ) / &994 ( LOG( z_mo / surf%z0(m) ) &995 - psi_m( z_mo / ol_m ) &996 + psi_m( surf%z0(m) / ol_m ) &997 )**2998 999 !1000 !-- Calculate df/dL1001 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / surf%z0h(m) ) &1002 - psi_h( z_mo / ol_u ) &1003 + psi_h( surf%z0h(m) / ol_u ) &1004 ) / &1005 ( LOG( z_mo / surf%z0(m) ) &1006 - psi_m( z_mo / ol_u ) &1007 + psi_m( surf%z0(m) / ol_u ) &1008 )**2 &1009 + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) &1010 - psi_h( z_mo / ol_l ) &1011 + psi_h( surf%z0h(m) / ol_l ) &1012 ) / &1013 ( LOG( z_mo / surf%z0(m) ) &1014 - psi_m( z_mo / ol_l ) &1015 + psi_m( surf%z0(m) / ol_l ) &1016 )**2 &1017 ) / ( ol_u - ol_l )1018 ELSE1019 !1020 !-- Calculate f = Ri - 1 /[...]^3 = 01021 f = surf%rib(m) - ( z_mo / ol_m ) / ( LOG( z_mo / surf%z0(m) ) &1022 - psi_m( z_mo / ol_m ) &1023 + psi_m( surf%z0(m) / ol_m ) &1024 )**31025 1026 !1027 !-- Calculate df/dL1028 f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) ) &1029 - psi_m( z_mo / ol_u ) &1030 + psi_m( surf%z0(m) / ol_u ) &1031 )**3 &1032 + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) ) &1033 - psi_m( z_mo / ol_l ) &1034 + psi_m( surf%z0(m) / ol_l ) &1035 )**3 &1036 ) / ( ol_u - ol_l )1037 ENDIF1038 !1039 !-- Calculate new L1040 surf%ol(m) = ol_m - f / f_d_ol1041 1042 !1043 !-- Ensure that the bulk Richardson number and the Obukhov1044 !-- length have the same sign and ensure convergence.1045 IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp1046 !1047 !-- If unrealistic value occurs, set L to the maximum1048 !-- value that is allowed1049 IF ( ABS( surf%ol(m) ) > ol_max ) THEN1050 surf%ol(m) = ol_max1051 EXIT1052 ENDIF1053 !1054 !-- Assure that Obukhov length does not become zero. If the limit is1055 !-- reached, exit the loop.1056 IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN1057 surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) )1058 EXIT1059 ENDIF1060 !1061 !-- Check for convergence1062 IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) EXIT1063 1064 ENDDO1065 ENDDO1066 1067 !1068 !-- Vector Version1069 ELSE1070 !1071 !-- Calculate the Obukhov length using Newton iteration1072 !-- First set arrays required for vectorization1073 !$ACC PARALLEL LOOP &1074 !$ACC PRESENT(surf)1075 DO m = 1, surf%ns1076 1077 z_mo_vec(m) = surf%z_mo(m)1078 1079 !1080 !-- Store current value in case the Newton iteration fails1081 ol_old_vec(m) = surf%ol(m)1082 1083 !1084 !-- Ensure that the bulk Richardson number and the Obukhov length have the same sign1085 IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. ABS( surf%ol(m) ) == ol_max ) THEN1086 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp1087 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp1088 ENDIF1089 !1090 !-- Initialize convergence flag1091 convergence_reached(m) = .FALSE.1092 1093 ENDDO1094 1095 962 ! 1096 963 !-- Iteration to find Obukhov length 1097 964 iter = 0 1098 965 DO 1099 1100 966 iter = iter + 1 1101 967 ! 1102 !-- In case of divergence, use the value (s)of the previous time step968 !-- In case of divergence, use the value of the previous time step 1103 969 IF ( iter > 1000 ) THEN 1104 !$ACC PARALLEL LOOP & 1105 !$ACC PRESENT(surf) 1106 DO m = 1, surf%ns 1107 IF ( .NOT. convergence_reached(m) ) surf%ol(m) = ol_old_vec(m) 1108 ENDDO 970 surf%ol(m) = ol_old 1109 971 EXIT 1110 972 ENDIF 1111 973 1112 !$ACC PARALLEL LOOP PRIVATE(ol_m, ol_l, ol_u, f, f_d_ol) & 1113 !$ACC PRESENT(surf) 1114 DO m = 1, surf%ns 1115 1116 IF ( convergence_reached(m) ) CYCLE 1117 1118 ol_m = surf%ol(m) 1119 ol_l = ol_m - 0.001_wp * ol_m 1120 ol_u = ol_m + 0.001_wp * ol_m 1121 1122 1123 IF ( ibc_pt_b /= 1 ) THEN 1124 ! 1125 !-- Calculate f = Ri - [...]/[...]^2 = 0 1126 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( LOG( z_mo_vec(m) / surf%z0h(m) ) & 1127 - psi_h( z_mo_vec(m) / ol_m ) & 1128 + psi_h( surf%z0h(m) / ol_m ) & 1129 ) / & 1130 ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1131 - psi_m( z_mo_vec(m) / ol_m ) & 1132 + psi_m( surf%z0(m) / ol_m ) & 1133 )**2 1134 1135 ! 1136 !-- Calculate df/dL 1137 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( LOG( z_mo_vec(m) / surf%z0h(m) ) & 1138 - psi_h( z_mo_vec(m) / ol_u ) & 1139 + psi_h( surf%z0h(m) / ol_u ) & 1140 ) / & 1141 ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1142 - psi_m( z_mo_vec(m) / ol_u ) & 1143 + psi_m( surf%z0(m) / ol_u ) & 1144 )**2 & 1145 + ( z_mo_vec(m) / ol_l ) * ( LOG( z_mo_vec(m) / surf%z0h(m) ) & 1146 - psi_h( z_mo_vec(m) / ol_l ) & 1147 + psi_h( surf%z0h(m) / ol_l ) & 1148 ) / & 1149 ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1150 - psi_m( z_mo_vec(m) / ol_l ) & 1151 + psi_m( surf%z0(m) / ol_l ) & 1152 )**2 & 1153 ) / ( ol_u - ol_l ) 1154 ELSE 1155 ! 1156 !-- Calculate f = Ri - 1 /[...]^3 = 0 1157 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1158 - psi_m( z_mo_vec(m) / ol_m ) & 1159 + psi_m( surf%z0(m) / ol_m ) & 1160 )**3 1161 1162 ! 1163 !-- Calculate df/dL 1164 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1165 - psi_m( z_mo_vec(m) / ol_u ) & 1166 + psi_m( surf%z0(m) / ol_u ) & 1167 )**3 & 1168 + ( z_mo_vec(m) / ol_l ) / ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1169 - psi_m( z_mo_vec(m) / ol_l ) & 1170 + psi_m( surf%z0(m) / ol_l ) & 1171 )**3 & 1172 ) / ( ol_u - ol_l ) 1173 ENDIF 1174 ! 1175 !-- Calculate new L 1176 surf%ol(m) = ol_m - f / f_d_ol 1177 1178 ! 1179 !-- Ensure that the bulk Richardson number and the Obukhov 1180 !-- length have the same sign and ensure convergence. 1181 IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp 1182 1183 ! 1184 !-- Check for convergence 1185 !-- This check does not modify surf%ol, therefore this is done first 1186 IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) THEN 1187 convergence_reached(m) = .TRUE. 1188 ENDIF 1189 ! 1190 !-- If unrealistic value occurs, set L to the maximum allowed value 1191 IF ( ABS( surf%ol(m) ) > ol_max ) THEN 1192 surf%ol(m) = ol_max 1193 convergence_reached(m) = .TRUE. 1194 ENDIF 1195 1196 ENDDO 1197 ! 1198 !-- Assure that Obukhov length does not become zero 974 ol_m = surf%ol(m) 975 ol_l = ol_m - 0.001_wp * ol_m 976 ol_u = ol_m + 0.001_wp * ol_m 977 978 979 IF ( ibc_pt_b /= 1 ) THEN 980 ! 981 !-- Calculate f = Ri - [...]/[...]^2 = 0 982 f = surf%rib(m) - ( z_mo / ol_m ) * ( LOG( z_mo / surf%z0h(m) ) & 983 - psi_h( z_mo / ol_m ) & 984 + psi_h( surf%z0h(m) / ol_m ) ) / & 985 ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_m ) & 986 + psi_m( surf%z0(m) / ol_m ) )**2 987 988 ! 989 !-- Calculate df/dL 990 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / surf%z0h(m) ) & 991 - psi_h( z_mo / ol_u ) & 992 + psi_h( surf%z0h(m) / ol_u ) ) / & 993 ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_u ) & 994 + psi_m( surf%z0(m) / ol_u ) )**2 & 995 + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / ol_l ) & 996 + psi_h( surf%z0h(m) / ol_l ) ) / & 997 ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_l ) & 998 + psi_m( surf%z0(m) / ol_l ) )**2 ) / ( ol_u - ol_l ) 999 ELSE 1000 ! 1001 !-- Calculate f = Ri - 1 /[...]^3 = 0 1002 f = surf%rib(m) - ( z_mo / ol_m ) / & 1003 ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / ol_m ) + psi_m( surf%z0(m) / ol_m ) & 1004 )**3 1005 1006 ! 1007 !-- Calculate df/dL 1008 f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) ) & 1009 - psi_m( z_mo / ol_u ) & 1010 + psi_m( surf%z0(m) / ol_u ) & 1011 )**3 & 1012 + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) ) & 1013 - psi_m( z_mo / ol_l ) & 1014 + psi_m( surf%z0(m) / ol_l ) & 1015 )**3 & 1016 ) / ( ol_u - ol_l ) 1017 ENDIF 1018 ! 1019 !-- Calculate new L 1020 surf%ol(m) = ol_m - f / f_d_ol 1021 1022 ! 1023 !-- Ensure that the bulk Richardson number and the Obukhov length have the same sign and 1024 !-- ensure convergence. 1025 IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp 1026 ! 1027 !-- If unrealistic value occurs, set L to the maximum value that is allowed 1028 IF ( ABS( surf%ol(m) ) > ol_max ) THEN 1029 surf%ol(m) = ol_max 1030 EXIT 1031 ENDIF 1032 ! 1033 !-- Assure that Obukhov length does not become zero. If the limit is reached, exit the loop. 1034 IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN 1035 surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) ) 1036 EXIT 1037 ENDIF 1038 ! 1039 !-- Check for convergence 1040 IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) EXIT 1041 ENDDO 1042 ENDDO 1043 1044 ! 1045 !-- Vector Version 1046 ELSE 1047 ! 1048 !-- Calculate the Obukhov length using Newton iteration 1049 !-- First set arrays required for vectorization 1050 !$ACC PARALLEL LOOP & 1051 !$ACC PRESENT(surf) 1052 DO m = 1, surf%ns 1053 z_mo_vec(m) = surf%z_mo(m) 1054 ! 1055 !-- Store current value in case the Newton iteration fails 1056 ol_old_vec(m) = surf%ol(m) 1057 ! 1058 !-- Ensure that the bulk Richardson number and the Obukhov length have the same sign 1059 IF ( surf%rib(m) * surf%ol(m) < 0.0_wp .OR. ABS( surf%ol(m) ) == ol_max ) THEN 1060 IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) = 0.01_wp 1061 IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp 1062 ENDIF 1063 ! 1064 !-- Initialize convergence flag 1065 convergence_reached(m) = .FALSE. 1066 ENDDO 1067 1068 ! 1069 !-- Iteration to find Obukhov length 1070 iter = 0 1071 DO 1072 iter = iter + 1 1073 ! 1074 !-- In case of divergence, use the value(s) of the previous time step 1075 IF ( iter > 1000 ) THEN 1199 1076 !$ACC PARALLEL LOOP & 1200 1077 !$ACC PRESENT(surf) 1201 1078 DO m = 1, surf%ns 1202 IF ( convergence_reached(m) ) CYCLE 1203 IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN 1204 surf%ol(m) = SIGN( 10E-6_wp, surf%ol(m) ) 1205 convergence_reached(m) = .TRUE. 1206 ENDIF 1079 IF ( .NOT. convergence_reached(m) ) surf%ol(m) = ol_old_vec(m) 1207 1080 ENDDO 1208 1209 IF ( ALL( convergence_reached ) ) EXIT 1210 1211 ENDDO ! end of iteration loop 1212 1213 ENDIF ! end of vector branch 1214 1215 END SUBROUTINE calc_ol 1216 1217 ! 1218 !-- Calculate friction velocity u* 1219 SUBROUTINE calc_us 1220 1221 IMPLICIT NONE 1222 1223 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1224 1225 ! 1226 !-- Compute u* at horizontal surfaces at the scalars' grid points 1227 IF ( .NOT. surf_vertical ) THEN 1228 ! 1229 !-- Compute u* at upward-facing surfaces 1230 IF ( .NOT. downward ) THEN 1231 !$OMP PARALLEL DO PRIVATE( z_mo ) 1232 !$ACC PARALLEL LOOP PRIVATE(z_mo) & 1233 !$ACC PRESENT(surf) 1234 DO m = 1, surf%ns 1235 1236 z_mo = surf%z_mo(m) 1237 ! 1238 !-- Compute u* at the scalars' grid points 1239 surf%us(m) = kappa * surf%uvw_abs(m) / & 1240 ( LOG( z_mo / surf%z0(m) ) & 1241 - psi_m( z_mo / surf%ol(m) ) & 1242 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1243 1244 ENDDO 1245 ! 1246 !-- Compute u* at downward-facing surfaces. This case, do not consider 1247 !-- any stability. 1248 ELSE 1249 !$OMP PARALLEL DO PRIVATE( z_mo ) 1250 !$ACC PARALLEL LOOP PRIVATE(z_mo) & 1251 !$ACC PRESENT(surf) 1252 DO m = 1, surf%ns 1253 1254 z_mo = surf%z_mo(m) 1255 ! 1256 !-- Compute u* at the scalars' grid points 1257 surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) ) 1258 1259 ENDDO 1081 EXIT 1260 1082 ENDIF 1261 ! 1262 !-- Compute u* at vertical surfaces at the u/v/v grid, respectively. 1263 !-- No stability is considered in this case. 1264 ELSE 1265 !$OMP PARALLEL DO PRIVATE( z_mo ) 1083 1084 !$ACC PARALLEL LOOP PRIVATE(ol_m, ol_l, ol_u, f, f_d_ol) & 1085 !$ACC PRESENT(surf) 1086 DO m = 1, surf%ns 1087 IF ( convergence_reached(m) ) CYCLE 1088 1089 ol_m = surf%ol(m) 1090 ol_l = ol_m - 0.001_wp * ol_m 1091 ol_u = ol_m + 0.001_wp * ol_m 1092 1093 1094 IF ( ibc_pt_b /= 1 ) THEN 1095 ! 1096 !-- Calculate f = Ri - [...]/[...]^2 = 0 1097 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( LOG( z_mo_vec(m) / surf%z0h(m) ) & 1098 - psi_h( z_mo_vec(m) / ol_m ) & 1099 + psi_h( surf%z0h(m) / ol_m ) & 1100 ) / & 1101 ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1102 - psi_m( z_mo_vec(m) / ol_m ) & 1103 + psi_m( surf%z0(m) / ol_m ) & 1104 )**2 1105 1106 ! 1107 !-- Calculate df/dL 1108 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( LOG( z_mo_vec(m) / surf%z0h(m) ) & 1109 - psi_h( z_mo_vec(m) / ol_u ) & 1110 + psi_h( surf%z0h(m) / ol_u ) & 1111 ) / & 1112 ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1113 - psi_m( z_mo_vec(m) / ol_u ) & 1114 + psi_m( surf%z0(m) / ol_u ) & 1115 )**2 & 1116 + ( z_mo_vec(m) / ol_l ) * ( LOG( z_mo_vec(m) / surf%z0h(m) ) & 1117 - psi_h( z_mo_vec(m) / ol_l ) & 1118 + psi_h( surf%z0h(m) / ol_l ) & 1119 ) / & 1120 ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1121 - psi_m( z_mo_vec(m) / ol_l ) & 1122 + psi_m( surf%z0(m) / ol_l ) & 1123 )**2 & 1124 ) / ( ol_u - ol_l ) 1125 ELSE 1126 ! 1127 !-- Calculate f = Ri - 1 /[...]^3 = 0 1128 f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1129 - psi_m( z_mo_vec(m) / ol_m ) & 1130 + psi_m( surf%z0(m) / ol_m ) & 1131 )**3 1132 1133 ! 1134 !-- Calculate df/dL 1135 f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1136 - psi_m( z_mo_vec(m) / ol_u ) & 1137 + psi_m( surf%z0(m) / ol_u ) & 1138 )**3 & 1139 + ( z_mo_vec(m) / ol_l ) / ( LOG( z_mo_vec(m) / surf%z0(m) ) & 1140 - psi_m( z_mo_vec(m) / ol_l ) & 1141 + psi_m( surf%z0(m) / ol_l ) & 1142 )**3 & 1143 ) / ( ol_u - ol_l ) 1144 ENDIF 1145 ! 1146 !-- Calculate new L 1147 surf%ol(m) = ol_m - f / f_d_ol 1148 1149 ! 1150 !-- Ensure that the bulk Richardson number and the Obukhov length have the same sign and 1151 !-- ensure convergence. 1152 IF ( surf%ol(m) * ol_m < 0.0_wp ) surf%ol(m) = ol_m * 0.5_wp 1153 1154 ! 1155 !-- Check for convergence 1156 !-- This check does not modify surf%ol, therefore this is done first 1157 IF ( ABS( ( surf%ol(m) - ol_m ) / surf%ol(m) ) < 1.0E-4_wp ) THEN 1158 convergence_reached(m) = .TRUE. 1159 ENDIF 1160 ! 1161 !-- If unrealistic value occurs, set L to the maximum allowed value 1162 IF ( ABS( surf%ol(m) ) > ol_max ) THEN 1163 surf%ol(m) = ol_max 1164 convergence_reached(m) = .TRUE. 1165 ENDIF 1166 ENDDO 1167 ! 1168 !-- Assure that Obukhov length does not become zero 1169 !$ACC PARALLEL LOOP & 1170 !$ACC PRESENT(surf) 1171 DO m = 1, surf%ns 1172 IF ( convergence_reached(m) ) CYCLE 1173 IF ( ABS( surf%ol(m) ) < 1E-5_wp ) THEN 1174 surf%ol(m) = SIGN( 10E-6_wp, surf%ol(m) ) 1175 convergence_reached(m) = .TRUE. 1176 ENDIF 1177 ENDDO 1178 1179 IF ( ALL( convergence_reached ) ) EXIT 1180 1181 ENDDO ! End of iteration loop 1182 1183 ENDIF ! End of vector branch 1184 1185 END SUBROUTINE calc_ol 1186 1187 1188 !--------------------------------------------------------------------------------------------------! 1189 ! Description: 1190 ! ------------ 1191 !> Calculate friction velocity u*. 1192 !--------------------------------------------------------------------------------------------------! 1193 SUBROUTINE calc_us 1194 1195 IMPLICIT NONE 1196 1197 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1198 1199 ! 1200 !-- Compute u* at horizontal surfaces at the scalars' grid points 1201 IF ( .NOT. surf_vertical ) THEN 1202 ! 1203 !-- Compute u* at upward-facing surfaces 1204 IF ( .NOT. downward ) THEN 1205 !$OMP PARALLEL DO PRIVATE( z_mo ) 1266 1206 !$ACC PARALLEL LOOP PRIVATE(z_mo) & 1267 1207 !$ACC PRESENT(surf) 1268 1208 DO m = 1, surf%ns 1269 1209 z_mo = surf%z_mo(m) 1270 1210 ! 1211 !-- Compute u* at the scalars' grid points 1212 surf%us(m) = kappa * surf%uvw_abs(m) / ( LOG( z_mo / surf%z0(m) ) & 1213 - psi_m( z_mo / surf%ol(m) ) & 1214 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1215 ENDDO 1216 ! 1217 !-- Compute u* at downward-facing surfaces. This case, do not consider any stability. 1218 ELSE 1219 !$OMP PARALLEL DO PRIVATE( z_mo ) 1220 !$ACC PARALLEL LOOP PRIVATE(z_mo) & 1221 !$ACC PRESENT(surf) 1222 DO m = 1, surf%ns 1223 z_mo = surf%z_mo(m) 1224 ! 1225 !-- Compute u* at the scalars' grid points 1271 1226 surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) ) 1272 1227 ENDDO 1273 1228 ENDIF 1274 1275 END SUBROUTINE calc_us 1276 1277 ! 1278 !-- Calculate potential temperature, specific humidity, and virtual potential 1279 !-- temperature at first grid level 1280 SUBROUTINE calc_pt_q 1281 1282 IMPLICIT NONE 1283 1284 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1229 ! 1230 !-- Compute u* at vertical surfaces at the u/v/v grid, respectively. 1231 !-- No stability is considered in this case. 1232 ELSE 1233 !$OMP PARALLEL DO PRIVATE( z_mo ) 1234 !$ACC PARALLEL LOOP PRIVATE(z_mo) & 1235 !$ACC PRESENT(surf) 1236 DO m = 1, surf%ns 1237 z_mo = surf%z_mo(m) 1238 surf%us(m) = kappa * surf%uvw_abs(m) / LOG( z_mo / surf%z0(m) ) 1239 ENDDO 1240 ENDIF 1241 1242 END SUBROUTINE calc_us 1243 1244 !--------------------------------------------------------------------------------------------------! 1245 ! Description: 1246 ! ------------ 1247 !> Calculate potential temperature, specific humidity, and virtual potential temperature at first 1248 !> grid level. 1249 !--------------------------------------------------------------------------------------------------! 1250 SUBROUTINE calc_pt_q 1251 1252 IMPLICIT NONE 1253 1254 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1255 1256 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1257 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 1258 !$ACC PRESENT(surf, pt) 1259 DO m = 1, surf%ns 1260 i = surf%i(m) 1261 j = surf%j(m) 1262 k = surf%k(m) 1263 1264 #ifndef _OPENACC 1265 IF ( bulk_cloud_model ) THEN 1266 surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i) 1267 surf%qv1(m) = q(k,j,i) - ql(k,j,i) 1268 ELSEIF( cloud_droplets ) THEN 1269 surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i) 1270 surf%qv1(m) = q(k,j,i) 1271 ELSE 1272 #endif 1273 surf%pt1(m) = pt(k,j,i) 1274 #ifndef _OPENACC 1275 IF ( humidity ) THEN 1276 surf%qv1(m) = q(k,j,i) 1277 ELSE 1278 #endif 1279 surf%qv1(m) = 0.0_wp 1280 #ifndef _OPENACC 1281 ENDIF 1282 ENDIF 1283 1284 IF ( humidity ) THEN 1285 surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) ) 1286 ENDIF 1287 #endif 1288 ENDDO 1289 1290 END SUBROUTINE calc_pt_q 1291 1292 1293 !--------------------------------------------------------------------------------------------------! 1294 ! Description: 1295 ! ------------ 1296 !> Set potential temperature at surface grid level( only for upward-facing surfs ). 1297 !--------------------------------------------------------------------------------------------------! 1298 SUBROUTINE calc_pt_surface 1299 1300 IMPLICIT NONE 1301 1302 INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 1303 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1304 1305 k_off = surf%koff 1306 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1307 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 1308 !$ACC PRESENT(surf, pt) 1309 DO m = 1, surf%ns 1310 i = surf%i(m) 1311 j = surf%j(m) 1312 k = surf%k(m) 1313 surf%pt_surface(m) = pt(k+k_off,j,i) 1314 ENDDO 1315 1316 END SUBROUTINE calc_pt_surface 1317 1318 ! 1319 !-- Set mixing ratio at surface grid level. ( Only for upward-facing surfs. ) 1320 SUBROUTINE calc_q_surface 1321 1322 IMPLICIT NONE 1323 1324 INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 1325 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1326 1327 k_off = surf%koff 1328 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1329 DO m = 1, surf%ns 1330 i = surf%i(m) 1331 j = surf%j(m) 1332 k = surf%k(m) 1333 surf%q_surface(m) = q(k+k_off,j,i) 1334 ENDDO 1335 1336 END SUBROUTINE calc_q_surface 1337 1338 !--------------------------------------------------------------------------------------------------! 1339 ! Description: 1340 ! ------------ 1341 !> Set virtual potential temperature at surface grid level ( only for upward-facing surfs ). 1342 !--------------------------------------------------------------------------------------------------! 1343 SUBROUTINE calc_vpt_surface 1344 1345 IMPLICIT NONE 1346 1347 INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 1348 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1349 1350 k_off = surf%koff 1351 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1352 DO m = 1, surf%ns 1353 i = surf%i(m) 1354 j = surf%j(m) 1355 k = surf%k(m) 1356 surf%vpt_surface(m) = vpt(k+k_off,j,i) 1357 1358 ENDDO 1359 1360 END SUBROUTINE calc_vpt_surface 1361 1362 !--------------------------------------------------------------------------------------------------! 1363 ! Description: 1364 ! ------------ 1365 !> Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*) 1366 !--------------------------------------------------------------------------------------------------! 1367 SUBROUTINE calc_scaling_parameters 1368 1369 IMPLICIT NONE 1370 1371 1372 INTEGER(iwp) :: lsp !< running index for chemical species 1373 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1374 ! 1375 !-- Compute theta* at horizontal surfaces 1376 IF ( constant_heatflux .AND. .NOT. surf_vertical ) THEN 1377 ! 1378 !-- For a given heat flux in the surface layer: 1285 1379 1286 1380 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1287 1381 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 1288 !$ACC PRESENT(surf, pt) 1289 DO m = 1, surf%ns 1290 1291 i = surf%i(m) 1292 j = surf%j(m) 1293 k = surf%k(m) 1294 1295 #ifndef _OPENACC 1296 IF ( bulk_cloud_model ) THEN 1297 surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i) 1298 surf%qv1(m) = q(k,j,i) - ql(k,j,i) 1299 ELSEIF( cloud_droplets ) THEN 1300 surf%pt1(m) = pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i) 1301 surf%qv1(m) = q(k,j,i) 1302 ELSE 1303 #endif 1304 surf%pt1(m) = pt(k,j,i) 1305 #ifndef _OPENACC 1306 IF ( humidity ) THEN 1307 surf%qv1(m) = q(k,j,i) 1308 ELSE 1309 #endif 1310 surf%qv1(m) = 0.0_wp 1311 #ifndef _OPENACC 1312 ENDIF 1313 ENDIF 1314 1315 IF ( humidity ) THEN 1316 surf%vpt1(m) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) ) 1317 ENDIF 1318 #endif 1319 1382 !$ACC PRESENT(surf, drho_air_zw) 1383 DO m = 1, surf%ns 1384 i = surf%i(m) 1385 j = surf%j(m) 1386 k = surf%k(m) 1387 surf%ts(m) = -surf%shf(m) * drho_air_zw(k-1) / ( surf%us(m) + 1E-30_wp ) 1388 ! 1389 !-- ts must be limited, because otherwise overflow may occur in case of us=0 when computing 1390 !-- ol further below 1391 IF ( surf%ts(m) < -1.05E5_wp ) surf%ts(m) = -1.0E5_wp 1392 IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp 1320 1393 ENDDO 1321 1394 1322 END SUBROUTINE calc_pt_q 1323 1324 1325 ! 1326 !-- Set potential temperature at surface grid level. 1327 !-- ( only for upward-facing surfs ) 1328 SUBROUTINE calc_pt_surface 1329 1330 IMPLICIT NONE 1331 1332 INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 1333 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1334 1335 k_off = surf%koff 1336 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1337 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 1338 !$ACC PRESENT(surf, pt) 1339 DO m = 1, surf%ns 1340 1341 i = surf%i(m) 1342 j = surf%j(m) 1343 k = surf%k(m) 1344 1345 surf%pt_surface(m) = pt(k+k_off,j,i) 1346 1347 ENDDO 1348 1349 END SUBROUTINE calc_pt_surface 1350 1351 ! 1352 !-- Set mixing ratio at surface grid level. ( Only for upward-facing surfs. ) 1353 SUBROUTINE calc_q_surface 1354 1355 IMPLICIT NONE 1356 1357 INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 1358 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1359 1360 k_off = surf%koff 1361 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1362 DO m = 1, surf%ns 1363 1364 i = surf%i(m) 1365 j = surf%j(m) 1366 k = surf%k(m) 1367 1368 surf%q_surface(m) = q(k+k_off,j,i) 1369 1370 ENDDO 1371 1372 END SUBROUTINE calc_q_surface 1373 1374 ! 1375 !-- Set virtual potential temperature at surface grid level. 1376 !-- ( only for upward-facing surfs ) 1377 SUBROUTINE calc_vpt_surface 1378 1379 IMPLICIT NONE 1380 1381 INTEGER(iwp) :: k_off !< index offset between surface and atmosphere grid point (-1 for upward-, +1 for downward-facing walls) 1382 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1383 1384 k_off = surf%koff 1385 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1386 DO m = 1, surf%ns 1387 1388 i = surf%i(m) 1389 j = surf%j(m) 1390 k = surf%k(m) 1391 1392 surf%vpt_surface(m) = vpt(k+k_off,j,i) 1393 1394 ENDDO 1395 1396 END SUBROUTINE calc_vpt_surface 1397 1398 ! 1399 !-- Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*) 1400 SUBROUTINE calc_scaling_parameters 1401 1402 IMPLICIT NONE 1403 1404 1405 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1406 INTEGER(iwp) :: lsp !< running index for chemical species 1407 ! 1408 !-- Compute theta* at horizontal surfaces 1409 IF ( constant_heatflux .AND. .NOT. surf_vertical ) THEN 1410 ! 1411 !-- For a given heat flux in the surface layer: 1395 ELSEIF ( .NOT. surf_vertical ) THEN 1396 ! 1397 !-- For a given surface temperature: 1398 IF ( large_scale_forcing .AND. lsf_surf ) THEN 1412 1399 1413 1400 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1414 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 1415 !$ACC PRESENT(surf, drho_air_zw) 1416 DO m = 1, surf%ns 1417 1418 i = surf%i(m) 1401 DO m = 1, surf%ns 1402 i = surf%i(m) 1419 1403 j = surf%j(m) 1420 1404 k = surf%k(m) 1421 1422 surf%ts(m) = -surf%shf(m) * drho_air_zw(k-1) / & 1423 ( surf%us(m) + 1E-30_wp ) 1424 1425 ! 1426 !-- ts must be limited, because otherwise overflow may occur in case 1427 !-- of us=0 when computing ol further below 1428 IF ( surf%ts(m) < -1.05E5_wp ) surf%ts(m) = -1.0E5_wp 1429 IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp 1430 1431 ENDDO 1432 1433 ELSEIF ( .NOT. surf_vertical ) THEN 1434 ! 1435 !-- For a given surface temperature: 1405 pt(k-1,j,i) = pt_surface 1406 ENDDO 1407 ENDIF 1408 1409 !$OMP PARALLEL DO PRIVATE( z_mo ) 1410 DO m = 1, surf%ns 1411 z_mo = surf%z_mo(m) 1412 surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) ) & 1413 / ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / surf%ol(m) ) & 1414 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1415 ENDDO 1416 1417 ENDIF 1418 ! 1419 !-- Compute theta* at vertical surfaces. This is only required in case of land-surface model, in 1420 !-- order to compute aerodynamical resistance. 1421 IF ( surf_vertical ) THEN 1422 !$OMP PARALLEL DO PRIVATE( i, j ) 1423 DO m = 1, surf%ns 1424 i = surf%i(m) 1425 j = surf%j(m) 1426 surf%ts(m) = -surf%shf(m) / ( surf%us(m) + 1E-30_wp ) 1427 ! 1428 !-- ts must be limited, because otherwise overflow may occur in case of us=0 when computing ol 1429 !-- further below 1430 IF ( surf%ts(m) < -1.05E5_wp ) surf%ts(m) = -1.0E5_wp 1431 IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp 1432 ENDDO 1433 ENDIF 1434 1435 ! 1436 !-- If required compute q* at horizontal surfaces 1437 IF ( humidity ) THEN 1438 IF ( constant_waterflux .AND. .NOT. surf_vertical ) THEN 1439 ! 1440 !-- For a given water flux in the surface layer 1441 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1442 DO m = 1, surf%ns 1443 i = surf%i(m) 1444 j = surf%j(m) 1445 k = surf%k(m) 1446 surf%qs(m) = -surf%qsws(m) * drho_air_zw(k-1) / ( surf%us(m) + 1E-30_wp ) 1447 ENDDO 1448 1449 ELSEIF ( .NOT. surf_vertical ) THEN 1450 coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled ) 1451 1436 1452 IF ( large_scale_forcing .AND. lsf_surf ) THEN 1437 1438 1453 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1439 DO m = 1, surf%ns1440 i = surf%i(m)1441 j = surf%j(m)1442 k = surf%k(m)1443 1444 pt(k-1,j,i) = pt_surface1445 ENDDO1446 ENDIF1447 1448 !$OMP PARALLEL DO PRIVATE( z_mo )1449 DO m = 1, surf%ns1450 1451 z_mo = surf%z_mo(m)1452 1453 surf%ts(m) = kappa * ( surf%pt1(m) - surf%pt_surface(m) ) &1454 / ( LOG( z_mo / surf%z0h(m) ) &1455 - psi_h( z_mo / surf%ol(m) ) &1456 + psi_h( surf%z0h(m) / surf%ol(m) ) )1457 1458 ENDDO1459 1460 ENDIF1461 !1462 !-- Compute theta* at vertical surfaces. This is only required in case of1463 !-- land-surface model, in order to compute aerodynamical resistance.1464 IF ( surf_vertical ) THEN1465 !$OMP PARALLEL DO PRIVATE( i, j )1466 DO m = 1, surf%ns1467 1468 i = surf%i(m)1469 j = surf%j(m)1470 surf%ts(m) = -surf%shf(m) / ( surf%us(m) + 1E-30_wp )1471 !1472 !-- ts must be limited, because otherwise overflow may occur in case1473 !-- of us=0 when computing ol further below1474 IF ( surf%ts(m) < -1.05E5_wp ) surf%ts(m) = -1.0E5_wp1475 IF ( surf%ts(m) > 1.0E5_wp ) surf%ts(m) = 1.0E5_wp1476 1477 ENDDO1478 ENDIF1479 1480 !1481 !-- If required compute q* at horizontal surfaces1482 IF ( humidity ) THEN1483 IF ( constant_waterflux .AND. .NOT. surf_vertical ) THEN1484 !1485 !-- For a given water flux in the surface layer1486 !$OMP PARALLEL DO PRIVATE( i, j, k )1487 DO m = 1, surf%ns1488 1489 i = surf%i(m)1490 j = surf%j(m)1491 k = surf%k(m)1492 surf%qs(m) = -surf%qsws(m) * drho_air_zw(k-1) / &1493 ( surf%us(m) + 1E-30_wp )1494 1495 ENDDO1496 1497 ELSEIF ( .NOT. surf_vertical ) THEN1498 coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. &1499 run_coupled )1500 1501 IF ( large_scale_forcing .AND. lsf_surf ) THEN1502 !$OMP PARALLEL DO PRIVATE( i, j, k )1503 DO m = 1, surf%ns1504 1505 i = surf%i(m)1506 j = surf%j(m)1507 k = surf%k(m)1508 q(k-1,j,i) = q_surface1509 1510 ENDDO1511 ENDIF1512 1513 !1514 !-- Assume saturation for atmosphere coupled to ocean (but not1515 !-- in case of precursor runs)1516 IF ( coupled_run ) THEN1517 !$OMP PARALLEL DO PRIVATE( i, j, k, e_s )1518 DO m = 1, surf%ns1519 i = surf%i(m)1520 j = surf%j(m)1521 k = surf%k(m)1522 e_s = 6.1_wp * &1523 EXP( 0.07_wp * ( MIN(pt(k-1,j,i),pt(k,j,i)) &1524 - 273.15_wp ) )1525 q(k-1,j,i) = rd_d_rv * e_s / ( surface_pressure - e_s )1526 ENDDO1527 ENDIF1528 1529 IF ( bulk_cloud_model .OR. cloud_droplets ) THEN1530 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1531 DO m = 1, surf%ns1532 1533 i = surf%i(m)1534 j = surf%j(m)1535 k = surf%k(m)1536 1537 z_mo = surf%z_mo(m)1538 1539 surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) ) &1540 / ( LOG( z_mo / surf%z0q(m) ) &1541 - psi_h( z_mo / surf%ol(m) ) &1542 + psi_h( surf%z0q(m) / &1543 surf%ol(m) ) )1544 ENDDO1545 ELSE1546 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1547 DO m = 1, surf%ns1548 1549 i = surf%i(m)1550 j = surf%j(m)1551 k = surf%k(m)1552 1553 z_mo = surf%z_mo(m)1554 1555 surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) ) &1556 / ( LOG( z_mo / surf%z0q(m) ) &1557 - psi_h( z_mo / surf%ol(m) ) &1558 + psi_h( surf%z0q(m) / &1559 surf%ol(m) ) )1560 ENDDO1561 ENDIF1562 ENDIF1563 !1564 !-- Compute q* at vertical surfaces1565 IF ( surf_vertical ) THEN1566 !$OMP PARALLEL DO PRIVATE( i, j )1567 DO m = 1, surf%ns1568 1569 i = surf%i(m)1570 j = surf%j(m)1571 surf%qs(m) = -surf%qsws(m) / ( surf%us(m) + 1E-30_wp )1572 1573 ENDDO1574 ENDIF1575 ENDIF1576 1577 !1578 !-- If required compute s*1579 IF ( passive_scalar ) THEN1580 !1581 !-- At horizontal surfaces1582 IF ( constant_scalarflux .AND. .NOT. surf_vertical ) THEN1583 !1584 !-- For a given scalar flux in the surface layer1585 !$OMP PARALLEL DO PRIVATE( i, j )1586 DO m = 1, surf%ns1587 i = surf%i(m)1588 j = surf%j(m)1589 surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )1590 ENDDO1591 ELSEIF ( .NOT. surf_vertical ) THEN1592 1593 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1594 DO m = 1, surf%ns1595 i = surf%i(m)1596 j = surf%j(m)1597 k = surf%k(m)1598 z_mo = surf%z_mo(m)1599 1600 surf%ss(m) = kappa * ( s(k,j,i) - s(k-1,j,i) ) &1601 / ( LOG( z_mo / surf%z0h(m) ) &1602 - psi_h( z_mo / surf%ol(m) ) &1603 + psi_h( surf%z0h(m) / surf%ol(m) ) )1604 1605 ENDDO1606 ENDIF1607 !1608 !-- At vertical surfaces1609 IF ( surf_vertical ) THEN1610 !$OMP PARALLEL DO PRIVATE( i, j )1611 DO m = 1, surf%ns1612 i = surf%i(m)1613 j = surf%j(m)1614 surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp )1615 ENDDO1616 ENDIF1617 ENDIF1618 1619 !1620 !-- If required compute cs* (chemical species)1621 IF ( air_chemistry ) THEN1622 !1623 !-- At horizontal surfaces1624 DO lsp = 1, nvar1625 IF ( constant_csflux(lsp) .AND. .NOT. surf_vertical ) THEN1626 !-- For a given chemical species' flux in the surface layer1627 !$OMP PARALLEL DO PRIVATE( i, j )1628 DO m = 1, surf%ns1629 i = surf%i(m)1630 j = surf%j(m)1631 surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )1632 ENDDO1633 ENDIF1634 ENDDO1635 !1636 !-- At vertical surfaces1637 IF ( surf_vertical ) THEN1638 DO lsp = 1, nvar1639 !$OMP PARALLEL DO PRIVATE( i, j )1640 DO m = 1, surf%ns1641 i = surf%i(m)1642 j = surf%j(m)1643 surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp )1644 ENDDO1645 ENDDO1646 ENDIF1647 ENDIF1648 1649 !1650 !-- If required compute qc* and nc*1651 IF ( bulk_cloud_model .AND. microphysics_morrison .AND. &1652 .NOT. surf_vertical ) THEN1653 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1654 DO m = 1, surf%ns1655 1656 i = surf%i(m)1657 j = surf%j(m)1658 k = surf%k(m)1659 1660 z_mo = surf%z_mo(m)1661 1662 surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) ) &1663 / ( LOG( z_mo / surf%z0q(m) ) &1664 - psi_h( z_mo / surf%ol(m) ) &1665 + psi_h( surf%z0q(m) / surf%ol(m) ) )1666 1667 surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) ) &1668 / ( LOG( z_mo / surf%z0q(m) ) &1669 - psi_h( z_mo / surf%ol(m) ) &1670 + psi_h( surf%z0q(m) / surf%ol(m) ) )1671 ENDDO1672 1673 ENDIF1674 1675 !1676 !-- If required compute qr* and nr*1677 IF ( bulk_cloud_model .AND. microphysics_seifert .AND. &1678 .NOT. surf_vertical ) THEN1679 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1680 DO m = 1, surf%ns1681 1682 i = surf%i(m)1683 j = surf%j(m)1684 k = surf%k(m)1685 1686 z_mo = surf%z_mo(m)1687 1688 surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) ) &1689 / ( LOG( z_mo / surf%z0q(m) ) &1690 - psi_h( z_mo / surf%ol(m) ) &1691 + psi_h( surf%z0q(m) / surf%ol(m) ) )1692 1693 surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) ) &1694 / ( LOG( z_mo / surf%z0q(m) ) &1695 - psi_h( z_mo / surf%ol(m) ) &1696 + psi_h( surf%z0q(m) / surf%ol(m) ) )1697 ENDDO1698 1699 ENDIF1700 1701 END SUBROUTINE calc_scaling_parameters1702 1703 1704 1705 !1706 !-- Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws)1707 SUBROUTINE calc_surface_fluxes1708 1709 IMPLICIT NONE1710 1711 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements1712 INTEGER(iwp) :: lsp !< running index for chemical species1713 1714 REAL(wp) :: dum !< dummy to precalculate logarithm1715 REAL(wp) :: flag_u !< flag indicating u-grid, used for calculation of horizontal momentum fluxes at vertical surfaces1716 REAL(wp) :: flag_v !< flag indicating v-grid, used for calculation of horizontal momentum fluxes at vertical surfaces1717 REAL(wp), DIMENSION(:), ALLOCATABLE :: u_i !< u-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces1718 REAL(wp), DIMENSION(:), ALLOCATABLE :: v_i !< v-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces1719 REAL(wp), DIMENSION(:), ALLOCATABLE :: w_i !< w-component interpolated onto scalar grid point, required for momentum fluxes at vertical surfaces1720 1721 !1722 !-- Calcuate surface fluxes at horizontal walls1723 IF ( .NOT. surf_vertical ) THEN1724 !1725 !-- Compute u'w' for the total model domain at upward-facing surfaces.1726 !-- First compute the corresponding component of u* and square it.1727 IF ( .NOT. downward ) THEN1728 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1729 !$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &1730 !$ACC PRESENT(surf, u, rho_air_zw)1731 DO m = 1, surf%ns1732 1733 i = surf%i(m)1734 j = surf%j(m)1735 k = surf%k(m)1736 1737 z_mo = surf%z_mo(m)1738 1739 surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) ) &1740 / ( LOG( z_mo / surf%z0(m) ) &1741 - psi_m( z_mo / surf%ol(m) ) &1742 + psi_m( surf%z0(m) / surf%ol(m) ) )1743 !1744 !-- Please note, the computation of usws is not fully accurate. Actually1745 !-- a further interpolation of us onto the u-grid, where usws is defined,1746 !-- is required. However, this is not done as this would require several1747 !-- data transfers between 2D-grid and the surf-type.1748 !-- The impact of the missing interpolation is negligible as several1749 !-- tests had shown.1750 !-- Same also for ol.1751 surf%usws(m) = -surf%usws(m) * surf%us(m) * rho_air_zw(k-1)1752 1753 ENDDO1754 !1755 !-- At downward-facing surfaces1756 ELSE1757 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1758 DO m = 1, surf%ns1759 1760 i = surf%i(m)1761 j = surf%j(m)1762 k = surf%k(m)1763 1764 z_mo = surf%z_mo(m)1765 1766 surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) )1767 surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k)1768 1769 ENDDO1770 ENDIF1771 1772 !1773 !-- Compute v'w' for the total model domain.1774 !-- First compute the corresponding component of u* and square it.1775 !-- Upward-facing surfaces1776 IF ( .NOT. downward ) THEN1777 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1778 !$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) &1779 !$ACC PRESENT(surf, v, rho_air_zw)1780 DO m = 1, surf%ns1781 i = surf%i(m)1782 j = surf%j(m)1783 k = surf%k(m)1784 1785 z_mo = surf%z_mo(m)1786 1787 surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) ) &1788 / ( LOG( z_mo / surf%z0(m) ) &1789 - psi_m( z_mo / surf%ol(m) ) &1790 + psi_m( surf%z0(m) / surf%ol(m) ) )1791 !1792 !-- Please note, the computation of vsws is not fully accurate. Actually1793 !-- a further interpolation of us onto the v-grid, where vsws is defined,1794 !-- is required. However, this is not done as this would require several1795 !-- data transfers between 2D-grid and the surf-type.1796 !-- The impact of the missing interpolation is negligible as several1797 !-- tests had shown.1798 !-- Same also for ol.1799 surf%vsws(m) = -surf%vsws(m) * surf%us(m) * rho_air_zw(k-1)1800 ENDDO1801 !1802 !-- Downward-facing surfaces1803 ELSE1804 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1805 DO m = 1, surf%ns1806 i = surf%i(m)1807 j = surf%j(m)1808 k = surf%k(m)1809 1810 z_mo = surf%z_mo(m)1811 1812 surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) )1813 surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k)1814 ENDDO1815 ENDIF1816 !1817 !-- Compute the vertical kinematic heat flux1818 IF ( .NOT. constant_heatflux .AND. ( ( time_since_reference_point&1819 <= skip_time_do_lsm .AND. simulated_time > 0.0_wp ) .OR. &1820 .NOT. land_surface ) .AND. .NOT. urban_surface .AND. &1821 .NOT. downward ) THEN1822 !$OMP PARALLEL DO PRIVATE( i, j, k )1823 DO m = 1, surf%ns1824 i = surf%i(m)1825 j = surf%j(m)1826 k = surf%k(m)1827 surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1)1828 ENDDO1829 ENDIF1830 !1831 !-- Compute the vertical water flux1832 IF ( .NOT. constant_waterflux .AND. humidity .AND. &1833 ( ( time_since_reference_point <= skip_time_do_lsm .AND. &1834 simulated_time > 0.0_wp ) .OR. .NOT. land_surface ) .AND. &1835 .NOT. urban_surface .AND. .NOT. downward ) THEN1836 !$OMP PARALLEL DO PRIVATE( i, j, k )1837 DO m = 1, surf%ns1838 i = surf%i(m)1839 j = surf%j(m)1840 k = surf%k(m)1841 surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1)1842 ENDDO1843 ENDIF1844 !1845 !-- Compute the vertical scalar flux1846 IF ( .NOT. constant_scalarflux .AND. passive_scalar .AND. &1847 .NOT. downward ) THEN1848 !$OMP PARALLEL DO PRIVATE( i, j )1849 DO m = 1, surf%ns1850 1851 i = surf%i(m)1852 j = surf%j(m)1853 surf%ssws(m) = -surf%ss(m) * surf%us(m)1854 1855 ENDDO1856 ENDIF1857 !1858 !-- Compute the vertical chemical species' flux1859 DO lsp = 1, nvar1860 IF ( .NOT. constant_csflux(lsp) .AND. air_chemistry .AND. &1861 .NOT. downward ) THEN1862 !$OMP PARALLEL DO PRIVATE( i, j )1863 DO m = 1, surf%ns1864 i = surf%i(m)1865 j = surf%j(m)1866 surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m)1867 ENDDO1868 ENDIF1869 ENDDO1870 1871 !1872 !-- Compute (turbulent) fluxes of cloud water content and cloud drop conc.1873 IF ( bulk_cloud_model .AND. microphysics_morrison .AND. &1874 .NOT. downward) THEN1875 !$OMP PARALLEL DO PRIVATE( i, j )1876 DO m = 1, surf%ns1877 1878 i = surf%i(m)1879 j = surf%j(m)1880 1881 surf%qcsws(m) = -surf%qcs(m) * surf%us(m)1882 surf%ncsws(m) = -surf%ncs(m) * surf%us(m)1883 ENDDO1884 ENDIF1885 !1886 !-- Compute (turbulent) fluxes of rain water content and rain drop conc.1887 IF ( bulk_cloud_model .AND. microphysics_seifert .AND. &1888 .NOT. downward) THEN1889 !$OMP PARALLEL DO PRIVATE( i, j )1890 DO m = 1, surf%ns1891 1892 i = surf%i(m)1893 j = surf%j(m)1894 1895 surf%qrsws(m) = -surf%qrs(m) * surf%us(m)1896 surf%nrsws(m) = -surf%nrs(m) * surf%us(m)1897 ENDDO1898 ENDIF1899 1900 !1901 !-- Bottom boundary condition for the TKE.1902 IF ( ibc_e_b == 2 ) THEN1903 !$OMP PARALLEL DO PRIVATE( i, j, k )1904 DO m = 1, surf%ns1905 1906 i = surf%i(m)1907 j = surf%j(m)1908 k = surf%k(m)1909 1910 e(k,j,i) = ( surf%us(m) / 0.1_wp )**21911 !1912 !-- As a test: cm = 0.41913 ! e(k,j,i) = ( us(j,i) / 0.4_wp )**21914 e(k-1,j,i) = e(k,j,i)1915 1916 ENDDO1917 ENDIF1918 !1919 !-- Calcuate surface fluxes at vertical surfaces. No stability is considered.1920 ELSE1921 !1922 !-- Compute usvs l={0,1} and vsus l={2,3}1923 IF ( mom_uv ) THEN1924 !1925 !-- Generalize computation by introducing flags. At north- and south-1926 !-- facing surfaces u-component is used, at east- and west-facing1927 !-- surfaces v-component is used.1928 flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0 .OR. l == 1 )1929 flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2 .OR. l == 3 )1930 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1931 DO m = 1, surf%ns1932 i = surf%i(m)1933 j = surf%j(m)1934 k = surf%k(m)1935 1936 z_mo = surf%z_mo(m)1937 1938 surf%mom_flux_uv(m) = kappa * &1939 ( flag_u * u(k,j,i) + flag_v * v(k,j,i) ) / &1940 LOG( z_mo / surf%z0(m) )1941 1942 surf%mom_flux_uv(m) = &1943 - surf%mom_flux_uv(m) * surf%us(m)1944 ENDDO1945 ENDIF1946 !1947 !-- Compute wsus l={0,1} and wsvs l={2,3}1948 IF ( mom_w ) THEN1949 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )1950 DO m = 1, surf%ns1951 i = surf%i(m)1952 j = surf%j(m)1953 k = surf%k(m)1954 1955 z_mo = surf%z_mo(m)1956 1957 surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) )1958 1959 surf%mom_flux_w(m) = &1960 - surf%mom_flux_w(m) * surf%us(m)1961 ENDDO1962 ENDIF1963 !1964 !-- Compute momentum fluxes used for subgrid-scale TKE production at1965 !-- vertical surfaces. In constrast to the calculated momentum fluxes at1966 !-- vertical surfaces before, which are defined on the u/v/w-grid,1967 !-- respectively), the TKE fluxes are defined at the scalar grid.1968 !--1969 IF ( mom_tke ) THEN1970 !1971 !-- Precalculate velocity components at scalar grid point.1972 ALLOCATE( u_i(1:surf%ns) )1973 ALLOCATE( v_i(1:surf%ns) )1974 ALLOCATE( w_i(1:surf%ns) )1975 1976 IF ( l == 0 .OR. l == 1 ) THEN1977 !$OMP PARALLEL DO PRIVATE( i, j, k )1978 DO m = 1, surf%ns1979 i = surf%i(m)1980 j = surf%j(m)1981 k = surf%k(m)1982 1983 u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )1984 v_i(m) = 0.0_wp1985 w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )1986 ENDDO1987 ELSE1988 !$OMP PARALLEL DO PRIVATE( i, j, k )1989 DO m = 1, surf%ns1990 i = surf%i(m)1991 j = surf%j(m)1992 k = surf%k(m)1993 1994 u_i(m) = 0.0_wp1995 v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )1996 w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )1997 ENDDO1998 ENDIF1999 2000 !$OMP PARALLEL DO PRIVATE( i, j, dum, z_mo )2001 DO m = 1, surf%ns2002 i = surf%i(m)2003 j = surf%j(m)2004 2005 z_mo = surf%z_mo(m)2006 2007 dum = kappa / LOG( z_mo / surf%z0(m) )2008 !2009 !-- usvs (l=0,1) and vsus (l=2,3)2010 surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) )2011 !2012 !-- wsvs (l=0,1) and wsus (l=2,3)2013 surf%mom_flux_tke(1,m) = dum * w_i(m)2014 2015 surf%mom_flux_tke(0:1,m) = &2016 - surf%mom_flux_tke(0:1,m) * surf%us(m)2017 ENDDO2018 !2019 !-- Deallocate temporary arrays2020 DEALLOCATE( u_i )2021 DEALLOCATE( v_i )2022 DEALLOCATE( w_i )2023 ENDIF2024 ENDIF2025 2026 END SUBROUTINE calc_surface_fluxes2027 2028 2029 !------------------------------------------------------------------------------!2030 ! Description:2031 ! ------------2032 !> Calculates temperature near surface (10 cm) for indoor model or 2 m2033 !> temperature for output2034 !------------------------------------------------------------------------------!2035 SUBROUTINE calc_pt_near_surface ( z_char )2036 2037 IMPLICIT NONE2038 2039 CHARACTER (LEN = *), INTENT(IN) :: z_char !< string identifier to identify z level2040 INTEGER(iwp) :: i !< grid index x-dimension2041 INTEGER(iwp) :: j !< grid index y-dimension2042 INTEGER(iwp) :: k !< grid index z-dimension2043 INTEGER(iwp) :: m !< running index for surface elements2044 2045 2046 SELECT CASE ( z_char)2047 2048 2049 CASE ( '10cm' )2050 2051 1454 DO m = 1, surf%ns 2052 2053 1455 i = surf%i(m) 2054 1456 j = surf%j(m) 2055 1457 k = surf%k(m) 2056 2057 surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa & 2058 * ( LOG( 0.1_wp / surf%z0h(m) ) & 2059 - psi_h( 0.1_wp / surf%ol(m) ) & 2060 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1458 q(k-1,j,i) = q_surface 2061 1459 2062 1460 ENDDO 2063 2064 END SELECT 2065 2066 END SUBROUTINE calc_pt_near_surface 2067 2068 2069 ! 2070 !-- Integrated stability function for momentum 2071 FUNCTION psi_m( zeta ) 2072 !$ACC ROUTINE SEQ 2073 2074 USE kinds 2075 2076 IMPLICIT NONE 2077 2078 REAL(wp) :: psi_m !< Integrated similarity function result 2079 REAL(wp) :: zeta !< Stability parameter z/L 2080 REAL(wp) :: x !< dummy variable 2081 2082 REAL(wp), PARAMETER :: a = 1.0_wp !< constant 2083 REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant 2084 REAL(wp), PARAMETER :: c = 5.0_wp !< constant 2085 REAL(wp), PARAMETER :: d = 0.35_wp !< constant 2086 REAL(wp), PARAMETER :: c_d_d = c / d !< constant 2087 REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant 2088 2089 2090 IF ( zeta < 0.0_wp ) THEN 2091 x = SQRT( SQRT( 1.0_wp - 16.0_wp * zeta ) ) 2092 psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2 & 2093 * ( 1.0_wp + x**2 ) * 0.125_wp ) 1461 ENDIF 1462 1463 ! 1464 !-- Assume saturation for atmosphere coupled to ocean (but not in case of precursor runs) 1465 IF ( coupled_run ) THEN 1466 !$OMP PARALLEL DO PRIVATE( i, j, k, e_s ) 1467 DO m = 1, surf%ns 1468 i = surf%i(m) 1469 j = surf%j(m) 1470 k = surf%k(m) 1471 e_s = 6.1_wp * EXP( 0.07_wp * ( MIN( pt(k-1,j,i), pt(k,j,i) ) - 273.15_wp ) ) 1472 q(k-1,j,i) = rd_d_rv * e_s / ( surface_pressure - e_s ) 1473 ENDDO 1474 ENDIF 1475 1476 IF ( bulk_cloud_model .OR. cloud_droplets ) THEN 1477 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1478 DO m = 1, surf%ns 1479 i = surf%i(m) 1480 j = surf%j(m) 1481 k = surf%k(m) 1482 z_mo = surf%z_mo(m) 1483 surf%qs(m) = kappa * ( surf%qv1(m) - surf%q_surface(m) ) & 1484 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) ) & 1485 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1486 ENDDO 1487 ELSE 1488 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1489 DO m = 1, surf%ns 1490 i = surf%i(m) 1491 j = surf%j(m) 1492 k = surf%k(m) 1493 z_mo = surf%z_mo(m) 1494 surf%qs(m) = kappa * ( q(k,j,i) - q(k-1,j,i) ) & 1495 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) ) & 1496 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1497 ENDDO 1498 ENDIF 1499 ENDIF 1500 ! 1501 !-- Compute q* at vertical surfaces 1502 IF ( surf_vertical ) THEN 1503 !$OMP PARALLEL DO PRIVATE( i, j ) 1504 DO m = 1, surf%ns 1505 1506 i = surf%i(m) 1507 j = surf%j(m) 1508 surf%qs(m) = -surf%qsws(m) / ( surf%us(m) + 1E-30_wp ) 1509 1510 ENDDO 1511 ENDIF 1512 ENDIF 1513 1514 ! 1515 !-- If required compute s* 1516 IF ( passive_scalar ) THEN 1517 ! 1518 !-- At horizontal surfaces 1519 IF ( constant_scalarflux .AND. .NOT. surf_vertical ) THEN 1520 ! 1521 !-- For a given scalar flux in the surface layer 1522 !$OMP PARALLEL DO PRIVATE( i, j ) 1523 DO m = 1, surf%ns 1524 i = surf%i(m) 1525 j = surf%j(m) 1526 surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp ) 1527 ENDDO 1528 ELSEIF ( .NOT. surf_vertical ) THEN 1529 1530 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1531 DO m = 1, surf%ns 1532 i = surf%i(m) 1533 j = surf%j(m) 1534 k = surf%k(m) 1535 z_mo = surf%z_mo(m) 1536 1537 surf%ss(m) = kappa * ( s(k,j,i) - s(k-1,j,i) ) & 1538 / ( LOG( z_mo / surf%z0h(m) ) - psi_h( z_mo / surf%ol(m) ) & 1539 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1540 ENDDO 1541 ENDIF 1542 ! 1543 !-- At vertical surfaces 1544 IF ( surf_vertical ) THEN 1545 !$OMP PARALLEL DO PRIVATE( i, j ) 1546 DO m = 1, surf%ns 1547 i = surf%i(m) 1548 j = surf%j(m) 1549 surf%ss(m) = -surf%ssws(m) / ( surf%us(m) + 1E-30_wp ) 1550 ENDDO 1551 ENDIF 1552 ENDIF 1553 1554 ! 1555 !-- If required compute cs* (chemical species) 1556 IF ( air_chemistry ) THEN 1557 ! 1558 !-- At horizontal surfaces 1559 DO lsp = 1, nvar 1560 IF ( constant_csflux(lsp) .AND. .NOT. surf_vertical ) THEN 1561 !-- For a given chemical species' flux in the surface layer 1562 !$OMP PARALLEL DO PRIVATE( i, j ) 1563 DO m = 1, surf%ns 1564 i = surf%i(m) 1565 j = surf%j(m) 1566 surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp ) 1567 ENDDO 1568 ENDIF 1569 ENDDO 1570 ! 1571 !-- At vertical surfaces 1572 IF ( surf_vertical ) THEN 1573 DO lsp = 1, nvar 1574 !$OMP PARALLEL DO PRIVATE( i, j ) 1575 DO m = 1, surf%ns 1576 i = surf%i(m) 1577 j = surf%j(m) 1578 surf%css(lsp,m) = -surf%cssws(lsp,m) / ( surf%us(m) + 1E-30_wp ) 1579 ENDDO 1580 ENDDO 1581 ENDIF 1582 ENDIF 1583 1584 ! 1585 !-- If required compute qc* and nc* 1586 IF ( bulk_cloud_model .AND. microphysics_morrison .AND. .NOT. surf_vertical ) THEN 1587 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1588 DO m = 1, surf%ns 1589 i = surf%i(m) 1590 j = surf%j(m) 1591 k = surf%k(m) 1592 1593 z_mo = surf%z_mo(m) 1594 1595 surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) ) & 1596 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) ) & 1597 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1598 1599 surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) ) & 1600 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) ) & 1601 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1602 ENDDO 1603 1604 ENDIF 1605 1606 ! 1607 !-- If required compute qr* and nr* 1608 IF ( bulk_cloud_model .AND. microphysics_seifert .AND. .NOT. surf_vertical ) THEN 1609 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1610 DO m = 1, surf%ns 1611 i = surf%i(m) 1612 j = surf%j(m) 1613 k = surf%k(m) 1614 1615 z_mo = surf%z_mo(m) 1616 1617 surf%qrs(m) = kappa * ( qr(k,j,i) - qr(k-1,j,i) ) & 1618 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) ) & 1619 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1620 1621 surf%nrs(m) = kappa * ( nr(k,j,i) - nr(k-1,j,i) ) & 1622 / ( LOG( z_mo / surf%z0q(m) ) - psi_h( z_mo / surf%ol(m) ) & 1623 + psi_h( surf%z0q(m) / surf%ol(m) ) ) 1624 ENDDO 1625 1626 ENDIF 1627 1628 END SUBROUTINE calc_scaling_parameters 1629 1630 1631 1632 !--------------------------------------------------------------------------------------------------! 1633 ! Description: 1634 ! ------------ 1635 !> Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws) 1636 !--------------------------------------------------------------------------------------------------! 1637 SUBROUTINE calc_surface_fluxes 1638 1639 IMPLICIT NONE 1640 1641 INTEGER(iwp) :: lsp !< running index for chemical species 1642 INTEGER(iwp) :: m !< loop variable over all horizontal surf elements 1643 1644 REAL(wp) :: dum !< dummy to precalculate logarithm 1645 REAL(wp) :: flag_u !< flag indicating u-grid, used for calculation of horizontal momentum fluxes at vertical surfaces 1646 REAL(wp) :: flag_v !< flag indicating v-grid, used for calculation of horizontal momentum fluxes at vertical surfaces 1647 1648 REAL(wp), DIMENSION(:), ALLOCATABLE :: u_i !< u-component interpolated onto scalar grid point, required for momentum fluxes 1649 !< at vertical surfaces 1650 REAL(wp), DIMENSION(:), ALLOCATABLE :: v_i !< v-component interpolated onto scalar grid point, required for momentum fluxes 1651 !< at vertical surfaces 1652 REAL(wp), DIMENSION(:), ALLOCATABLE :: w_i !< w-component interpolated onto scalar grid point, required for momentum fluxes 1653 !< at vertical surfaces 1654 1655 ! 1656 !-- Calcuate surface fluxes at horizontal walls 1657 IF ( .NOT. surf_vertical ) THEN 1658 ! 1659 !-- Compute u'w' for the total model domain at upward-facing surfaces. First compute the 1660 !-- corresponding component of u* and square it. 1661 IF ( .NOT. downward ) THEN 1662 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1663 !$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) & 1664 !$ACC PRESENT(surf, u, rho_air_zw) 1665 DO m = 1, surf%ns 1666 i = surf%i(m) 1667 j = surf%j(m) 1668 k = surf%k(m) 1669 1670 z_mo = surf%z_mo(m) 1671 1672 surf%usws(m) = kappa * ( u(k,j,i) - u(k-1,j,i) ) & 1673 / ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / surf%ol(m) ) & 1674 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1675 ! 1676 !-- Please note, the computation of usws is not fully accurate. Actually a further 1677 !-- interpolation of us onto the u-grid, where usws is defined, is required. However, this 1678 !-- is not done as this would require several data transfers between 2D-grid and the 1679 !-- surf-type. The impact of the missing interpolation is negligible as several tests have 1680 !-- shown. Same also for ol. 1681 surf%usws(m) = -surf%usws(m) * surf%us(m) * rho_air_zw(k-1) 1682 ENDDO 1683 ! 1684 !-- At downward-facing surfaces 2094 1685 ELSE 2095 2096 psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta & 2097 - bc_d_d 2098 ! 2099 !-- Old version for stable conditions (only valid for z/L < 0.5) 2100 !-- psi_m = - 5.0_wp * zeta 2101 2102 ENDIF 2103 2104 END FUNCTION psi_m 2105 2106 2107 ! 2108 !-- Integrated stability function for heat and moisture 2109 FUNCTION psi_h( zeta ) 2110 !$ACC ROUTINE SEQ 2111 2112 USE kinds 2113 2114 IMPLICIT NONE 2115 2116 REAL(wp) :: psi_h !< Integrated similarity function result 2117 REAL(wp) :: zeta !< Stability parameter z/L 2118 REAL(wp) :: x !< dummy variable 2119 2120 REAL(wp), PARAMETER :: a = 1.0_wp !< constant 2121 REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant 2122 REAL(wp), PARAMETER :: c = 5.0_wp !< constant 2123 REAL(wp), PARAMETER :: d = 0.35_wp !< constant 2124 REAL(wp), PARAMETER :: c_d_d = c / d !< constant 2125 REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant 2126 2127 2128 IF ( zeta < 0.0_wp ) THEN 2129 x = SQRT( 1.0_wp - 16.0_wp * zeta ) 2130 psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp ) 1686 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1687 DO m = 1, surf%ns 1688 i = surf%i(m) 1689 j = surf%j(m) 1690 k = surf%k(m) 1691 1692 z_mo = surf%z_mo(m) 1693 1694 surf%usws(m) = kappa * u(k,j,i) / LOG( z_mo / surf%z0(m) ) 1695 surf%usws(m) = surf%usws(m) * surf%us(m) * rho_air_zw(k) 1696 ENDDO 1697 ENDIF 1698 1699 ! 1700 !-- Compute v'w' for the total model domain. First compute the corresponding component of u* and 1701 !-- square it. 1702 !-- Upward-facing surfaces 1703 IF ( .NOT. downward ) THEN 1704 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1705 !$ACC PARALLEL LOOP PRIVATE(i, j, k, z_mo) & 1706 !$ACC PRESENT(surf, v, rho_air_zw) 1707 DO m = 1, surf%ns 1708 i = surf%i(m) 1709 j = surf%j(m) 1710 k = surf%k(m) 1711 1712 z_mo = surf%z_mo(m) 1713 1714 surf%vsws(m) = kappa * ( v(k,j,i) - v(k-1,j,i) ) & 1715 / ( LOG( z_mo / surf%z0(m) ) - psi_m( z_mo / surf%ol(m) ) & 1716 + psi_m( surf%z0(m) / surf%ol(m) ) ) 1717 ! 1718 !-- Please note, the computation of vsws is not fully accurate. Actually a further 1719 !-- interpolation of us onto the v-grid, where vsws is defined, is required. However, this 1720 !-- is not done as this would require several data transfers between 2D-grid and the 1721 !-- surf-type. The impact of the missing interpolation is negligible as several tests have 1722 !-- shown. Same also for ol. 1723 surf%vsws(m) = -surf%vsws(m) * surf%us(m) * rho_air_zw(k-1) 1724 ENDDO 1725 ! 1726 !-- Downward-facing surfaces 2131 1727 ELSE 2132 psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp & 2133 + 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d & 2134 + 1.0_wp 2135 ! 2136 !-- Old version for stable conditions (only valid for z/L < 0.5) 2137 !-- psi_h = - 5.0_wp * zeta 2138 ENDIF 2139 2140 END FUNCTION psi_h 2141 2142 2143 !------------------------------------------------------------------------------! 1728 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1729 DO m = 1, surf%ns 1730 i = surf%i(m) 1731 j = surf%j(m) 1732 k = surf%k(m) 1733 1734 z_mo = surf%z_mo(m) 1735 1736 surf%vsws(m) = kappa * v(k,j,i) / LOG( z_mo / surf%z0(m) ) 1737 surf%vsws(m) = surf%vsws(m) * surf%us(m) * rho_air_zw(k) 1738 ENDDO 1739 ENDIF 1740 ! 1741 !-- Compute the vertical kinematic heat flux 1742 IF ( .NOT. constant_heatflux .AND. ( ( time_since_reference_point <= skip_time_do_lsm & 1743 .AND. simulated_time > 0.0_wp ) .OR. .NOT. land_surface ) .AND. & 1744 .NOT. urban_surface .AND. .NOT. downward ) THEN 1745 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1746 DO m = 1, surf%ns 1747 i = surf%i(m) 1748 j = surf%j(m) 1749 k = surf%k(m) 1750 surf%shf(m) = -surf%ts(m) * surf%us(m) * rho_air_zw(k-1) 1751 ENDDO 1752 ENDIF 1753 ! 1754 !-- Compute the vertical water flux 1755 IF ( .NOT. constant_waterflux .AND. humidity .AND. & 1756 ( ( time_since_reference_point <= skip_time_do_lsm .AND. simulated_time > 0.0_wp ) & 1757 .OR. .NOT. land_surface ) .AND. .NOT. urban_surface .AND. .NOT. downward ) & 1758 THEN 1759 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1760 DO m = 1, surf%ns 1761 i = surf%i(m) 1762 j = surf%j(m) 1763 k = surf%k(m) 1764 surf%qsws(m) = -surf%qs(m) * surf%us(m) * rho_air_zw(k-1) 1765 ENDDO 1766 ENDIF 1767 ! 1768 !-- Compute the vertical scalar flux 1769 IF ( .NOT. constant_scalarflux .AND. passive_scalar .AND. .NOT. downward ) THEN 1770 !$OMP PARALLEL DO PRIVATE( i, j ) 1771 DO m = 1, surf%ns 1772 i = surf%i(m) 1773 j = surf%j(m) 1774 surf%ssws(m) = -surf%ss(m) * surf%us(m) 1775 ENDDO 1776 ENDIF 1777 ! 1778 !-- Compute the vertical chemical species' flux 1779 DO lsp = 1, nvar 1780 IF ( .NOT. constant_csflux(lsp) .AND. air_chemistry .AND. .NOT. downward ) THEN 1781 !$OMP PARALLEL DO PRIVATE( i, j ) 1782 DO m = 1, surf%ns 1783 i = surf%i(m) 1784 j = surf%j(m) 1785 surf%cssws(lsp,m) = -surf%css(lsp,m) * surf%us(m) 1786 ENDDO 1787 ENDIF 1788 ENDDO 1789 1790 ! 1791 !-- Compute (turbulent) fluxes of cloud water content and cloud drop conc. 1792 IF ( bulk_cloud_model .AND. microphysics_morrison .AND. .NOT. downward) THEN 1793 !$OMP PARALLEL DO PRIVATE( i, j ) 1794 DO m = 1, surf%ns 1795 i = surf%i(m) 1796 j = surf%j(m) 1797 1798 surf%qcsws(m) = -surf%qcs(m) * surf%us(m) 1799 surf%ncsws(m) = -surf%ncs(m) * surf%us(m) 1800 ENDDO 1801 ENDIF 1802 ! 1803 !-- Compute (turbulent) fluxes of rain water content and rain drop conc. 1804 IF ( bulk_cloud_model .AND. microphysics_seifert .AND. .NOT. downward) THEN 1805 !$OMP PARALLEL DO PRIVATE( i, j ) 1806 DO m = 1, surf%ns 1807 i = surf%i(m) 1808 j = surf%j(m) 1809 1810 surf%qrsws(m) = -surf%qrs(m) * surf%us(m) 1811 surf%nrsws(m) = -surf%nrs(m) * surf%us(m) 1812 ENDDO 1813 ENDIF 1814 1815 ! 1816 !-- Bottom boundary condition for the TKE. 1817 IF ( ibc_e_b == 2 ) THEN 1818 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1819 DO m = 1, surf%ns 1820 i = surf%i(m) 1821 j = surf%j(m) 1822 k = surf%k(m) 1823 1824 e(k,j,i) = ( surf%us(m) / 0.1_wp )**2 1825 ! 1826 !-- As a test: cm = 0.4 1827 ! e(k,j,i) = ( us(j,i) / 0.4_wp )**2 1828 e(k-1,j,i) = e(k,j,i) 1829 1830 ENDDO 1831 ENDIF 1832 ! 1833 !-- Calcuate surface fluxes at vertical surfaces. No stability is considered. 1834 ELSE 1835 ! 1836 !-- Compute usvs l={0,1} and vsus l={2,3} 1837 IF ( mom_uv ) THEN 1838 ! 1839 !-- Generalize computation by introducing flags. At north- and south-facing surfaces 1840 !-- u-component is used, at east- and west-facing surfaces v-component is used. 1841 flag_u = MERGE( 1.0_wp, 0.0_wp, l == 0 .OR. l == 1 ) 1842 flag_v = MERGE( 1.0_wp, 0.0_wp, l == 2 .OR. l == 3 ) 1843 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1844 DO m = 1, surf%ns 1845 i = surf%i(m) 1846 j = surf%j(m) 1847 k = surf%k(m) 1848 1849 z_mo = surf%z_mo(m) 1850 1851 surf%mom_flux_uv(m) = kappa * ( flag_u * u(k,j,i) + flag_v * v(k,j,i) ) / & 1852 LOG( z_mo / surf%z0(m) ) 1853 1854 surf%mom_flux_uv(m) = - surf%mom_flux_uv(m) * surf%us(m) 1855 ENDDO 1856 ENDIF 1857 ! 1858 !-- Compute wsus l={0,1} and wsvs l={2,3} 1859 IF ( mom_w ) THEN 1860 !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo ) 1861 DO m = 1, surf%ns 1862 i = surf%i(m) 1863 j = surf%j(m) 1864 k = surf%k(m) 1865 1866 z_mo = surf%z_mo(m) 1867 1868 surf%mom_flux_w(m) = kappa * w(k,j,i) / LOG( z_mo / surf%z0(m) ) 1869 1870 surf%mom_flux_w(m) = - surf%mom_flux_w(m) * surf%us(m) 1871 ENDDO 1872 ENDIF 1873 ! 1874 !-- Compute momentum fluxes used for subgrid-scale TKE production at vertical surfaces. In 1875 !-- constrast to the calculated momentum fluxes at vertical surfaces before, which are defined on 1876 !-- the u/v/w-grid, respectively), the TKE fluxes are defined at the scalar grid. 1877 !-- 1878 IF ( mom_tke ) THEN 1879 ! 1880 !-- Precalculate velocity components at scalar grid point. 1881 ALLOCATE( u_i(1:surf%ns) ) 1882 ALLOCATE( v_i(1:surf%ns) ) 1883 ALLOCATE( w_i(1:surf%ns) ) 1884 1885 IF ( l == 0 .OR. l == 1 ) THEN 1886 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1887 DO m = 1, surf%ns 1888 i = surf%i(m) 1889 j = surf%j(m) 1890 k = surf%k(m) 1891 1892 u_i(m) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) 1893 v_i(m) = 0.0_wp 1894 w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) 1895 ENDDO 1896 ELSE 1897 !$OMP PARALLEL DO PRIVATE( i, j, k ) 1898 DO m = 1, surf%ns 1899 i = surf%i(m) 1900 j = surf%j(m) 1901 k = surf%k(m) 1902 1903 u_i(m) = 0.0_wp 1904 v_i(m) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) 1905 w_i(m) = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) ) 1906 ENDDO 1907 ENDIF 1908 1909 !$OMP PARALLEL DO PRIVATE( i, j, dum, z_mo ) 1910 DO m = 1, surf%ns 1911 i = surf%i(m) 1912 j = surf%j(m) 1913 1914 z_mo = surf%z_mo(m) 1915 1916 dum = kappa / LOG( z_mo / surf%z0(m) ) 1917 ! 1918 !-- usvs (l=0,1) and vsus (l=2,3) 1919 surf%mom_flux_tke(0,m) = dum * ( u_i(m) + v_i(m) ) 1920 ! 1921 !-- wsvs (l=0,1) and wsus (l=2,3) 1922 surf%mom_flux_tke(1,m) = dum * w_i(m) 1923 1924 surf%mom_flux_tke(0:1,m) = - surf%mom_flux_tke(0:1,m) * surf%us(m) 1925 ENDDO 1926 ! 1927 !-- Deallocate temporary arrays 1928 DEALLOCATE( u_i ) 1929 DEALLOCATE( v_i ) 1930 DEALLOCATE( w_i ) 1931 ENDIF 1932 ENDIF 1933 1934 END SUBROUTINE calc_surface_fluxes 1935 1936 1937 !--------------------------------------------------------------------------------------------------! 1938 ! Description: 1939 ! ------------ 1940 !> Calculates temperature near surface (10 cm) for indoor model or 2 m temperature for output. 1941 !--------------------------------------------------------------------------------------------------! 1942 SUBROUTINE calc_pt_near_surface ( z_char ) 1943 1944 IMPLICIT NONE 1945 1946 CHARACTER(LEN = *), INTENT(IN) :: z_char !< string identifier to identify z level 1947 1948 INTEGER(iwp) :: i !< grid index x-dimension 1949 INTEGER(iwp) :: j !< grid index y-dimension 1950 INTEGER(iwp) :: k !< grid index z-dimension 1951 INTEGER(iwp) :: m !< running index for surface elements 1952 1953 1954 SELECT CASE ( z_char) 1955 1956 CASE ( '10cm' ) 1957 1958 DO m = 1, surf%ns 1959 1960 i = surf%i(m) 1961 j = surf%j(m) 1962 k = surf%k(m) 1963 1964 surf%pt_10cm(m) = surf%pt_surface(m) + surf%ts(m) / kappa & 1965 * ( LOG( 0.1_wp / surf%z0h(m) ) - psi_h( 0.1_wp / surf%ol(m) ) & 1966 + psi_h( surf%z0h(m) / surf%ol(m) ) ) 1967 ENDDO 1968 1969 END SELECT 1970 1971 END SUBROUTINE calc_pt_near_surface 1972 1973 1974 !--------------------------------------------------------------------------------------------------! 1975 ! Description: 1976 ! ------------ 1977 !> Integrated stability function for momentum. 1978 !--------------------------------------------------------------------------------------------------! 1979 FUNCTION psi_m( zeta ) 1980 !$ACC ROUTINE SEQ 1981 1982 USE kinds 1983 1984 IMPLICIT NONE 1985 1986 REAL(wp) :: psi_m !< Integrated similarity function result 1987 REAL(wp) :: zeta !< Stability parameter z/L 1988 REAL(wp) :: x !< dummy variable 1989 1990 REAL(wp), PARAMETER :: a = 1.0_wp !< constant 1991 REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant 1992 REAL(wp), PARAMETER :: c = 5.0_wp !< constant 1993 REAL(wp), PARAMETER :: d = 0.35_wp !< constant 1994 REAL(wp), PARAMETER :: c_d_d = c / d !< constant 1995 REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant 1996 1997 1998 IF ( zeta < 0.0_wp ) THEN 1999 x = SQRT( SQRT( 1.0_wp - 16.0_wp * zeta ) ) 2000 psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2 & 2001 * ( 1.0_wp + x**2 ) * 0.125_wp ) 2002 ELSE 2003 2004 psi_m = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - a * zeta - bc_d_d 2005 ! 2006 !-- Old version for stable conditions (only valid for z/L < 0.5) psi_m = - 5.0_wp * zeta 2007 2008 ENDIF 2009 2010 END FUNCTION psi_m 2011 2012 2013 !--------------------------------------------------------------------------------------------------! 2014 ! Description: 2015 !------------ 2016 !> Integrated stability function for heat and moisture. 2017 !--------------------------------------------------------------------------------------------------! 2018 FUNCTION psi_h( zeta ) 2019 !$ACC ROUTINE SEQ 2020 2021 USE kinds 2022 2023 IMPLICIT NONE 2024 2025 REAL(wp) :: psi_h !< Integrated similarity function result 2026 REAL(wp) :: zeta !< Stability parameter z/L 2027 REAL(wp) :: x !< dummy variable 2028 2029 REAL(wp), PARAMETER :: a = 1.0_wp !< constant 2030 REAL(wp), PARAMETER :: b = 0.66666666666_wp !< constant 2031 REAL(wp), PARAMETER :: c = 5.0_wp !< constant 2032 REAL(wp), PARAMETER :: d = 0.35_wp !< constant 2033 REAL(wp), PARAMETER :: c_d_d = c / d !< constant 2034 REAL(wp), PARAMETER :: bc_d_d = b * c / d !< constant 2035 2036 2037 IF ( zeta < 0.0_wp ) THEN 2038 x = SQRT( 1.0_wp - 16.0_wp * zeta ) 2039 psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp ) 2040 ELSE 2041 psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp & 2042 + 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d + 1.0_wp 2043 ! 2044 !-- Old version for stable conditions (only valid for z/L < 0.5) 2045 !-- psi_h = - 5.0_wp * zeta 2046 ENDIF 2047 2048 END FUNCTION psi_h 2049 2050 2051 !--------------------------------------------------------------------------------------------------! 2144 2052 ! Description: 2145 2053 ! ------------ … … 2147 2055 !> 2148 2056 !> @author Hauke Wurps 2149 !------------------------------------------------------------------------------ !2150 FUNCTION phi_m( zeta )2151 2152 2153 IMPLICIT NONE2154 2155 REAL(wp) :: phi_m!< Value of the function2156 REAL(wp) :: zeta!< Stability parameter z/L2157 2158 REAL(wp), PARAMETER :: a = 16.0_wp!< constant2159 REAL(wp), PARAMETER :: c = 5.0_wp!< constant2160 2161 2162 2163 2164 2165 2166 2167 2057 !--------------------------------------------------------------------------------------------------! 2058 FUNCTION phi_m( zeta ) 2059 !$ACC ROUTINE SEQ 2060 2061 IMPLICIT NONE 2062 2063 REAL(wp) :: phi_m !< Value of the function 2064 REAL(wp) :: zeta !< Stability parameter z/L 2065 2066 REAL(wp), PARAMETER :: a = 16.0_wp !< constant 2067 REAL(wp), PARAMETER :: c = 5.0_wp !< constant 2068 2069 IF ( zeta < 0.0_wp ) THEN 2070 phi_m = 1.0_wp / SQRT( SQRT( 1.0_wp - a * zeta ) ) 2071 ELSE 2072 phi_m = 1.0_wp + c * zeta 2073 ENDIF 2074 2075 END FUNCTION phi_m 2168 2076 2169 2077 END MODULE surface_layer_fluxes_mod -
palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
r4559 r4562 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Parts of r4559 re-formatted 28 ! 29 ! 4559 2020-06-11 08:51:48Z raasch 27 30 ! File re-formatted to follow the PALM coding standard 28 31 ! … … 1375 1378 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1376 1379 !-- Additional factors are added to improve the variance of v and w 1377 dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) * MERGE( 1.0_wp, 0.0_wp,&1378 BTEST( wall_flags_total_0(k,j,i), 1 ) )1380 dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) * & 1381 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 1379 1382 ENDDO 1380 1383 ENDDO … … 1393 1396 dist_yz(k,j,3) = MIN( ( SQRT( a33(k) / MAXVAL( a33 ) ) * 1.3_wp ) * & 1394 1397 ( a31(k) * fu_yz(k,j) + a32(k) * fv_yz(k,j) + a33(k) & 1395 * fw_yz(k,j) ), 3.0_wp ) * MERGE( 1.0_wp, 0.0_wp,&1396 1398 * fw_yz(k,j) ), 3.0_wp ) * & 1399 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 1397 1400 ENDDO 1398 1401 ENDDO … … 1412 1415 IF ( myidx == id_stg_right ) i = nxr+1 1413 1416 1414 mc_factor_l(1) = SUM( dist_yz(nzb:nzt,nys:nyn,1) * MERGE( 1.0_wp, 0.0_wp,&1415 1417 mc_factor_l(1) = SUM( dist_yz(nzb:nzt,nys:nyn,1) * & 1418 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 1 ) ) ) 1416 1419 1417 1420 IF ( myidx == id_stg_left ) i = nxl-1 1418 1421 IF ( myidx == id_stg_right ) i = nxr+1 1419 1422 1420 mc_factor_l(2) = SUM( dist_yz(nzb:nzt,nysv:nyn,2) * MERGE( 1.0_wp, 0.0_wp,&1421 1422 mc_factor_l(3) = SUM( dist_yz(nzb:nzt,nys:nyn,3) * MERGE( 1.0_wp, 0.0_wp,&1423 1423 mc_factor_l(2) = SUM( dist_yz(nzb:nzt,nysv:nyn,2) * & 1424 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,nysv:nyn,i), 2 ) ) ) 1425 mc_factor_l(3) = SUM( dist_yz(nzb:nzt,nys:nyn,3) * & 1426 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 3 ) ) ) 1424 1427 1425 1428 #if defined( __parallel ) … … 1435 1438 IF ( myidx == id_stg_right ) i = nxr+1 1436 1439 1437 dist_yz(:,nys:nyn,1) = ( dist_yz(:,nys:nyn,1) - mc_factor(1) ) * MERGE( 1.0_wp, 0.0_wp,&1438 BTEST( wall_flags_total_0(:,nys:nyn,i), 1 ) )1440 dist_yz(:,nys:nyn,1) = ( dist_yz(:,nys:nyn,1) - mc_factor(1) ) * & 1441 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,nys:nyn,i), 1 ) ) 1439 1442 1440 1443 … … 1442 1445 IF ( myidx == id_stg_right ) i = nxr+1 1443 1446 1444 dist_yz(:,nys:nyn,2) = ( dist_yz(:,nys:nyn,2) - mc_factor(2) ) * MERGE( 1.0_wp, 0.0_wp,&1445 BTEST( wall_flags_total_0(:,nys:nyn,i), 2 ) )1446 1447 dist_yz(:,nys:nyn,3) = ( dist_yz(:,nys:nyn,3) - mc_factor(3) ) * MERGE( 1.0_wp, 0.0_wp,&1448 BTEST( wall_flags_total_0(:,nys:nyn,i), 3 ) )1447 dist_yz(:,nys:nyn,2) = ( dist_yz(:,nys:nyn,2) - mc_factor(2) ) * & 1448 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,nys:nyn,i), 2 ) ) 1449 1450 dist_yz(:,nys:nyn,3) = ( dist_yz(:,nys:nyn,3) - mc_factor(3) ) * & 1451 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,nys:nyn,i), 3 ) ) 1449 1452 ! 1450 1453 !-- Add disturbances … … 1459 1462 DO j = nys, nyn 1460 1463 DO k = nzb, nzt+1 1461 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) ) 1462 *MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,0), 1 ) )1463 v(k,j,-nbgp:-1) = ( mean_inflow_profiles(k,2) + dist_yz(k,j,2) )&1464 *MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,-1), 2 ) )1465 w(k,j,-nbgp:-1) = dist_yz(k,j,3) * MERGE( 1.0_wp, 0.0_wp,&1466 1464 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) ) * & 1465 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,0), 1 ) ) 1466 v(k,j,-nbgp:-1) = ( mean_inflow_profiles(k,2) + dist_yz(k,j,2) ) * & 1467 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,-1), 2 ) ) 1468 w(k,j,-nbgp:-1) = dist_yz(k,j,3) * & 1469 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,-1), 3 ) ) 1467 1470 ENDDO 1468 1471 ENDDO … … 1471 1474 DO j = nys, nyn 1472 1475 DO k = nzb+1, nzt 1473 u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) ) * MERGE( 1.0_wp, 0.0_wp,&1474 BTEST( wall_flags_total_0(k,j,0), 1 ) )1476 u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) ) * & 1477 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,0), 1 ) ) 1475 1478 u(k,j,-1) = u(k,j,0) 1476 v(k,j,-1) = ( v(k,j,-1) + dist_yz(k,j,2) ) * MERGE( 1.0_wp, 0.0_wp,&1477 BTEST( wall_flags_total_0(k,j,-1), 2 ) )1478 w(k,j,-1) = ( w(k,j,-1) + dist_yz(k,j,3) ) * MERGE( 1.0_wp, 0.0_wp,&1479 BTEST( wall_flags_total_0(k,j,-1), 3 ) )1479 v(k,j,-1) = ( v(k,j,-1) + dist_yz(k,j,2) ) * & 1480 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,-1), 2 ) ) 1481 w(k,j,-1) = ( w(k,j,-1) + dist_yz(k,j,3) ) * & 1482 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,-1), 3 ) ) 1480 1483 ENDDO 1481 1484 ENDDO … … 1485 1488 DO j = nys, nyn 1486 1489 DO k = nzb+1, nzt 1487 u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) ) * MERGE( 1.0_wp, 0.0_wp,&1488 BTEST( wall_flags_total_0(k,j,nxr+1), 1 ) )1489 v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) ) * MERGE( 1.0_wp, 0.0_wp,&1490 1491 w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) ) * MERGE( 1.0_wp, 0.0_wp,&1492 BTEST( wall_flags_total_0(k,j,nxr+1), 3 ) )1490 u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) ) * & 1491 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 1 ) ) 1492 v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) ) * & 1493 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 2 ) ) 1494 w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) ) * & 1495 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr+1), 3 ) ) 1493 1496 ENDDO 1494 1497 ENDDO … … 1556 1559 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1557 1560 !-- Additional factors are added to improve the variance of v and w 1558 dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) * MERGE( 1.0_wp, 0.0_wp,&1559 BTEST( wall_flags_total_0(k,j,i), 1 ) )1561 dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) * & 1562 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 1560 1563 1561 1564 dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL( a33 ) ) * 1.3_wp ) & 1562 1565 * ( a31(k) * fu_xz(k,i) + a32(k) * fv_xz(k,i) + a33(k) & 1563 * fw_xz(k,i) ), 3.0_wp ) * MERGE( 1.0_wp, 0.0_wp,&1564 BTEST( wall_flags_total_0(k,j,i), 3 ) )1566 * fw_xz(k,i) ), 3.0_wp ) * & 1567 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 1565 1568 ENDDO 1566 1569 ENDDO … … 1579 1582 IF ( myidy == id_stg_north ) j = nyn+1 1580 1583 1581 mc_factor_l(2) = SUM( dist_xz(nzb:nzt,nxl:nxr,2) * MERGE( 1.0_wp, 0.0_wp,&1582 1584 mc_factor_l(2) = SUM( dist_xz(nzb:nzt,nxl:nxr,2) * & 1585 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 2 ) ) ) 1583 1586 1584 1587 IF ( myidy == id_stg_south ) j = nys-1 1585 1588 IF ( myidy == id_stg_north ) j = nyn+1 1586 1589 1587 mc_factor_l(1) = SUM( dist_xz(nzb:nzt,nxlu:nxr,1) * MERGE( 1.0_wp, 0.0_wp,&1588 1589 mc_factor_l(3) = SUM( dist_xz(nzb:nzt,nxl:nxr,3) * MERGE( 1.0_wp, 0.0_wp,&1590 1590 mc_factor_l(1) = SUM( dist_xz(nzb:nzt,nxlu:nxr,1) * & 1591 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,j,nxlu:nxr), 1 ) ) ) 1592 mc_factor_l(3) = SUM( dist_xz(nzb:nzt,nxl:nxr,3) * & 1593 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 3 ) ) ) 1591 1594 1592 1595 #if defined( __parallel ) … … 1601 1604 IF ( myidy == id_stg_north ) j = nyn+1 1602 1605 1603 dist_xz(:,nxl:nxr,2) = ( dist_xz(:,nxl:nxr,2) - mc_factor(2) ) * MERGE( 1.0_wp, 0.0_wp,&1604 BTEST( wall_flags_total_0(:,j,nxl:nxr), 2 ) )1606 dist_xz(:,nxl:nxr,2) = ( dist_xz(:,nxl:nxr,2) - mc_factor(2) ) * & 1607 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,j,nxl:nxr), 2 ) ) 1605 1608 1606 1609 … … 1608 1611 IF ( myidy == id_stg_north ) j = nyn+1 1609 1612 1610 dist_xz(:,nxl:nxr,1) = ( dist_xz(:,nxl:nxr,1) - mc_factor(1) ) * MERGE( 1.0_wp, 0.0_wp,&1611 BTEST( wall_flags_total_0(:,j,nxl:nxr), 1 ) )1612 1613 dist_xz(:,nxl:nxr,3) = ( dist_xz(:,nxl:nxr,3) - mc_factor(3) ) * MERGE( 1.0_wp, 0.0_wp,&1614 BTEST( wall_flags_total_0(:,j,nxl:nxr), 3 ) )1613 dist_xz(:,nxl:nxr,1) = ( dist_xz(:,nxl:nxr,1) - mc_factor(1) ) * & 1614 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,j,nxl:nxr), 1 ) ) 1615 1616 dist_xz(:,nxl:nxr,3) = ( dist_xz(:,nxl:nxr,3) - mc_factor(3) ) * & 1617 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(:,j,nxl:nxr), 3 ) ) 1615 1618 ! 1616 1619 !-- Add disturbances … … 1632 1635 DO i = nxl, nxr 1633 1636 DO k = nzb+1, nzt 1634 u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) ) * MERGE( 1.0_wp, 0.0_wp,&1635 BTEST( wall_flags_total_0(k,nyn+1,i), 1 ) )1636 v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) ) * MERGE( 1.0_wp, 0.0_wp,&1637 BTEST( wall_flags_total_0(k,nyn+1,i), 2 ) )1638 w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) ) * MERGE( 1.0_wp, 0.0_wp,&1639 BTEST( wall_flags_total_0(k,nyn+1,i), 3 ) )1637 u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) ) * & 1638 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 1 ) ) 1639 v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) ) * & 1640 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 2 ) ) 1641 w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) ) * & 1642 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn+1,i), 3 ) ) 1640 1643 ENDDO 1641 1644 ENDDO … … 2022 2025 !-- u'u' and v'v'. Assume isotropy. Note, add a small negative number to the denominator, else 2023 2026 !-- the mergpe-function can crash if scale_l is zero. 2024 r11(k) = scale_us**2 * ( MERGE( 0.35_wp & 2025 * ABS( - zi_ribulk / ( kappa * scale_l - 10E-4_wp ) )**( 2.0_wp / 3.0_wp ), & 2026 0.0_wp, scale_l < 0.0_wp ) + 5.0_wp - 4.0_wp * zzi ) * blend 2027 r11(k) = scale_us**2 * ( MERGE( 0.35_wp * & 2028 ABS( - zi_ribulk / ( kappa * scale_l - 10E-4_wp ) )**( 2.0_wp / 3.0_wp ), & 2029 0.0_wp, scale_l < 0.0_wp ) + 5.0_wp - 4.0_wp * zzi & 2030 ) * blend 2027 2031 2028 2032 r22(k) = r11(k)
Note: See TracChangeset
for help on using the changeset viewer.
