Changeset 4559 for palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
- Timestamp:
- Jun 11, 2020 8:51:48 AM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
r4535 r4559 1 1 !> @synthetic_turbulence_generator_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 Software7 ! Foundation, either version 3 of the License, or (at your option) any later8 ! version.9 ! 10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR12 ! 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 with15 ! PALM. If not, see <http://www.gnu.org/licenses/>.16 ! 17 ! Copyright 2017-2020 Leibniz Universitaet Hannover18 ! ------------------------------------------------------------------------------!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/>. 15 ! 16 ! Copyright 1997-2020 Leibniz Universitaet Hannover 17 !--------------------------------------------------------------------------------------------------! 18 ! 19 19 ! 20 20 ! Current revisions: … … 25 25 ! ----------------- 26 26 ! $Id$ 27 ! bugfix for restart data format query 28 ! 27 ! File re-formatted to follow the PALM coding standard 28 ! 29 ! 4535 2020-05-15 12:07:23Z raasch 30 ! Bugfix for restart data format query 31 ! 29 32 ! 4495 2020-04-13 20:11:20Z raasch 30 ! restart data handling with MPI-IO added31 ! 33 ! Restart data handling with MPI-IO added 34 ! 32 35 ! 4481 2020-03-31 18:55:54Z maronga 33 ! bugfix: cpp-directives for serial mode added, dummy statements to prevent compile errors added34 ! 36 ! Bugfix: cpp-directives for serial mode added, dummy statements to prevent compile errors added 37 ! 35 38 ! 4442 2020-03-04 19:21:13Z suehring 36 39 ! Set back turbulent length scale to 8 x grid spacing in the parametrized mode 37 40 ! (was accidantly changed). 38 ! 41 ! 39 42 ! 4441 2020-03-04 19:20:35Z suehring 40 43 ! Correct misplaced preprocessor directive 41 ! 44 ! 42 45 ! 4438 2020-03-03 20:49:28Z suehring 43 46 ! Performance optimizations in velocity-seed calculation: 44 ! - random number array is only defined and computed locally (except for 45 ! normalization to zero mean and unit variance) 46 ! - parallel random number generator is applied independent on the 2D random 47 ! numbers in other routines 48 ! - option to decide wheter velocity seeds are computed locally without any 49 ! further communication or are computed by all processes along the 50 ! communicator 51 ! 47 ! - Random number array is only defined and computed locally (except for normalization to zero mean 48 ! and unit variance) 49 ! - Parallel random number generator is applied independent on the 2D random numbers in other 50 ! routines 51 ! - Option to decide wheter velocity seeds are computed locally without any further communication 52 ! or are computed by all processes along the communicator 53 ! 52 54 ! 4346 2019-12-18 11:55:56Z motisi 53 ! Introduction of wall_flags_total_0, which currently sets bits based on static 54 ! topographyinformation used in wall_flags_static_055 ! 55 ! Introduction of wall_flags_total_0, which currently sets bits based on static topography 56 ! information used in wall_flags_static_0 57 ! 56 58 ! 4335 2019-12-12 16:39:05Z suehring 57 59 ! Commentation of last commit 58 ! 60 ! 59 61 ! 4332 2019-12-10 19:44:12Z suehring 60 ! Limit initial velocity seeds in restart runs, if not the seed calculation 61 ! may become unstable. Further, minor bugfix in initial velocity seed 62 ! calculation. 63 ! 62 ! Limit initial velocity seeds in restart runs, if not the seed calculation may become unstable. 63 ! Further, minor bugfix in initial velocity seed calculation. 64 ! 64 65 ! 4329 2019-12-10 15:46:36Z motisi 65 66 ! Renamed wall_flags_0 to wall_flags_static_0 66 ! 67 ! 67 68 ! 4309 2019-11-26 18:49:59Z suehring 68 ! Computation of velocity seeds optimized. This implies that random numbers 69 ! are computed now using the parallel random number generator. Random numbers 70 ! are now only computed and normalized locally, while distributed over all 71 ! mpi ranks afterwards, instead of computing random numbers on a global array. 72 ! Further, the number of calls for the time-consuming velocity-seed generation 73 ! is reduced - now the left and right, as well as the north and south boundary 74 ! share the same velocity-seed matrices. 75 ! 69 ! Computation of velocity seeds optimized. This implies that random numbers are computed now using 70 ! the parallel random number generator. Random numbers are now only computed and normalized locally, 71 ! while distributed over all mpi ranks afterwards, instead of computing random numbers on a global 72 ! array. 73 ! Further, the number of calls for the time-consuming velocity-seed generation is reduced - now the 74 ! left and right, as well as the north and south boundary share the same velocity-seed matrices. 75 ! 76 76 ! 4182 2019-08-22 15:20:23Z scharf 77 77 ! Corrected "Former revisions" section 78 ! 78 ! 79 79 ! 4148 2019-08-08 11:26:00Z suehring 80 80 ! Remove unused variable 81 ! 81 ! 82 82 ! 4144 2019-08-06 09:11:47Z raasch 83 ! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.84 ! 83 ! Relational operators .EQ., .NE., etc. replaced by ==, /=, etc. 84 ! 85 85 ! 4071 2019-07-03 20:02:00Z suehring 86 ! Bugfix, initialize mean_inflow_profiles in case turbulence and inflow 87 ! information is not read from file.88 ! 86 ! Bugfix, initialize mean_inflow_profiles in case turbulence and inflow information is not read from 87 ! file. 88 ! 89 89 ! 4022 2019-06-12 11:52:39Z suehring 90 90 ! Several bugfixes and improvements 91 ! - revise bias correction of the imposed perturbations (correction via volume 92 ! flow can create instabilities in case the mean volume flow is close to zero) 93 ! - introduce lower limits in calculation of coefficient matrix, else the 94 ! calculation may become numerically unstable 95 ! - impose perturbations every timestep, even though no new set of perturbations 96 ! is generated in case dt_stg_call /= dt_3d 97 ! - Implement a gradual decrease of Reynolds stress and length scales above 98 ! ABL height (within 1 length scale above ABL depth to 1/10) rather than a 99 ! discontinuous decrease 91 ! - Revise bias correction of the imposed perturbations (correction via volume flow can create 92 ! instabilities in case the mean volume flow is close to zero) 93 ! - Introduce lower limits in calculation of coefficient matrix, else the calculation may become 94 ! numerically unstable 95 ! - Impose perturbations every timestep, even though no new set of perturbations is generated in 96 ! case dt_stg_call /= dt_3d 97 ! - Implement a gradual decrease of Reynolds stress and length scales above ABL height (within 1 98 ! length scale above ABL depth to 1/10) rather than a discontinuous decrease 100 99 ! - Bugfix in non-nested case: use ABL height for parametrized turbulence 101 ! 100 ! 102 101 ! 3987 2019-05-22 09:52:13Z kanani 103 102 ! Introduce alternative switch for debug output during timestepping 104 ! 103 ! 105 104 ! 3938 2019-04-29 16:06:25Z suehring 106 105 ! Remove unused variables 107 ! 106 ! 108 107 ! 3937 2019-04-29 15:09:07Z suehring 109 ! Minor bugfix in case of a very early restart where mc_factor is sill not 110 ! present. 111 ! Some modification and fixing of potential bugs in the calculation of scaling 112 ! parameters used for synthetic turbulence parametrization. 113 ! 108 ! Minor bugfix in case of a very early restart where mc_factor is sill not present. 109 ! Some modification and fixing of potential bugs in the calculation of scaling parameters used for 110 ! synthetic turbulence parametrization. 111 ! 114 112 ! 3909 2019-04-17 09:13:25Z suehring 115 113 ! Minor bugfix for last commit 116 ! 114 ! 117 115 ! 3900 2019-04-16 15:17:43Z suehring 118 116 ! Missing re-calculation of perturbation seeds in case of restarts 119 ! 117 ! 120 118 ! 3891 2019-04-12 17:52:01Z suehring 121 ! Bugfix in initialization in case of restart runs. 122 ! 119 ! Bugfix in initialization in case of restart runs. 120 ! 123 121 ! 3885 2019-04-11 11:29:34Z kanani 124 ! Changes related to global restructuring of location messages and introduction 125 ! of additional debugmessages126 ! 127 ! 128 ! removed unused variables129 ! 122 ! Changes related to global restructuring of location messages and introduction of additional debug 123 ! messages 124 ! 125 ! 126 ! Removed unused variables 127 ! 130 128 ! 3719 2019-02-06 13:10:18Z kanani 131 ! Removed log_point measurement from stg_init, since this part is counted to 132 ! log_point(2) 'initialisation' already. Moved other log_points to calls of133 ! the subroutines in time_integrationfor better overview.134 ! 129 ! Removed log_point measurement from stg_init, since this part is counted to log_point(2) 130 ! 'initialization' already. Moved other log_points to calls of the subroutines in time_integration 131 ! for better overview. 132 ! 135 133 ! 2259 2017-06-08 09:09:11Z gronemeier 136 134 ! Initial revision … … 143 141 ! Description: 144 142 ! ------------ 145 !> The module generates turbulence at the inflow boundary based on a method by 146 !> Xie and Castro (2008) utilizing a Lund rotation (Lund, 1998) and a mass-flux 147 !> correction by Kim et al. (2013). 148 !> The turbulence is correlated based on length scales in y- and z-direction and 149 !> a time scale for each velocity component. The profiles of length and time 150 !> scales, mean u, v, w, e and pt, and all components of the Reynolds stress 151 !> tensor can be either read from file STG_PROFILES, or will be parametrized 152 !> within the boundary layer. 143 !> The module generates turbulence at the inflow boundary based on a method by Xie and Castro (2008) 144 !> utilizing a Lund rotation (Lund, 1998) and a mass-flux correction by Kim et al. (2013). 145 !> The turbulence is correlated based on length scales in y- and z-direction and a time scale for 146 !> each velocity component. The profiles of length and time scales, mean u, v, w, e and pt, and all 147 !> components of the Reynolds stress tensor can be either read from file STG_PROFILES, or will be 148 !> parametrized within the boundary layer. 153 149 !> 154 150 !> @todo test restart 155 !> enable cyclic_fill156 !> implement turbulence generation for e and pt151 !> Enable cyclic_fill 152 !> Implement turbulence generation for e and pt 157 153 !> @todo Input of height-constant length scales via namelist 158 154 !> @note <Enter notes on the module> 159 !> @bug Height information from input file is not used. Profiles from input 160 !> must match with current PALM grid. 161 !> In case of restart, velocity seeds differ from precursor run if a11, 162 !> a22, or a33 are zero. 163 !------------------------------------------------------------------------------! 155 !> @bug Height information from input file is not used. Profiles from input must match with current 156 !> PALM grid. 157 !> In case of restart, velocity seeds differ from precursor run if a11, a22, or a33 are zero. 158 !--------------------------------------------------------------------------------------------------! 164 159 MODULE synthetic_turbulence_generator_mod 165 160 166 161 167 USE arrays_3d, &168 ONLY: dzw, &169 ddzw, &170 drho_air, &171 mean_inflow_profiles, &172 q,&173 q_init,&174 pt,&175 pt_init,&176 u, &177 u_init, &178 v, &179 v_init, &180 w, &181 zu, &162 USE arrays_3d, & 163 ONLY: dzw, & 164 ddzw, & 165 drho_air, & 166 mean_inflow_profiles, & 167 pt, & 168 pt_init, & 169 q, & 170 q_init, & 171 u, & 172 u_init, & 173 v, & 174 v_init, & 175 w, & 176 zu, & 182 177 zw 183 178 184 USE basic_constants_and_equations_mod, &185 ONLY: g, &186 kappa, &179 USE basic_constants_and_equations_mod, & 180 ONLY: g, & 181 kappa, & 187 182 pi 188 183 189 USE control_parameters, &190 ONLY: bc_lr, &191 bc_ns, &192 child_domain, &193 coupling_char, &194 debug_output_timestep, &195 dt_3d, &196 e_init, &197 humidity, &198 initializing_actions, &199 intermediate_timestep_count, &200 intermediate_timestep_count_max, &201 length, &202 message_string, &203 nesting_offline, &204 neutral, &205 num_mean_inflow_profiles, &206 random_generator, &207 rans_mode, &208 restart_data_format_output, &209 restart_string, &210 syn_turb_gen, &211 time_since_reference_point, &184 USE control_parameters, & 185 ONLY: bc_lr, & 186 bc_ns, & 187 child_domain, & 188 coupling_char, & 189 debug_output_timestep, & 190 dt_3d, & 191 e_init, & 192 humidity, & 193 initializing_actions, & 194 intermediate_timestep_count, & 195 intermediate_timestep_count_max, & 196 length, & 197 message_string, & 198 nesting_offline, & 199 neutral, & 200 num_mean_inflow_profiles, & 201 random_generator, & 202 rans_mode, & 203 restart_data_format_output, & 204 restart_string, & 205 syn_turb_gen, & 206 time_since_reference_point, & 212 207 turbulent_inflow 213 208 214 USE cpulog, &215 ONLY: cpu_log, &209 USE cpulog, & 210 ONLY: cpu_log, & 216 211 log_point_s 217 212 218 USE grid_variables, &219 ONLY: ddx, &220 ddy, &221 dx, &213 USE grid_variables, & 214 ONLY: ddx, & 215 ddy, & 216 dx, & 222 217 dy 223 218 224 USE indices, &225 ONLY: nbgp, &226 nz, &227 nzb, &228 nzt, &229 nx, &230 nxl, &231 nxlu, &232 nxr, &233 ny, &234 nys, &235 nysv, &236 nyn, &219 USE indices, & 220 ONLY: nbgp, & 221 nz, & 222 nzb, & 223 nzt, & 224 nx, & 225 nxl, & 226 nxlu, & 227 nxr, & 228 ny, & 229 nys, & 230 nysv, & 231 nyn, & 237 232 wall_flags_total_0 238 233 … … 243 238 #endif 244 239 245 USE nesting_offl_mod, &246 ONLY: nesting_offl_calc_zi, &240 USE nesting_offl_mod, & 241 ONLY: nesting_offl_calc_zi, & 247 242 zi_ribulk 248 243 249 USE pegrid, &250 ONLY: comm1dx, &251 comm1dy, &252 comm2d, &253 ierr, &254 myidx, &255 myidy, &244 USE pegrid, & 245 ONLY: comm1dx, & 246 comm1dy, & 247 comm2d, & 248 ierr, & 249 myidx, & 250 myidy, & 256 251 pdims 257 252 258 USE pmc_interface, &253 USE pmc_interface, & 259 254 ONLY : rans_mode_parent 260 255 261 USE random_generator_parallel, &262 ONLY: init_parallel_random_generator, &263 random_dummy, &264 random_number_parallel, &256 USE random_generator_parallel, & 257 ONLY: init_parallel_random_generator, & 258 random_dummy, & 259 random_number_parallel, & 265 260 random_seed_parallel 266 261 … … 269 264 wrd_mpi_io 270 265 271 USE transpose_indices, &272 ONLY: nzb_x,&273 nzt_x274 275 USE surface_mod, &276 ONLY: surf_def_h, &277 surf_lsm_h, &266 USE transpose_indices, & 267 ONLY: nzb_x, & 268 nzt_x 269 270 USE surface_mod, & 271 ONLY: surf_def_h, & 272 surf_lsm_h, & 278 273 surf_usm_h 279 274 … … 284 279 #endif 285 280 286 287 LOGICAL :: velocity_seed_initialized = .FALSE. !< true after first call of stg_main288 LOGICAL :: parametrize_inflow_turbulence = .FALSE. !< flag indicating that inflow turbulence is either read from file (.FALSE.) or if it parametrized289 LOGICAL :: use_syn_turb_gen = .FALSE. !< switch to use synthetic turbulence generator290 LOGICAL :: compute_velocity_seeds_local = .TRUE. !< switch to decide whether velocity seeds are computed locally or if computation291 !< is distributed over several processes292 281 293 282 INTEGER(iwp) :: id_stg_left !< left lateral boundary core id in case of turbulence generator … … 307 296 #endif 308 297 309 INTEGER(iwp), DIMENSION(3) :: nr_non_topo_xz = 0 !< number of non-topography grid points at xz cross-sections,310 !< required for bias correction of imposed perturbations311 INTEGER(iwp), DIMENSION(3) :: nr_non_topo_yz = 0 !< number of non-topography grid points at yz cross-sections,312 !< required for bias correction of imposed perturbations313 298 INTEGER(iwp), DIMENSION(3) :: nr_non_topo_xz = 0 !< number of non-topography grid points at xz cross-sections, 299 !< required for bias correction of imposed perturbations 300 INTEGER(iwp), DIMENSION(3) :: nr_non_topo_yz = 0 !< number of non-topography grid points at yz cross-sections, 301 !< required for bias correction of imposed perturbations 302 314 303 #if defined( __parallel ) 315 304 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_xz !< displacement for MPI_GATHERV … … 328 317 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz !< length scale of w in z direction (in gp) 329 318 330 INTEGER(isp), DIMENSION(:), ALLOCATABLE :: id_rand_xz !< initial random IDs at xz inflow boundary 331 INTEGER(isp), DIMENSION(:), ALLOCATABLE :: id_rand_yz !< initial random IDs at yz inflow boundary 332 INTEGER(isp), DIMENSION(:,:), ALLOCATABLE :: seq_rand_xz !< initial random seeds at xz inflow boundary 333 INTEGER(isp), DIMENSION(:,:), ALLOCATABLE :: seq_rand_yz !< initial random seeds at yz inflow boundary 334 335 REAL(wp) :: blend !< value to create gradually and smooth blending of Reynolds stress and length 336 !< scales above the boundary layer 337 REAL(wp) :: blend_coeff = -2.3_wp !< coefficient used to ensure that blending functions decreases to 1/10 after 338 !< one length scale above ABL top 339 REAL(wp) :: d_l !< blend_coeff/length_scale 340 REAL(wp) :: length_scale !< length scale, default is 8 x minimum grid spacing 341 REAL(wp) :: dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics 342 REAL(wp) :: dt_stg_call = 0.0_wp !< time interval for calling synthetic turbulence generator 343 REAL(wp) :: scale_l !< scaling parameter used for turbulence parametrization - Obukhov length 344 REAL(wp) :: scale_us !< scaling parameter used for turbulence parametrization - friction velocity 345 REAL(wp) :: scale_wm !< scaling parameter used for turbulence parametrization - momentum scale 346 REAL(wp) :: time_stg_adjust = 0.0_wp !< time counter for adjusting turbulence information 347 REAL(wp) :: time_stg_call = 0.0_wp !< time counter for calling generator 348 349 REAL(wp), DIMENSION(3) :: mc_factor = 1.0_wp !< correction factor for the u,v,w-components to maintain original mass flux 350 351 352 REAL(wp),DIMENSION(:), ALLOCATABLE :: r11 !< Reynolds parameter 353 REAL(wp),DIMENSION(:), ALLOCATABLE :: r21 !< Reynolds parameter 354 REAL(wp),DIMENSION(:), ALLOCATABLE :: r22 !< Reynolds parameter 355 REAL(wp),DIMENSION(:), ALLOCATABLE :: r31 !< Reynolds parameter 356 REAL(wp),DIMENSION(:), ALLOCATABLE :: r32 !< Reynolds parameter 357 REAL(wp),DIMENSION(:), ALLOCATABLE :: r33 !< Reynolds parameter 358 359 REAL(wp), DIMENSION(:), ALLOCATABLE :: a11 !< coefficient for Lund rotation 360 REAL(wp), DIMENSION(:), ALLOCATABLE :: a21 !< coefficient for Lund rotation 361 REAL(wp), DIMENSION(:), ALLOCATABLE :: a22 !< coefficient for Lund rotation 362 REAL(wp), DIMENSION(:), ALLOCATABLE :: a31 !< coefficient for Lund rotation 363 REAL(wp), DIMENSION(:), ALLOCATABLE :: a32 !< coefficient for Lund rotation 364 REAL(wp), DIMENSION(:), ALLOCATABLE :: a33 !< coefficient for Lund rotation 365 REAL(wp), DIMENSION(:), ALLOCATABLE :: tu !< Lagrangian time scale of u 366 REAL(wp), DIMENSION(:), ALLOCATABLE :: tv !< Lagrangian time scale of v 367 REAL(wp), DIMENSION(:), ALLOCATABLE :: tw !< Lagrangian time scale of w 368 369 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bux !< filter function for u in x direction 370 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy !< filter function for u in y direction 371 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz !< filter function for u in z direction 372 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvx !< filter function for v in x direction 373 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy !< filter function for v in y direction 374 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz !< filter function for v in z direction 375 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwx !< filter function for w in y direction 376 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy !< filter function for w in y direction 377 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz !< filter function for w in z direction 378 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_xz !< velocity seed for u at xz plane 379 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_xz !< velocity seed for u at xz plane with new random number 380 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_yz !< velocity seed for u at yz plane 381 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_yz !< velocity seed for u at yz plane with new random number 382 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_xz !< velocity seed for v at xz plane 383 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_xz !< velocity seed for v at xz plane with new random number 384 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_yz !< velocity seed for v at yz plane 385 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_yz !< velocity seed for v at yz plane with new random number 386 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_xz !< velocity seed for w at xz plane 387 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_xz !< velocity seed for w at xz plane with new random number 388 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_yz !< velocity seed for w at yz plane 389 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz !< velocity seed for w at yz plane with new random number 390 391 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dist_xz !< imposed disturbances at north/south boundary 392 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dist_yz !< imposed disturbances at north/south boundary 319 INTEGER(isp), DIMENSION(:), ALLOCATABLE :: id_rand_xz !< initial random IDs at xz inflow boundary 320 INTEGER(isp), DIMENSION(:), ALLOCATABLE :: id_rand_yz !< initial random IDs at yz inflow boundary 321 INTEGER(isp), DIMENSION(:,:), ALLOCATABLE :: seq_rand_xz !< initial random seeds at xz inflow boundary 322 INTEGER(isp), DIMENSION(:,:), ALLOCATABLE :: seq_rand_yz !< initial random seeds at yz inflow boundary 323 324 325 LOGICAL :: compute_velocity_seeds_local = .TRUE. !< switch to decide whether velocity seeds are computed locally 326 !< or if computation is distributed over several processes 327 LOGICAL :: parametrize_inflow_turbulence = .FALSE. !< flag indicating that inflow turbulence is either read from file 328 !< (.FALSE.) or if it parametrized 329 LOGICAL :: use_syn_turb_gen = .FALSE. !< switch to use synthetic turbulence generator 330 LOGICAL :: velocity_seed_initialized = .FALSE. !< true after first call of stg_main 331 332 333 REAL(wp) :: blend !< value to create gradually and smooth blending of Reynolds stress and length 334 !< scales above the boundary layer 335 REAL(wp) :: blend_coeff = -2.3_wp !< coefficient used to ensure that blending functions decreases to 1/10 after 336 !< one length scale above ABL top 337 REAL(wp) :: d_l !< blend_coeff/length_scale 338 REAL(wp) :: dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics 339 REAL(wp) :: dt_stg_call = 0.0_wp !< time interval for calling synthetic turbulence generator 340 REAL(wp) :: length_scale !< length scale, default is 8 x minimum grid spacing 341 REAL(wp) :: scale_l !< scaling parameter used for turbulence parametrization - Obukhov length 342 REAL(wp) :: scale_us !< scaling parameter used for turbulence parametrization - friction velocity 343 REAL(wp) :: scale_wm !< scaling parameter used for turbulence parametrization - momentum scale 344 REAL(wp) :: time_stg_adjust = 0.0_wp !< time counter for adjusting turbulence information 345 REAL(wp) :: time_stg_call = 0.0_wp !< time counter for calling generator 346 347 REAL(wp), DIMENSION(3) :: mc_factor = 1.0_wp !< correction factor for the u,v,w-components to maintain original mass flux 348 349 350 REAL(wp),DIMENSION(:), ALLOCATABLE :: r11 !< Reynolds parameter 351 REAL(wp),DIMENSION(:), ALLOCATABLE :: r21 !< Reynolds parameter 352 REAL(wp),DIMENSION(:), ALLOCATABLE :: r22 !< Reynolds parameter 353 REAL(wp),DIMENSION(:), ALLOCATABLE :: r31 !< Reynolds parameter 354 REAL(wp),DIMENSION(:), ALLOCATABLE :: r32 !< Reynolds parameter 355 REAL(wp),DIMENSION(:), ALLOCATABLE :: r33 !< Reynolds parameter 356 357 REAL(wp), DIMENSION(:), ALLOCATABLE :: a11 !< coefficient for Lund rotation 358 REAL(wp), DIMENSION(:), ALLOCATABLE :: a21 !< coefficient for Lund rotation 359 REAL(wp), DIMENSION(:), ALLOCATABLE :: a22 !< coefficient for Lund rotation 360 REAL(wp), DIMENSION(:), ALLOCATABLE :: a31 !< coefficient for Lund rotation 361 REAL(wp), DIMENSION(:), ALLOCATABLE :: a32 !< coefficient for Lund rotation 362 REAL(wp), DIMENSION(:), ALLOCATABLE :: a33 !< coefficient for Lund rotation 363 REAL(wp), DIMENSION(:), ALLOCATABLE :: tu !< Lagrangian time scale of u 364 REAL(wp), DIMENSION(:), ALLOCATABLE :: tv !< Lagrangian time scale of v 365 REAL(wp), DIMENSION(:), ALLOCATABLE :: tw !< Lagrangian time scale of w 366 367 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bux !< filter function for u in x direction 368 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy !< filter function for u in y direction 369 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz !< filter function for u in z direction 370 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvx !< filter function for v in x direction 371 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy !< filter function for v in y direction 372 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz !< filter function for v in z direction 373 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwx !< filter function for w in y direction 374 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy !< filter function for w in y direction 375 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz !< filter function for w in z direction 376 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_xz !< velocity seed for u at xz plane 377 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_xz !< velocity seed for u at xz plane with new random number 378 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_yz !< velocity seed for u at yz plane 379 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_yz !< velocity seed for u at yz plane with new random number 380 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_xz !< velocity seed for v at xz plane 381 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_xz !< velocity seed for v at xz plane with new random number 382 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_yz !< velocity seed for v at yz plane 383 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_yz !< velocity seed for v at yz plane with new random number 384 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_xz !< velocity seed for w at xz plane 385 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_xz !< velocity seed for w at xz plane with new random number 386 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_yz !< velocity seed for w at yz plane 387 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz !< velocity seed for w at yz plane with new random number 388 389 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dist_xz !< imposed disturbances at north/south boundary 390 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dist_yz !< imposed disturbances at north/south boundary 393 391 394 392 ! … … 464 462 ! 465 463 !-- Public interfaces 466 PUBLIC stg_adjust, stg_check_parameters, stg_header, stg_init, stg_main, & 467 stg_parin, stg_rrd_global, stg_wrd_global 464 PUBLIC stg_adjust, & 465 stg_check_parameters, & 466 stg_header, & 467 stg_init, & 468 stg_main, & 469 stg_parin, & 470 stg_rrd_global, & 471 stg_wrd_global 468 472 469 473 ! 470 474 !-- Public variables 471 PUBLIC dt_stg_call, dt_stg_adjust, id_stg_left, id_stg_north, & 472 id_stg_right, id_stg_south, parametrize_inflow_turbulence, & 473 time_stg_adjust, time_stg_call, use_syn_turb_gen 475 PUBLIC dt_stg_call, & 476 dt_stg_adjust, & 477 id_stg_left, & 478 id_stg_north, & 479 id_stg_right, & 480 id_stg_south, & 481 parametrize_inflow_turbulence, & 482 time_stg_adjust, & 483 time_stg_call, & 484 use_syn_turb_gen 474 485 475 486 … … 477 488 478 489 479 !------------------------------------------------------------------------------ !490 !--------------------------------------------------------------------------------------------------! 480 491 ! Description: 481 492 ! ------------ 482 493 !> Check parameters routine for synthetic turbulence generator 483 !------------------------------------------------------------------------------ !494 !--------------------------------------------------------------------------------------------------! 484 495 SUBROUTINE stg_check_parameters 485 496 486 IF ( .NOT. use_syn_turb_gen .AND. .NOT. rans_mode .AND. &497 IF ( .NOT. use_syn_turb_gen .AND. .NOT. rans_mode .AND. & 487 498 nesting_offline ) THEN 488 message_string = 'Synthetic turbulence generator is required ' // & 489 'if offline nesting is applied and PALM operates ' // & 490 'in LES mode.' 499 message_string = 'Synthetic turbulence generator is required ' // & 500 'if offline nesting is applied and PALM operates in LES mode.' 491 501 CALL message( 'stg_check_parameters', 'PA0520', 0, 0, 0, 6, 0 ) 492 502 ENDIF 493 503 494 IF ( .NOT. use_syn_turb_gen .AND. child_domain &504 IF ( .NOT. use_syn_turb_gen .AND. child_domain & 495 505 .AND. rans_mode_parent .AND. .NOT. rans_mode ) THEN 496 message_string = 'Synthetic turbulence generator is required ' // & 497 'when nesting is applied and parent operates in ' // & 498 'RANS-mode but current child in LES mode.' 506 message_string = 'Synthetic turbulence generator is required when nesting is applied ' // & 507 'and parent operates in RANS-mode but current child in LES mode.' 499 508 CALL message( 'stg_check_parameters', 'PA0524', 1, 2, 0, 6, 0 ) 500 509 ENDIF … … 502 511 IF ( use_syn_turb_gen ) THEN 503 512 504 IF ( child_domain .AND. .NOT. rans_mode .AND. & 505 .NOT. rans_mode_parent ) THEN 506 message_string = 'Using synthetic turbulence generator ' // & 507 'is not allowed in LES-LES nesting.' 513 IF ( child_domain .AND. .NOT. rans_mode .AND. .NOT. rans_mode_parent ) THEN 514 message_string = 'Using synthetic turbulence generator is not allowed in LES-LES nesting.' 508 515 CALL message( 'stg_check_parameters', 'PA0620', 1, 2, 0, 6, 0 ) 509 516 510 517 ENDIF 511 512 IF ( child_domain .AND. rans_mode .AND. & 513 rans_mode_parent ) THEN 514 message_string = 'Using synthetic turbulence generator ' // & 515 'is not allowed in RANS-RANS nesting.' 518 519 IF ( child_domain .AND. rans_mode .AND. rans_mode_parent ) THEN 520 message_string = 'Using synthetic turbulence generator is not allowed in RANS-RANS nesting.' 516 521 CALL message( 'stg_check_parameters', 'PA0621', 1, 2, 0, 6, 0 ) 517 522 518 523 ENDIF 519 520 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 521 522 IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0 & 523 .AND. INDEX( initializing_actions, 'read_restart_data' ) == 0 ) THEN 524 message_string = 'Using synthetic turbulence generator ' // & 525 'requires %initializing_actions = ' // & 526 '"set_constant_profiles" or "read_restart_data"' //& 527 ', if not offline nesting is applied.' 524 525 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 526 527 IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0 & 528 .AND. INDEX( initializing_actions, 'read_restart_data' ) == 0 ) THEN 529 message_string = 'Using synthetic turbulence generator requires ' // & 530 '%initializing_actions = "set_constant_profiles" or ' // & 531 ' "read_restart_data", if not offline nesting is applied.' 528 532 CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 ) 529 533 ENDIF 530 534 531 535 IF ( bc_lr /= 'dirichlet/radiation' ) THEN 532 message_string = 'Using synthetic turbulence generator ' // & 533 'requires &bc_lr = "dirichlet/radiation", ' // & 534 'if not offline nesting is applied.' 536 message_string = 'Using synthetic turbulence generator requires &bc_lr = ' // & 537 ' "dirichlet/radiation", if not offline nesting is applied.' 535 538 CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 ) 536 539 ENDIF 537 540 IF ( bc_ns /= 'cyclic' ) THEN 538 message_string = 'Using synthetic turbulence generator ' // & 539 'requires &bc_ns = "cyclic", ' // & 540 'if not offline nesting is applied.' 541 message_string = 'Using synthetic turbulence generator requires &bc_ns = ' // & 542 ' "cyclic", if not offline nesting is applied.' 541 543 CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 ) 542 544 ENDIF … … 545 547 546 548 IF ( turbulent_inflow ) THEN 547 message_string = 'Using synthetic turbulence generator ' // & 548 'in combination &with turbulent_inflow = .T. '// & 549 'is not allowed' 549 message_string = 'Using synthetic turbulence generator in combination ' // & 550 '&with turbulent_inflow = .T. is not allowed' 550 551 CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 ) 551 552 ENDIF … … 553 554 !-- Synthetic turbulence generator requires the parallel random generator 554 555 IF ( random_generator /= 'random-parallel' ) THEN 555 message_string = 'Using synthetic turbulence generator ' //&556 'r equires random_generator = random-parallel.'556 message_string = 'Using synthetic turbulence generator requires random_generator = ' // & 557 'random-parallel.' 557 558 CALL message( 'stg_check_parameters', 'PA0421', 1, 2, 0, 6, 0 ) 558 559 ENDIF … … 563 564 564 565 565 !------------------------------------------------------------------------------ !566 !--------------------------------------------------------------------------------------------------! 566 567 ! Description: 567 568 ! ------------ 568 569 !> Header output for synthetic turbulence generator 569 !------------------------------------------------------------------------------ !570 !--------------------------------------------------------------------------------------------------! 570 571 SUBROUTINE stg_header ( io ) 571 572 572 INTEGER(iwp), INTENT(IN) :: io 573 INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file 573 574 574 575 ! … … 580 581 WRITE( io, 3 ) 581 582 ENDIF 582 583 583 584 IF ( parametrize_inflow_turbulence ) THEN 584 585 WRITE( io, 4 ) dt_stg_adjust … … 587 588 ENDIF 588 589 589 1 FORMAT (//' Synthetic turbulence generator information:'/ &590 1 FORMAT (//' Synthetic turbulence generator information:'/ & 590 591 ' ------------------------------------------'/) 591 592 2 FORMAT (' synthetic turbulence generator is switched on') … … 597 598 598 599 599 !------------------------------------------------------------------------------ !600 !--------------------------------------------------------------------------------------------------! 600 601 ! Description: 601 602 ! ------------ 602 603 !> Initialization of the synthetic turbulence generator 603 !------------------------------------------------------------------------------ !604 !--------------------------------------------------------------------------------------------------! 604 605 SUBROUTINE stg_init 605 606 606 LOGICAL :: file_stg_exist = .FALSE. !< flag indicating whether parameter file for Reynolds stress and length scales exist607 608 607 #if defined( __parallel ) 609 INTEGER(KIND=MPI_ADDRESS_KIND) :: extent !< extent of new MPI type610 INTEGER(KIND=MPI_ADDRESS_KIND) :: tob =0 !< dummy variable608 INTEGER(KIND=MPI_ADDRESS_KIND) :: extent !< extent of new MPI type 609 INTEGER(KIND=MPI_ADDRESS_KIND) :: tob = 0 !< dummy variable 611 610 #endif 612 611 613 INTEGER(iwp) :: i 614 INTEGER(iwp) :: j 615 INTEGER(iwp) :: k 612 INTEGER(iwp) :: i !> grid index in x-direction 613 INTEGER(iwp) :: j !> loop index 614 INTEGER(iwp) :: k !< index 616 615 #if defined( __parallel ) 617 INTEGER(iwp) :: newtype 618 INTEGER(iwp) :: realsize 616 INTEGER(iwp) :: newtype !< dummy MPI type 617 INTEGER(iwp) :: realsize !< size of REAL variables 619 618 #endif 620 619 621 620 INTEGER(iwp), DIMENSION(3) :: nr_non_topo_xz_l = 0 !< number of non-topography grid points at xz-cross-section on subdomain 622 621 INTEGER(iwp), DIMENSION(3) :: nr_non_topo_yz_l = 0 !< number of non-topography grid points at yz-cross-section on subdomain 622 623 624 LOGICAL :: file_stg_exist = .FALSE. !< flag indicating whether parameter file for Reynolds stress and length scales exist 625 623 626 ! 624 627 !-- Dummy variables used for reading profiles 625 REAL(wp) :: d1 626 REAL(wp) :: d2 627 REAL(wp) :: d3 628 REAL(wp) :: d5 629 REAL(wp) :: luy 630 REAL(wp) :: luz 631 REAL(wp) :: lvy 632 REAL(wp) :: lvz 633 REAL(wp) :: lwy 634 REAL(wp) :: lwz 628 REAL(wp) :: d1 !< u profile 629 REAL(wp) :: d2 !< v profile 630 REAL(wp) :: d3 !< w profile 631 REAL(wp) :: d5 !< e profile 632 REAL(wp) :: luy !< length scale for u in y direction 633 REAL(wp) :: luz !< length scale for u in z direction 634 REAL(wp) :: lvy !< length scale for v in y direction 635 REAL(wp) :: lvz !< length scale for v in z direction 636 REAL(wp) :: lwy !< length scale for w in y direction 637 REAL(wp) :: lwz !< length scale for w in z direction 635 638 #if defined( __parallel ) 636 REAL(wp) :: nnz 639 REAL(wp) :: nnz !< increment used to determine processor decomposition of z-axis along x and y direction 637 640 #endif 638 REAL(wp) :: zz 641 REAL(wp) :: zz !< height 639 642 640 643 … … 643 646 #endif 644 647 ! 645 !-- Create mpi-datatypes for exchange in case of non-local but distributed 646 !-- computation of the velocity seeds. This option is useful in 647 !-- case large turbulent length scales are present, where the computational 648 !-- effort becomes large and need to be parallelized. For parameterized 649 !-- turbulence the length scales are small and computing the velocity seeds 650 !-- locally is faster (no overhead by communication). 648 !-- Create mpi-datatypes for exchange in case of non-local but distributed computation of the 649 !-- velocity seeds. This option is useful in case large turbulent length scales are present, where 650 !-- the computational effort becomes large and needs to be parallelized. For parameterized turbulence 651 !-- the length scales are small and computing the velocity seeds locally is faster (no overhead by 652 !-- communication). 651 653 IF ( .NOT. compute_velocity_seeds_local ) THEN 652 654 #if defined( __parallel ) 653 ! 655 ! 654 656 !-- Determine processor decomposition of z-axis along x- and y-direction 655 657 nnz = nz / pdims(1) 656 658 nzb_x_stg = 1 + myidx * INT( nnz ) 657 659 nzt_x_stg = ( myidx + 1 ) * INT( nnz ) 658 659 IF ( MOD( nz , pdims(1) ) /= 0 .AND. myidx == id_stg_right ) &660 661 IF ( MOD( nz , pdims(1) ) /= 0 .AND. myidx == id_stg_right ) & 660 662 nzt_x_stg = nzt_x_stg + myidx * ( nnz - INT( nnz ) ) 661 662 IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent &663 664 IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent & 663 665 .AND. .NOT. rans_mode ) ) THEN 664 666 nnz = nz / pdims(2) 665 667 nzb_y_stg = 1 + myidy * INT( nnz ) 666 668 nzt_y_stg = ( myidy + 1 ) * INT( nnz ) 667 668 IF ( MOD( nz , pdims(2) ) /= 0 .AND. myidy == id_stg_north ) &669 670 IF ( MOD( nz , pdims(2) ) /= 0 .AND. myidy == id_stg_north ) & 669 671 nzt_y_stg = nzt_y_stg + myidy * ( nnz - INT( nnz ) ) 670 672 ENDIF 671 672 ! 673 !-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted 674 !-- velocity seeds 673 674 ! 675 !-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted velocity seeds 675 676 CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr ) 676 677 extent = 1 * realsize 677 ! 678 !-- Set-up MPI datatyp to involve all cores for turbulence generation at yz 679 !-- layer 678 ! 679 !-- Set-up MPI datatyp to involve all cores for turbulence generation at yz layer 680 680 !-- stg_type_yz: yz-slice with vertical bounds nzb:nzt+1 681 CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyn-nys+1], &681 CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyn-nys+1], & 682 682 [1,nyn-nys+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr ) 683 683 CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz, ierr ) 684 684 CALL MPI_TYPE_COMMIT( stg_type_yz, ierr ) 685 685 CALL MPI_TYPE_FREE( newtype, ierr ) 686 686 687 687 ! stg_type_yz_small: yz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1 688 CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyn-nys+1], &688 CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyn-nys+1], & 689 689 [1,nyn-nys+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr ) 690 690 CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr ) 691 691 CALL MPI_TYPE_COMMIT( stg_type_yz_small, ierr ) 692 692 CALL MPI_TYPE_FREE( newtype, ierr ) 693 694 ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz693 694 ! Receive count and displacement for MPI_GATHERV in stg_generate_seed_yz 695 695 ALLOCATE( recv_count_yz(pdims(1)), displs_yz(pdims(1)) ) 696 696 697 697 recv_count_yz = nzt_x_stg-nzb_x_stg + 1 698 698 recv_count_yz(pdims(1)) = recv_count_yz(pdims(1)) + 1 699 699 700 700 DO j = 1, pdims(1) 701 701 displs_yz(j) = 0 + (nzt_x_stg-nzb_x_stg+1) * (j-1) 702 702 ENDDO 703 ! 704 !-- Set-up MPI datatyp to involve all cores for turbulence generation at xz 705 !-- layer 703 ! 704 !-- Set-up MPI datatyp to involve all cores for turbulence generation at xz layer 706 705 !-- stg_type_xz: xz-slice with vertical bounds nzb:nzt+1 707 706 IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent & … … 712 711 CALL MPI_TYPE_COMMIT( stg_type_xz, ierr ) 713 712 CALL MPI_TYPE_FREE( newtype, ierr ) 714 713 715 714 ! stg_type_yz_small: xz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1 716 715 CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_y_stg-nzb_y_stg+2,nxr-nxl+1], & … … 719 718 CALL MPI_TYPE_COMMIT( stg_type_xz_small, ierr ) 720 719 CALL MPI_TYPE_FREE( newtype, ierr ) 721 722 ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz720 721 ! Receive count and displacement for MPI_GATHERV in stg_generate_seed_yz 723 722 ALLOCATE( recv_count_xz(pdims(2)), displs_xz(pdims(2)) ) 724 723 725 724 recv_count_xz = nzt_y_stg-nzb_y_stg + 1 726 725 recv_count_xz(pdims(2)) = recv_count_xz(pdims(2)) + 1 727 726 728 727 DO j = 1, pdims(2) 729 728 displs_xz(j) = 0 + (nzt_y_stg-nzb_y_stg+1) * (j-1) 730 729 ENDDO 731 730 732 731 ENDIF 733 732 #endif … … 735 734 ! 736 735 !-- Allocate required arrays. 737 !-- In case no offline nesting or self nesting is used, the arrary 738 !-- mean_inflow profiles is required. Check if it is already allocated, else 739 !-- allocate and initialize it appropriately. Note, in case turbulence and 740 !-- inflow information is read from file, mean_inflow_profiles is already 741 !-- allocated and initialized appropriately. 742 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 736 !-- In case no offline nesting or self nesting is used, the arrary mean_inflow profiles is required. 737 !-- Check if it is already allocated, else allocate and initialize it appropriately. Note, in case 738 !-- turbulence and inflow information is read from file, mean_inflow_profiles is already allocated 739 !-- and initialized appropriately. 740 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 743 741 IF ( .NOT. ALLOCATED( mean_inflow_profiles ) ) THEN 744 742 ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) ) … … 747 745 mean_inflow_profiles(:,2) = v_init 748 746 ! 749 !-- Even though potential temperature and humidity are not perturbed, 750 !-- they need to be initialized appropriately.751 IF ( .NOT. neutral ) &747 !-- Even though potential temperature and humidity are not perturbed, they need to be 748 !-- initialized appropriately. 749 IF ( .NOT. neutral ) & 752 750 mean_inflow_profiles(:,4) = pt_init 753 IF ( humidity ) &754 mean_inflow_profiles(:,6) = q_init 755 ENDIF 751 IF ( humidity ) & 752 mean_inflow_profiles(:,6) = q_init 753 ENDIF 756 754 ENDIF 757 755 758 ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1), &759 a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1), &760 nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1), &761 nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1), &762 nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1), &763 r11(nzb:nzt+1), r21(nzb:nzt+1), r22(nzb:nzt+1), &764 r31(nzb:nzt+1), r32(nzb:nzt+1), r33(nzb:nzt+1), &765 tu(nzb:nzt+1), tv(nzb:nzt+1), tw(nzb:nzt+1) 766 756 ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1), & 757 a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1), & 758 nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1), & 759 nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1), & 760 nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1), & 761 r11(nzb:nzt+1), r21(nzb:nzt+1), r22(nzb:nzt+1), & 762 r31(nzb:nzt+1), r32(nzb:nzt+1), r33(nzb:nzt+1), & 763 tu(nzb:nzt+1), tv(nzb:nzt+1), tw(nzb:nzt+1) ) 764 767 765 ALLOCATE ( dist_xz(nzb:nzt+1,nxl:nxr,3) ) 768 766 ALLOCATE ( dist_yz(nzb:nzt+1,nys:nyn,3) ) … … 772 770 !-- Read inflow profile 773 771 !-- Height levels of profiles in input profile is as follows: 774 !-- zu: luy, luz, tu, lvy, lvz, tv, r11, r21, r22, d1, d2, d5 775 !-- zw: lwy, lwz, tw, r31, r32, r33, d3 772 !-- zu: luy, luz, tu, lvy, lvz, tv, r11, r21, r22, d1, d2, d5 zw: lwy, lwz, tw, r31, r32, r33, d3 776 773 !-- WARNING: zz is not used at the moment 777 INQUIRE( FILE = 'STG_PROFILES' // TRIM( coupling_char ), & 778 EXIST = file_stg_exist ) 774 INQUIRE( FILE = 'STG_PROFILES' // TRIM( coupling_char ), EXIST = file_stg_exist ) 779 775 780 776 IF ( file_stg_exist ) THEN 781 777 782 OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD', & 783 FORM='FORMATTED') 778 OPEN( 90, FILE = 'STG_PROFILES' // TRIM( coupling_char ), STATUS = 'OLD', FORM = 'FORMATTED' ) 784 779 ! 785 780 !-- Skip header … … 787 782 788 783 DO k = nzb+1, nzt+1 789 READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k), & 790 r11(k), r21(k), r22(k), r31(k), r32(k), r33(k), & 791 d1, d2, d3, d5 792 793 ! 794 !-- Convert length scales from meter to number of grid points. 784 READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k), r11(k), r21(k), & 785 r22(k), r31(k), r32(k), r33(k), d1, d2, d3, d5 786 787 ! 788 !-- Convert length scales from meter to number of grid points. 795 789 nuy(k) = INT( luy * ddy ) 796 790 nuz(k) = INT( luz * ddzw(k) ) … … 815 809 ! 816 810 !-- Set lenght scales at surface grid point 817 nuy(nzb) = nuy(nzb+1) 811 nuy(nzb) = nuy(nzb+1) 818 812 nuz(nzb) = nuz(nzb+1) 819 813 nvy(nzb) = nvy(nzb+1) … … 821 815 nwy(nzb) = nwy(nzb+1) 822 816 nwz(nzb) = nwz(nzb+1) 823 817 824 818 CLOSE( 90 ) 825 819 ! 826 !-- Calculate coefficient matrix from Reynolds stress tensor 827 !-- (Lund rotation) 820 !-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation) 828 821 CALL calc_coeff_matrix 829 822 ! 830 !-- No information about turbulence and its length scales are available. 831 !-- Instead, parametrize turbulence which is imposed at the boundaries. 832 !-- Set flag which indicates that turbulence is parametrized, which is done 833 !-- later when energy-balance models are already initialized. This is 834 !-- because the STG needs information about surface properties such as 835 !-- roughness to build 'realistic' turbulence profiles. 823 !-- No information about turbulence and its length scales are available. Instead, parametrize 824 !-- turbulence which is imposed at the boundaries. Set flag which indicates that turbulence is 825 !-- parametrized, which is done later when energy-balance models are already initialized. This is 826 !-- because the STG needs information about surface properties such as roughness to build 827 !-- 'realistic' turbulence profiles. 836 828 ELSE 837 829 ! 838 !-- Define length scale for the imposed turbulence, which is defined as 839 !-- 8 times the minimum gridspacing830 !-- Define length scale for the imposed turbulence, which is defined as 8 times the minimum grid 831 !-- spacing 840 832 length_scale = 8.0_wp * MIN( dx, dy, MINVAL( dzw ) ) 841 833 ! 842 !-- Define constant to gradually decreas e length scales and Reynolds stress843 !-- above ABL top. Assure that no zero length scales are used.834 !-- Define constant to gradually decreasing length scales and Reynolds stress above ABL top. 835 !-- Assure that no zero length scales are used. 844 836 d_l = blend_coeff / MAX( length_scale, dx, dy, MINVAL( dzw ) ) 845 837 ! … … 847 839 parametrize_inflow_turbulence = .TRUE. 848 840 ! 849 !-- In case of dirichlet inflow boundary conditions only at one lateral 850 !-- boundary, i.e. in the case no offline or self nesting is applied but 851 !-- synthetic turbulence shall be parametrized nevertheless, the 852 !-- boundary-layer depth need to determined first. 853 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) & 841 !-- In case of dirichlet inflow boundary conditions only at one lateral boundary, i.e. in the 842 !-- case no offline or self nesting is applied but synthetic turbulence shall be parametrized 843 !-- nevertheless, the boundary-layer depth needs to be determined first. 844 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) & 854 845 CALL nesting_offl_calc_zi 855 846 ! 856 !-- Determine boundary-layer depth, which is used to initialize lenght 857 !-- scales 847 !-- Determine boundary-layer depth, which is used to initialize lenght scales 858 848 CALL calc_scaling_variables 859 849 ! 860 !-- Initialize lenght and time scales, which in turn are used 861 !-- to initialize the filter functions. 850 !-- Initialize lenght and time scales, which in turn are used to initialize the filter functions. 862 851 CALL calc_length_and_time_scale 863 852 ! 864 !-- Parametrize Reynolds-stress tensor, diagonal elements as well 865 !-- as r21 (v'u'), r31 (w'u'), r32 (w'v'). Parametrization follows 866 !-- Rotach et al. (1996) and is based on boundary-layer depth, 867 !-- friction velocity and velocity scale. 853 !-- Parametrize Reynolds-stress tensor, diagonal elements as well as r21 (v'u'), r31 (w'u'), 854 !-- r32 (w'v'). Parametrization follows Rotach et al. (1996) and is based on boundary-layer depth, 855 !-- friction velocity and velocity scale. 868 856 CALL parametrize_reynolds_stress 869 ! 870 !-- Calculate coefficient matrix from Reynolds stress tensor 871 !-- (Lund rotation) 857 ! 858 !-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation) 872 859 CALL calc_coeff_matrix 873 860 874 861 ENDIF 875 862 876 863 ! 877 !-- Assign initial profiles. Note, this is only required if turbulent 878 !-- inflow from the left is desired, not in case of any of the 879 !-- nesting (offline or self nesting) approaches. 880 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 864 !-- Assign initial profiles. Note, this is only required if turbulent inflow from the left is 865 !-- desired, not in case of any of the nesting (offline or self nesting) approaches. 866 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 881 867 u_init = mean_inflow_profiles(:,1) 882 868 v_init = mean_inflow_profiles(:,2) … … 884 870 e_init = MAXVAL( mean_inflow_profiles(:,5) ) 885 871 ENDIF 886 872 887 873 ! 888 874 !-- Define the size of the filter functions and allocate them. … … 891 877 ! arrays must be large enough to cover the largest length scale 892 878 DO k = nzb, nzt+1 893 j = MAX( ABS( nux(k)), ABS(nuy(k)), ABS(nuz(k)),&894 ABS( nvx(k)), ABS(nvy(k)), ABS(nvz(k)),&895 ABS( nwx(k)), ABS(nwy(k)), ABS(nwz(k)))879 j = MAX( ABS( nux(k) ), ABS( nuy(k) ), ABS( nuz(k) ), & 880 ABS( nvx(k) ), ABS( nvy(k) ), ABS( nvz(k) ), & 881 ABS( nwx(k) ), ABS( nwy(k) ), ABS( nwz(k) ) ) 896 882 IF ( j > mergp ) mergp = j 897 883 ENDDO … … 900 886 ! mergp = mergp 901 887 902 ALLOCATE ( bux(-mergp:mergp,nzb:nzt+1), &903 buy(-mergp:mergp,nzb:nzt+1), &904 buz(-mergp:mergp,nzb:nzt+1), &905 bvx(-mergp:mergp,nzb:nzt+1), &906 bvy(-mergp:mergp,nzb:nzt+1), &907 bvz(-mergp:mergp,nzb:nzt+1), &908 bwx(-mergp:mergp,nzb:nzt+1), &909 bwy(-mergp:mergp,nzb:nzt+1), &910 bwz(-mergp:mergp,nzb:nzt+1) 888 ALLOCATE ( bux(-mergp:mergp,nzb:nzt+1), & 889 buy(-mergp:mergp,nzb:nzt+1), & 890 buz(-mergp:mergp,nzb:nzt+1), & 891 bvx(-mergp:mergp,nzb:nzt+1), & 892 bvy(-mergp:mergp,nzb:nzt+1), & 893 bvz(-mergp:mergp,nzb:nzt+1), & 894 bwx(-mergp:mergp,nzb:nzt+1), & 895 bwy(-mergp:mergp,nzb:nzt+1), & 896 bwz(-mergp:mergp,nzb:nzt+1) ) 911 897 912 898 ! 913 899 !-- Allocate velocity seeds for turbulence at xz-layer 914 ALLOCATE ( fu_xz( nzb:nzt+1,nxl:nxr), fuo_xz(nzb:nzt+1,nxl:nxr),&915 fv_xz( nzb:nzt+1,nxl:nxr), fvo_xz(nzb:nzt+1,nxl:nxr),&916 fw_xz( nzb:nzt+1,nxl:nxr), fwo_xz(nzb:nzt+1,nxl:nxr))900 ALLOCATE ( fu_xz(nzb:nzt+1,nxl:nxr), fuo_xz(nzb:nzt+1,nxl:nxr), & 901 fv_xz(nzb:nzt+1,nxl:nxr), fvo_xz(nzb:nzt+1,nxl:nxr), & 902 fw_xz(nzb:nzt+1,nxl:nxr), fwo_xz(nzb:nzt+1,nxl:nxr) ) 917 903 918 904 ! 919 905 !-- Allocate velocity seeds for turbulence at yz-layer 920 ALLOCATE ( fu_yz( nzb:nzt+1,nys:nyn), fuo_yz(nzb:nzt+1,nys:nyn),&921 fv_yz( nzb:nzt+1,nys:nyn), fvo_yz(nzb:nzt+1,nys:nyn),&922 fw_yz( nzb:nzt+1,nys:nyn), fwo_yz(nzb:nzt+1,nys:nyn))906 ALLOCATE ( fu_yz(nzb:nzt+1,nys:nyn), fuo_yz(nzb:nzt+1,nys:nyn), & 907 fv_yz(nzb:nzt+1,nys:nyn), fvo_yz(nzb:nzt+1,nys:nyn), & 908 fw_yz(nzb:nzt+1,nys:nyn), fwo_yz(nzb:nzt+1,nys:nyn) ) 923 909 924 910 fu_xz = 0.0_wp … … 953 939 954 940 ! 955 !-- In case of restart, calculate velocity seeds fu, fv, fw from former 956 ! time step. 957 ! Bug: fu, fv, fw are different in those heights where a11, a22, a33 958 ! are 0 compared to the prerun. This is mostly for k=nzt+1. 941 !-- In case of restart, calculate velocity seeds fu, fv, fw from former time step. 942 ! Bug: fu, fv, fw are different in those heights where a11, a22, a33 are 0 compared to the prerun. 943 !-- This is mostly for k=nzt+1. 959 944 IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN 960 945 IF ( myidx == id_stg_left .OR. myidx == id_stg_right ) THEN … … 972 957 973 958 IF ( a22(k) > 10E-8_wp ) THEN 974 fv_yz(k,j) = ( v(k,j,i) - & 975 a21(k) * fu_yz(k,j) - v_init(k) ) / a22(k) 959 fv_yz(k,j) = ( v(k,j,i) - a21(k) * fu_yz(k,j) - v_init(k) ) / a22(k) 976 960 ELSE 977 961 fv_yz(k,j) = 10E-8_wp … … 979 963 980 964 IF ( a33(k) > 10E-8_wp ) THEN 981 fw_yz(k,j) = ( w(k,j,i) - & 982 a31(k) * fu_yz(k,j) - a32(k) * & 983 fv_yz(k,j) ) / a33(k) 965 fw_yz(k,j) = ( w(k,j,i) - a31(k) * fu_yz(k,j) - a32(k) * fv_yz(k,j) ) / a33(k) 984 966 ELSE 985 967 fw_yz(k,j) = 10E-8_wp … … 988 970 ENDDO 989 971 ENDIF 990 972 991 973 IF ( myidy == id_stg_south .OR. myidy == id_stg_north ) THEN 992 974 … … 997 979 DO k = nzb, nzt+1 998 980 ! 999 !-- In case the correlation coefficients are very small, the 1000 !-- velocity seeds may become very large finally creating 1001 !-- numerical instabilities in the synthetic turbulence generator. 1002 !-- Empirically, a value of 10E-8 seems to be sufficient. 981 !-- In case the correlation coefficients are very small, the velocity seeds may become 982 !-- very large finally creating numerical instabilities in the synthetic turbulence 983 !-- generator. Empirically, a value of 10E-8 seems to be sufficient. 1003 984 IF ( a11(k) > 10E-8_wp ) THEN 1004 985 fu_xz(k,i) = ( u(k,j,i) - u_init(k) ) / a11(k) … … 1008 989 1009 990 IF ( a22(k) > 10E-8_wp ) THEN 1010 fv_xz(k,i) = ( v(k,j,i) - & 1011 a21(k) * fu_xz(k,i) - v_init(k) ) / a22(k) 991 fv_xz(k,i) = ( v(k,j,i) - a21(k) * fu_xz(k,i) - v_init(k) ) / a22(k) 1012 992 ELSE 1013 993 fv_xz(k,i) = 10E-8_wp … … 1015 995 1016 996 IF ( a33(k) > 10E-8_wp ) THEN 1017 fw_xz(k,i) = ( w(k,j,i) - & 1018 a31(k) * fu_xz(k,i) - & 1019 a32(k) * fv_xz(k,i) ) / a33(k) 997 fw_xz(k,i) = ( w(k,j,i) - a31(k) * fu_xz(k,i) - a32(k) * fv_xz(k,i) ) / a33(k) 1020 998 ELSE 1021 999 fw_xz(k,i) = 10E-8_wp … … 1027 1005 ENDIF 1028 1006 ! 1029 !-- Count the number of non-topography grid points at the boundaries where 1030 !-- perturbations are imposed. This number is later used for bias corrections 1031 !-- of the perturbations, i.e. to force that their mean is zero. Please note, 1032 !-- due to the asymetry of u and v along x and y direction, respectively, 1033 !-- different cases must be distinguished. 1007 !-- Count the number of non-topography grid points at the boundaries where perturbations are imposed. 1008 !-- This number is later used for bias corrections of the perturbations, i.e. to force their mean to 1009 !-- be zero. Please note, due to the asymetry of u and v along x and y direction, respectively, 1010 !-- different cases must be distinguished. 1034 1011 IF ( myidx == id_stg_left .OR. myidx == id_stg_right ) THEN 1035 1012 ! … … 1037 1014 IF ( myidx == id_stg_left ) i = nxl 1038 1015 IF ( myidx == id_stg_right ) i = nxr+1 1039 1040 nr_non_topo_yz_l(1) = SUM( & 1041 MERGE( 1, 0, & 1042 BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 1 ) ) ) 1043 ! 1044 !-- Number of grid points where perturbations are imposed on v and w 1016 1017 nr_non_topo_yz_l(1) = SUM( MERGE( 1, 0, BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 1 ) ) ) 1018 ! 1019 !-- Number of grid points where perturbations are imposed on v and w 1045 1020 IF ( myidx == id_stg_left ) i = nxl-1 1046 1021 IF ( myidx == id_stg_right ) i = nxr+1 1047 1048 nr_non_topo_yz_l(2) = SUM( & 1049 MERGE( 1, 0, & 1050 BTEST( wall_flags_total_0(nzb:nzt,nysv:nyn,i), 2 ) ) ) 1051 nr_non_topo_yz_l(3) = SUM( & 1052 MERGE( 1, 0, & 1053 BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 3 ) ) ) 1054 1022 1023 nr_non_topo_yz_l(2) = SUM( MERGE( 1, 0, BTEST( wall_flags_total_0(nzb:nzt,nysv:nyn,i), 2 ) ) ) 1024 nr_non_topo_yz_l(3) = SUM( MERGE( 1, 0, BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 3 ) ) ) 1025 1055 1026 #if defined( __parallel ) 1056 CALL MPI_ALLREDUCE( nr_non_topo_yz_l, nr_non_topo_yz, 3, MPI_INTEGER, & 1057 MPI_SUM, comm1dy, ierr ) 1027 CALL MPI_ALLREDUCE( nr_non_topo_yz_l, nr_non_topo_yz, 3, MPI_INTEGER, MPI_SUM, comm1dy, ierr ) 1058 1028 #else 1059 1029 nr_non_topo_yz = nr_non_topo_yz_l 1060 #endif 1030 #endif 1061 1031 ENDIF 1062 1032 1063 1033 IF ( myidy == id_stg_south .OR. myidy == id_stg_north ) THEN 1064 1034 ! … … 1066 1036 IF ( myidy == id_stg_south ) j = nys 1067 1037 IF ( myidy == id_stg_north ) j = nyn+1 1068 1069 nr_non_topo_xz_l(2) = SUM( & 1070 MERGE( 1, 0, & 1071 BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 2 ) ) ) 1072 ! 1073 !-- Number of grid points where perturbations are imposed on u and w 1038 1039 nr_non_topo_xz_l(2) = SUM( MERGE( 1, 0, BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 2 ) ) ) 1040 ! 1041 !-- Number of grid points where perturbations are imposed on u and w 1074 1042 IF ( myidy == id_stg_south ) j = nys-1 1075 1043 IF ( myidy == id_stg_north ) j = nyn+1 1076 1077 nr_non_topo_xz_l(1) = SUM( & 1078 MERGE( 1, 0, & 1079 BTEST( wall_flags_total_0(nzb:nzt,j,nxlu:nxr), 1 ) ) ) 1080 nr_non_topo_xz_l(3) = SUM( & 1081 MERGE( 1, 0, & 1082 BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 3 ) ) ) 1083 1044 1045 nr_non_topo_xz_l(1) = SUM( MERGE( 1, 0, BTEST( wall_flags_total_0(nzb:nzt,j,nxlu:nxr), 1 ) ) ) 1046 nr_non_topo_xz_l(3) = SUM( MERGE( 1, 0, BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 3 ) ) ) 1047 1084 1048 #if defined( __parallel ) 1085 CALL MPI_ALLREDUCE( nr_non_topo_xz_l, nr_non_topo_xz, 3, MPI_INTEGER, & 1086 MPI_SUM, comm1dx, ierr ) 1049 CALL MPI_ALLREDUCE( nr_non_topo_xz_l, nr_non_topo_xz, 3, MPI_INTEGER, MPI_SUM, comm1dx, ierr ) 1087 1050 #else 1088 1051 nr_non_topo_xz = nr_non_topo_xz_l 1089 #endif 1052 #endif 1090 1053 ENDIF 1091 1054 ! 1092 !-- Initialize random number generator at xz- and yz-layers. Random numbers 1093 !-- are initialized at each core. In case there is only inflow from the left, 1094 !-- it is sufficient to generate only random numbers for the yz-layer, else 1095 !-- random numbers for the xz-layer are also required. 1055 !-- Initialize random number generator at xz- and yz-layers. Random numbers are initialized at each 1056 !-- core. In case there is only inflow from the left, it is sufficient to generate only random 1057 !-- numbers for the yz-layer, else random numbers for the xz-layer are also required. 1096 1058 ALLOCATE ( id_rand_yz(-mergp+nys:nyn+mergp) ) 1097 1059 ALLOCATE ( seq_rand_yz(5,-mergp+nys:nyn+mergp) ) … … 1099 1061 seq_rand_yz = 0 1100 1062 1101 CALL init_parallel_random_generator( ny, -mergp+nys, nyn+mergp, & 1102 id_rand_yz, seq_rand_yz ) 1103 1104 IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent & 1105 .AND. .NOT. rans_mode ) ) THEN 1063 CALL init_parallel_random_generator( ny, -mergp+nys, nyn+mergp, id_rand_yz, seq_rand_yz ) 1064 1065 IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent & 1066 .AND. .NOT. rans_mode ) ) THEN 1106 1067 ALLOCATE ( id_rand_xz(-mergp+nxl:nxr+mergp) ) 1107 1068 ALLOCATE ( seq_rand_xz(5,-mergp+nxl:nxr+mergp) ) … … 1109 1070 seq_rand_xz = 0 1110 1071 1111 CALL init_parallel_random_generator( nx, -mergp+nxl, nxr+mergp, & 1112 id_rand_xz, seq_rand_xz ) 1072 CALL init_parallel_random_generator( nx, -mergp+nxl, nxr+mergp, id_rand_xz, seq_rand_xz ) 1113 1073 ENDIF 1114 1074 … … 1118 1078 1119 1079 1120 !------------------------------------------------------------------------------ !1080 !--------------------------------------------------------------------------------------------------! 1121 1081 ! Description: 1122 1082 ! ------------ 1123 !> Calculate filter function bxx from length scale nxx following Eg.9 and 10 1124 !> (Xie and Castro, 2008) 1125 !------------------------------------------------------------------------------! 1083 !> Calculate filter function bxx from length scale nxx following Eg.9 and 10 (Xie and Castro, 2008) 1084 !--------------------------------------------------------------------------------------------------! 1126 1085 SUBROUTINE stg_filter_func( nxx, bxx ) 1127 1086 1128 INTEGER(iwp) :: k !< loop index 1129 INTEGER(iwp) :: n_k !< length scale nXX in height k 1130 INTEGER(iwp) :: nf !< index for length scales 1087 INTEGER(iwp) :: k !< loop index 1088 INTEGER(iwp) :: n_k !< length scale nXX in height k 1089 INTEGER(iwp) :: nf !< index for length scales 1090 1091 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: nxx !< length scale (in gp) 1131 1092 1132 1093 REAL(wp) :: bdenom !< denominator for filter functions bXX 1133 1094 REAL(wp) :: qsi = 1.0_wp !< minimization factor 1134 1135 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: nxx !< length scale (in gp)1136 1095 1137 1096 REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) :: bxx !< filter function … … 1148 1107 !-- ( Eq.10 )^2 1149 1108 DO nf = -n_k, n_k 1150 bdenom = bdenom + EXP( -qsi * pi * ABS( nf) / n_k )**21109 bdenom = bdenom + EXP( -qsi * pi * ABS( nf ) / n_k )**2 1151 1110 ENDDO 1152 1111 … … 1155 1114 bdenom = SQRT( bdenom ) 1156 1115 DO nf = -n_k, n_k 1157 bxx(nf,k) = EXP( -qsi * pi * ABS( nf) / n_k ) / bdenom1116 bxx(nf,k) = EXP( -qsi * pi * ABS( nf ) / n_k ) / bdenom 1158 1117 ENDDO 1159 1118 ENDIF … … 1163 1122 1164 1123 1165 !------------------------------------------------------------------------------ !1124 !--------------------------------------------------------------------------------------------------! 1166 1125 ! Description: 1167 1126 ! ------------ 1168 1127 !> Parin for &stg_par for synthetic turbulence generator 1169 !------------------------------------------------------------------------------ !1128 !--------------------------------------------------------------------------------------------------! 1170 1129 SUBROUTINE stg_parin 1171 1130 1172 CHARACTER (LEN=80) :: line 1173 1174 1175 NAMELIST /stg_par/ dt_stg_adjust, &1176 dt_stg_call, &1177 use_syn_turb_gen, &1131 CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file 1132 1133 1134 NAMELIST /stg_par/ dt_stg_adjust, & 1135 dt_stg_call, & 1136 use_syn_turb_gen, & 1178 1137 compute_velocity_seeds_local 1179 1138 … … 1184 1143 line = ' ' 1185 1144 DO WHILE ( INDEX( line, '&stg_par' ) == 0 ) 1186 READ ( 11, '(A)', END =20 ) line1145 READ ( 11, '(A)', END = 20 ) line 1187 1146 ENDDO 1188 1147 BACKSPACE ( 11 ) … … 1193 1152 1194 1153 ! 1195 !-- Set flag that indicates that the synthetic turbulence generator is switched 1196 !-- on 1154 !-- Set flag that indicates that the synthetic turbulence generator is switched on 1197 1155 syn_turb_gen = .TRUE. 1198 1156 GOTO 20 … … 1207 1165 1208 1166 1209 !------------------------------------------------------------------------------ !1167 !--------------------------------------------------------------------------------------------------! 1210 1168 ! Description: 1211 1169 ! ------------ 1212 1170 !> Read module-specific global restart data (Fortran binary format). 1213 !------------------------------------------------------------------------------ !1171 !--------------------------------------------------------------------------------------------------! 1214 1172 SUBROUTINE stg_rrd_global_ftn( found ) 1215 1173 1216 LOGICAL, INTENT(OUT) 1174 LOGICAL, INTENT(OUT) :: found !< flag indicating if variable was found 1217 1175 1218 1176 found = .TRUE. … … 1220 1178 1221 1179 SELECT CASE ( restart_string(1:length) ) 1222 1180 1223 1181 CASE ( 'time_stg_adjust' ) 1224 1182 READ ( 13 ) time_stg_adjust 1225 1183 1226 1184 CASE ( 'time_stg_call' ) 1227 1185 READ ( 13 ) time_stg_call 1228 1186 1229 1187 CASE ( 'use_syn_turb_gen' ) 1230 1188 READ ( 13 ) use_syn_turb_gen … … 1240 1198 1241 1199 1242 !------------------------------------------------------------------------------ !1200 !--------------------------------------------------------------------------------------------------! 1243 1201 ! Description: 1244 1202 ! ------------ 1245 1203 !> Read module-specific global restart data (MPI-IO). 1246 !------------------------------------------------------------------------------ !1204 !--------------------------------------------------------------------------------------------------! 1247 1205 SUBROUTINE stg_rrd_global_mpi 1248 1206 … … 1254 1212 1255 1213 1256 !------------------------------------------------------------------------------ !1214 !--------------------------------------------------------------------------------------------------! 1257 1215 ! Description: 1258 1216 ! ------------ 1259 1217 !> This routine writes the respective restart data. 1260 !------------------------------------------------------------------------------ !1218 !--------------------------------------------------------------------------------------------------! 1261 1219 SUBROUTINE stg_wrd_global 1262 1220 … … 1265 1223 CALL wrd_write_string( 'time_stg_adjust' ) 1266 1224 WRITE ( 14 ) time_stg_adjust 1267 1225 1268 1226 CALL wrd_write_string( 'time_stg_call' ) 1269 1227 WRITE ( 14 ) time_stg_call … … 1283 1241 1284 1242 1285 !------------------------------------------------------------------------------ !1243 !--------------------------------------------------------------------------------------------------! 1286 1244 ! Description: 1287 1245 ! ------------ 1288 !> Create turbulent inflow fields for u, v, w with prescribed length scales and 1289 !> Reynolds stress tensor after a method of Xie and Castro (2008), modified 1290 !> following suggestions of Kim et al. (2013), and using a Lund rotation 1291 !> (Lund, 1998). 1292 !------------------------------------------------------------------------------! 1246 !> Create turbulent inflow fields for u, v, w with prescribed length scales and Reynolds stress 1247 !> tensor after a method of Xie and Castro (2008), modified following suggestions of Kim et al. 1248 !> (2013), and using a Lund rotation (Lund, 1998). 1249 !--------------------------------------------------------------------------------------------------! 1293 1250 SUBROUTINE stg_main 1294 1251 1295 USE exchange_horiz_mod, &1252 USE exchange_horiz_mod, & 1296 1253 ONLY: exchange_horiz 1297 1254 1298 INTEGER(iwp) :: i!< grid index in x-direction1299 INTEGER(iwp) :: j!< loop index in y-direction1300 INTEGER(iwp) :: k!< loop index in z-direction1301 1302 LOGICAL :: stg_call = .FALSE.!< control flag indicating whether turbulence was updated or only restored from last call1303 1304 REAL(wp) :: dt_stg!< wheighted subtimestep1305 1306 REAL(wp), DIMENSION(3) :: mc_factor_l!< local mass flux correction factor1255 INTEGER(iwp) :: i !< grid index in x-direction 1256 INTEGER(iwp) :: j !< loop index in y-direction 1257 INTEGER(iwp) :: k !< loop index in z-direction 1258 1259 LOGICAL :: stg_call = .FALSE. !< control flag indicating whether turbulence was updated or only restored from last call 1260 1261 REAL(wp) :: dt_stg !< wheighted subtimestep 1262 1263 REAL(wp), DIMENSION(3) :: mc_factor_l !< local mass flux correction factor 1307 1264 1308 1265 IF ( debug_output_timestep ) CALL debug_message( 'stg_main', 'start' ) … … 1311 1268 dt_stg = MAX( dt_3d, dt_stg_call ) 1312 1269 ! 1313 !-- Check if synthetic turbulence generator needs to be called and new 1314 !-- perturbations need to be created or if old disturbances can be imposed 1315 !-- again. 1316 IF ( time_stg_call >= dt_stg_call .AND. & 1270 !-- Check if synthetic turbulence generator needs to be called and new perturbations need to be 1271 !-- created or if old disturbances can be imposed again. 1272 IF ( time_stg_call >= dt_stg_call .AND. & 1317 1273 intermediate_timestep_count == intermediate_timestep_count_max ) THEN 1318 1274 stg_call = .TRUE. … … 1324 1280 IF ( time_since_reference_point == 0.0_wp .AND. .NOT. velocity_seed_initialized ) THEN 1325 1281 ! 1326 !-- Generate turbulence at the left and right boundary. Random numbers 1327 !-- for the yz-planes at the left/right boundary are generated by the 1328 !-- left-sided mpi ranks only. After random numbers are calculated, they 1329 !-- are distributed to all other mpi ranks in the model, so that the 1330 !-- velocity seed matrices are available on all mpi ranks (i.e. also on the 1331 !-- right-sided boundary mpi ranks). In case of offline nesting, this implies, 1332 !-- that the left- and the right-sided lateral boundary have the same initial 1333 !-- seeds. 1334 !-- Note, in case of inflow from the right only, only turbulence at the left 1335 !-- boundary is required. 1336 IF ( .NOT. ( nesting_offline .OR. & 1337 ( child_domain .AND. rans_mode_parent & 1338 .AND. .NOT. rans_mode ) ) ) THEN 1282 !-- Generate turbulence at the left and right boundary. Random numbers for the yz-planes at the 1283 !-- left/right boundary are generated by the left-sided mpi ranks only. After random numbers are 1284 !-- calculated, they are distributed to all other mpi ranks in the model, so that the velocity 1285 !-- seed matrices are available on all mpi ranks (i.e. also on the right-sided boundary mpi ranks). 1286 !-- In case of offline nesting, this implies, that the left- and the right-sided lateral boundary 1287 !-- have the same initial seeds. 1288 !-- Note, in case of inflow from the right only, only turbulence at the left boundary is required. 1289 IF ( .NOT. ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent & 1290 .AND. .NOT. rans_mode ) ) ) THEN 1339 1291 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left ) 1340 1292 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left ) 1341 1293 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left ) 1342 1294 ELSE 1343 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, & 1344 id_stg_left, id_stg_right ) 1345 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, & 1346 id_stg_left, id_stg_right ) 1347 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, & 1348 id_stg_left, id_stg_right ) 1349 ! 1350 !-- Generate turbulence at the south and north boundary. Random numbers 1351 !-- for the xz-planes at the south/north boundary are generated by the 1352 !-- south-sided mpi ranks only. Please see also comment above. 1353 CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, & 1354 id_stg_south, id_stg_north ) 1355 CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, & 1356 id_stg_south, id_stg_north ) 1357 CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, & 1358 id_stg_south, id_stg_north ) 1295 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left, id_stg_right ) 1296 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left, id_stg_right ) 1297 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left, id_stg_right ) 1298 ! 1299 !-- Generate turbulence at the south and north boundary. Random numbers for the xz-planes at 1300 !-- the south/north boundary are generated by the south-sided mpi ranks only. Please see also 1301 !-- comment above. 1302 CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_south, id_stg_north ) 1303 CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_south, id_stg_north ) 1304 CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_south, id_stg_north ) 1359 1305 ENDIF 1360 1306 velocity_seed_initialized = .TRUE. 1361 1307 ENDIF 1362 1308 ! 1363 !-- New set of fu, fv, fw. Note, for inflow from the left side only, velocity 1364 !-- seeds are only required at the left boundary, while in case of offline 1365 !-- nesting or RANS-LES nesting velocity seeds are required also at the 1366 !-- right, south and north boundaries. 1309 !-- New set of fu, fv, fw. Note, for inflow from the left side only, velocity seeds are only 1310 !-- required at the left boundary, while in case of offline nesting or RANS-LES nesting velocity 1311 !-- seeds are required also at the right, south and north boundaries. 1367 1312 IF ( stg_call ) THEN 1368 IF ( .NOT. ( nesting_offline .OR. &1369 ( child_domain .AND. rans_mode_parent &1313 IF ( .NOT. ( nesting_offline .OR. & 1314 ( child_domain .AND. rans_mode_parent & 1370 1315 .AND. .NOT. rans_mode ) ) ) THEN 1371 1316 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left ) … … 1374 1319 1375 1320 ELSE 1376 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, & 1377 id_stg_left, id_stg_right ) 1378 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, & 1379 id_stg_left, id_stg_right ) 1380 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, & 1381 id_stg_left, id_stg_right ) 1382 1383 CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, & 1384 id_stg_south, id_stg_north ) 1385 CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, & 1386 id_stg_south, id_stg_north ) 1387 CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, & 1388 id_stg_south, id_stg_north ) 1321 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left, id_stg_right ) 1322 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left, id_stg_right ) 1323 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left, id_stg_right ) 1324 1325 CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south, id_stg_north ) 1326 CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south, id_stg_north ) 1327 CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south, id_stg_north ) 1389 1328 ENDIF 1390 1329 ENDIF 1391 1330 1392 1331 ! 1393 1332 !-- Turbulence generation at left and/or right boundary 1394 1333 IF ( myidx == id_stg_left .OR. myidx == id_stg_right ) THEN 1395 1334 ! 1396 !-- Calculate new set of perturbations. Do this only at last RK3-substep and 1397 !-- when dt_stg_call is exceeded. Else the old set of perturbations is 1398 !-- imposed 1335 !-- Calculate new set of perturbations. Do this only at last RK3-substep and when dt_stg_call is 1336 !-- exceeded. Else the old set of perturbations is imposed 1399 1337 IF ( stg_call ) THEN 1400 1338 1401 1339 DO j = nys, nyn 1402 1340 DO k = nzb, nzt + 1 1403 ! 1341 ! 1404 1342 !-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008) 1405 1343 IF ( tu(k) == 0.0_wp ) THEN 1406 1344 fu_yz(k,j) = fuo_yz(k,j) 1407 1345 ELSE 1408 fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &1409 fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )1346 fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + & 1347 fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) ) 1410 1348 ENDIF 1411 1349 1412 1350 IF ( tv(k) == 0.0_wp ) THEN 1413 1351 fv_yz(k,j) = fvo_yz(k,j) 1414 1352 ELSE 1415 fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &1416 fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )1353 fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + & 1354 fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) ) 1417 1355 ENDIF 1418 1356 1419 1357 IF ( tw(k) == 0.0_wp ) THEN 1420 1358 fw_yz(k,j) = fwo_yz(k,j) 1421 1359 ELSE 1422 fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &1423 fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )1360 fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + & 1361 fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) ) 1424 1362 ENDIF 1425 1363 ENDDO 1426 1364 ENDDO 1427 1365 1428 1366 dist_yz(nzb,:,1) = 0.0_wp 1429 1367 dist_yz(nzb,:,2) = 0.0_wp 1430 1368 dist_yz(nzb,:,3) = 0.0_wp 1431 1369 1432 1370 IF ( myidx == id_stg_left ) i = nxl 1433 IF ( myidx == id_stg_right ) i = nxr+1 1371 IF ( myidx == id_stg_right ) i = nxr+1 1434 1372 DO j = nys, nyn 1435 1373 DO k = nzb+1, nzt + 1 1436 ! 1374 ! 1437 1375 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1438 1376 !-- Additional factors are added to improve the variance of v and w 1439 dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) * & 1440 MERGE( 1.0_wp, 0.0_wp, & 1441 BTEST( wall_flags_total_0(k,j,i), 1 ) ) 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 ) ) 1442 1379 ENDDO 1443 1380 ENDDO … … 1447 1384 DO j = nys, nyn 1448 1385 DO k = nzb+1, nzt + 1 1449 ! 1386 ! 1450 1387 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1451 !-- Additional factors are added to improve the variance of v and w 1452 !-- experimental test of 1.2 1453 dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) ) & 1454 * 1.2_wp ) & 1455 * ( a21(k) * fu_yz(k,j) & 1456 + a22(k) * fv_yz(k,j) ), 3.0_wp ) * & 1457 MERGE( 1.0_wp, 0.0_wp, & 1458 BTEST( wall_flags_total_0(k,j,i), 2 ) ) 1459 dist_yz(k,j,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) ) & 1460 * 1.3_wp ) & 1461 * ( a31(k) * fu_yz(k,j) & 1462 + a32(k) * fv_yz(k,j) & 1463 + a33(k) * fw_yz(k,j) ), 3.0_wp ) * & 1464 MERGE( 1.0_wp, 0.0_wp, & 1388 !-- Additional factors are added to improve the variance of v and w experimental test 1389 !-- of 1.2 1390 dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL( a22 ) ) * 1.2_wp ) & 1391 * ( a21(k) * fu_yz(k,j) + a22(k) * fv_yz(k,j) ), 3.0_wp ) * & 1392 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) ) 1393 dist_yz(k,j,3) = MIN( ( SQRT( a33(k) / MAXVAL( a33 ) ) * 1.3_wp ) * & 1394 ( 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, & 1465 1396 BTEST( wall_flags_total_0(k,j,i), 3 ) ) 1466 1397 ENDDO … … 1469 1400 ! 1470 1401 !-- Mass flux correction following Kim et al. (2013) 1471 !-- This correction factor insures that the mass flux is preserved at the 1472 !-- inflow boundary. First, calculate mean value of the imposed 1473 !-- perturbations at yz boundary. 1474 !-- Note, this needs to be done only after the last RK3-substep, else 1475 !-- the boundary values won't be accessed. 1402 !-- This correction factor ensures that the mass flux is preserved at the inflow boundary. First, 1403 !-- calculate mean value of the imposed perturbations at yz boundary. Note, this needs to be done 1404 !-- only after the last RK3-substep, else the boundary values won't be accessed. 1476 1405 IF ( intermediate_timestep_count == intermediate_timestep_count_max ) THEN 1477 mc_factor_l = 0.0_wp 1478 mc_factor = 0.0_wp 1479 ! 1480 !-- Sum up the original volume flows (with and without perturbations). 1481 !-- Note, for non-normal components (here v and w) it is actually no 1482 !-- volume flow. 1406 mc_factor_l = 0.0_wp 1407 mc_factor = 0.0_wp 1408 ! 1409 !-- Sum up the original volume flows (with and without perturbations). 1410 !-- Note, for non-normal components (here v and w) it is actually no volume flow. 1483 1411 IF ( myidx == id_stg_left ) i = nxl 1484 1412 IF ( myidx == id_stg_right ) i = nxr+1 1485 1486 mc_factor_l(1) = SUM( dist_yz(nzb:nzt,nys:nyn,1) & 1487 * MERGE( 1.0_wp, 0.0_wp, & 1413 1414 mc_factor_l(1) = SUM( dist_yz(nzb:nzt,nys:nyn,1) * MERGE( 1.0_wp, 0.0_wp, & 1488 1415 BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 1 ) ) ) 1489 1416 1490 1417 IF ( myidx == id_stg_left ) i = nxl-1 1491 1418 IF ( myidx == id_stg_right ) i = nxr+1 1492 1493 mc_factor_l(2) = SUM( dist_yz(nzb:nzt,nysv:nyn,2) & 1494 * MERGE( 1.0_wp, 0.0_wp, & 1495 BTEST( wall_flags_total_0(nzb:nzt,nysv:nyn,i), 2 ) ) ) 1496 mc_factor_l(3) = SUM( dist_yz(nzb:nzt,nys:nyn,3) & 1497 * MERGE( 1.0_wp, 0.0_wp, & 1498 BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 3 ) ) ) 1499 1419 1420 mc_factor_l(2) = SUM( dist_yz(nzb:nzt,nysv:nyn,2) * MERGE( 1.0_wp, 0.0_wp, & 1421 BTEST( wall_flags_total_0(nzb:nzt,nysv:nyn,i), 2 ) ) ) 1422 mc_factor_l(3) = SUM( dist_yz(nzb:nzt,nys:nyn,3) * MERGE( 1.0_wp, 0.0_wp, & 1423 BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 3 ) ) ) 1424 1500 1425 #if defined( __parallel ) 1501 CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, & 1502 3, MPI_REAL, MPI_SUM, comm1dy, ierr ) 1503 #else 1504 mc_factor = mc_factor_l 1426 CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 3, MPI_REAL, MPI_SUM, comm1dy, ierr ) 1427 #else 1428 mc_factor = mc_factor_l 1505 1429 #endif 1506 ! 1507 !-- Calculate correction factor and force zero mean perturbations. 1508 mc_factor = mc_factor / REAL( nr_non_topo_yz, KIND = wp ) 1509 1510 IF ( myidx == id_stg_left ) i = nxl 1511 IF ( myidx == id_stg_right ) i = nxr+1 1512 1513 dist_yz(:,nys:nyn,1) = ( dist_yz(:,nys:nyn,1) - mc_factor(1) ) & 1514 * MERGE( 1.0_wp, 0.0_wp, & 1515 BTEST( wall_flags_total_0(:,nys:nyn,i), 1 ) ) 1516 1517 1518 IF ( myidx == id_stg_left ) i = nxl-1 1519 IF ( myidx == id_stg_right ) i = nxr+1 1520 1521 dist_yz(:,nys:nyn,2) = ( dist_yz(:,nys:nyn,2) - mc_factor(2) ) & 1522 * MERGE( 1.0_wp, 0.0_wp, & 1523 BTEST( wall_flags_total_0(:,nys:nyn,i), 2 ) ) 1524 1525 dist_yz(:,nys:nyn,3) = ( dist_yz(:,nys:nyn,3) - mc_factor(3) ) & 1526 * MERGE( 1.0_wp, 0.0_wp, & 1527 BTEST( wall_flags_total_0(:,nys:nyn,i), 3 ) ) 1528 ! 1430 ! 1431 !-- Calculate correction factor and force zero mean perturbations. 1432 mc_factor = mc_factor / REAL( nr_non_topo_yz, KIND = wp ) 1433 1434 IF ( myidx == id_stg_left ) i = nxl 1435 IF ( myidx == id_stg_right ) i = nxr+1 1436 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 ) ) 1439 1440 1441 IF ( myidx == id_stg_left ) i = nxl-1 1442 IF ( myidx == id_stg_right ) i = nxr+1 1443 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 ) ) 1449 ! 1529 1450 !-- Add disturbances 1530 1451 IF ( myidx == id_stg_left ) THEN 1531 ! 1532 !-- For the left boundary distinguish between mesoscale offline / self 1533 !-- nesting and turbulent inflow at the left boundary only. In the latter1534 !-- case turbulence is imposed on the mean inflow profiles.1535 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 1536 ! 1452 ! 1453 !-- For the left boundary distinguish between mesoscale offline / self nesting and 1454 !-- turbulent inflow at the left boundary only. In the latter case turbulence is imposed on 1455 !-- the mean inflow profiles. 1456 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 1457 ! 1537 1458 !-- Add disturbance at the inflow 1538 1459 DO j = nys, nyn 1539 1460 DO k = nzb, nzt+1 1540 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) + & 1541 dist_yz(k,j,1) ) & 1542 * MERGE( 1.0_wp, 0.0_wp, & 1543 BTEST( wall_flags_total_0(k,j,0), 1 ) ) 1544 v(k,j,-nbgp:-1) = ( mean_inflow_profiles(k,2) + & 1545 dist_yz(k,j,2) ) & 1546 * MERGE( 1.0_wp, 0.0_wp, & 1547 BTEST( wall_flags_total_0(k,j,-1), 2 ) ) 1548 w(k,j,-nbgp:-1) = dist_yz(k,j,3) & 1549 * MERGE( 1.0_wp, 0.0_wp, & 1550 BTEST( wall_flags_total_0(k,j,-1), 3 ) ) 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 BTEST( wall_flags_total_0(k,j,-1), 3 ) ) 1551 1467 ENDDO 1552 1468 ENDDO 1553 1469 ELSE 1554 1470 1555 1471 DO j = nys, nyn 1556 1472 DO k = nzb+1, nzt 1557 u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) ) & 1558 * MERGE( 1.0_wp, 0.0_wp, & 1473 u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) ) * MERGE( 1.0_wp, 0.0_wp, & 1559 1474 BTEST( wall_flags_total_0(k,j,0), 1 ) ) 1560 1475 u(k,j,-1) = u(k,j,0) 1561 v(k,j,-1) = ( v(k,j,-1) + dist_yz(k,j,2) ) & 1562 * MERGE( 1.0_wp, 0.0_wp, & 1476 v(k,j,-1) = ( v(k,j,-1) + dist_yz(k,j,2) ) * MERGE( 1.0_wp, 0.0_wp, & 1563 1477 BTEST( wall_flags_total_0(k,j,-1), 2 ) ) 1564 w(k,j,-1) = ( w(k,j,-1) + dist_yz(k,j,3) ) & 1565 * MERGE( 1.0_wp, 0.0_wp, & 1478 w(k,j,-1) = ( w(k,j,-1) + dist_yz(k,j,3) ) * MERGE( 1.0_wp, 0.0_wp, & 1566 1479 BTEST( wall_flags_total_0(k,j,-1), 3 ) ) 1567 1480 ENDDO … … 1572 1485 DO j = nys, nyn 1573 1486 DO k = nzb+1, nzt 1574 u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) ) & 1575 * MERGE( 1.0_wp, 0.0_wp, & 1576 BTEST( wall_flags_total_0(k,j,nxr+1), 1 ) ) 1577 v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) ) & 1578 * MERGE( 1.0_wp, 0.0_wp, & 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, & 1579 1490 BTEST( wall_flags_total_0(k,j,nxr+1), 2 ) ) 1580 w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) ) & 1581 * MERGE( 1.0_wp, 0.0_wp, & 1582 BTEST( wall_flags_total_0(k,j,nxr+1), 3 ) ) 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 ) ) 1583 1493 ENDDO 1584 1494 ENDDO … … 1590 1500 IF ( myidy == id_stg_north .OR. myidy == id_stg_south ) THEN 1591 1501 ! 1592 !-- Calculate new set of perturbations. Do this only at last RK3-substep and 1593 !-- when dt_stg_call is exceeded. Else the old set of perturbations is 1594 !-- imposed 1502 !-- Calculate new set of perturbations. Do this only at last RK3-substep and when dt_stg_call is 1503 !-- exceeded. Else the old set of perturbations is imposed 1595 1504 IF ( stg_call ) THEN 1596 1505 DO i = nxl, nxr 1597 1506 DO k = nzb, nzt + 1 1598 ! 1507 ! 1599 1508 !-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008) 1600 1509 IF ( tu(k) == 0.0_wp ) THEN 1601 1510 fu_xz(k,i) = fuo_xz(k,i) 1602 1511 ELSE 1603 fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &1604 fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )1512 fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + & 1513 fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) ) 1605 1514 ENDIF 1606 1515 1607 1516 IF ( tv(k) == 0.0_wp ) THEN 1608 1517 fv_xz(k,i) = fvo_xz(k,i) 1609 1518 ELSE 1610 fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &1611 fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )1519 fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + & 1520 fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) ) 1612 1521 ENDIF 1613 1522 1614 1523 IF ( tw(k) == 0.0_wp ) THEN 1615 1524 fw_xz(k,i) = fwo_xz(k,i) 1616 1525 ELSE 1617 fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &1618 fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )1526 fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + & 1527 fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) ) 1619 1528 ENDIF 1620 1529 ENDDO 1621 1530 ENDDO 1622 1623 1531 1532 1624 1533 dist_xz(nzb,:,1) = 0.0_wp 1625 1534 dist_xz(nzb,:,2) = 0.0_wp 1626 1535 dist_xz(nzb,:,3) = 0.0_wp 1627 1536 1628 1537 IF ( myidy == id_stg_south ) j = nys 1629 1538 IF ( myidy == id_stg_north ) j = nyn+1 1630 1539 DO i = nxl, nxr 1631 1540 DO k = nzb+1, nzt + 1 1632 ! 1541 ! 1633 1542 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1634 !-- Additional factors are added to improve the variance of v and w 1635 !experimental test of 1.2 1636 dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) ) & 1637 * 1.2_wp ) & 1638 * ( a21(k) * fu_xz(k,i) & 1639 + a22(k) * fv_xz(k,i) ), 3.0_wp ) * & 1640 MERGE( 1.0_wp, 0.0_wp, & 1641 BTEST( wall_flags_total_0(k,j,i), 2 ) ) 1543 !-- Additional factors are added to improve the variance of v and w 1544 !experimental test of 1.2 1545 dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL( a22 ) ) * 1.2_wp ) & 1546 * ( a21(k) * fu_xz(k,i) + a22(k) * fv_xz(k,i) ), 3.0_wp ) * & 1547 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) ) 1642 1548 ENDDO 1643 1549 ENDDO 1644 1550 1645 1551 IF ( myidy == id_stg_south ) j = nys-1 1646 1552 IF ( myidy == id_stg_north ) j = nyn+1 1647 1553 DO i = nxl, nxr 1648 1554 DO k = nzb+1, nzt + 1 1649 ! 1555 ! 1650 1556 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1651 1557 !-- Additional factors are added to improve the variance of v and w 1652 dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) * & 1653 MERGE( 1.0_wp, 0.0_wp, & 1654 BTEST( wall_flags_total_0(k,j,i), 1 ) ) 1655 1656 dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) ) & 1657 * 1.3_wp ) & 1658 * ( a31(k) * fu_xz(k,i) & 1659 + a32(k) * fv_xz(k,i) & 1660 + a33(k) * fw_xz(k,i) ), 3.0_wp ) * & 1661 MERGE( 1.0_wp, 0.0_wp, & 1662 BTEST( wall_flags_total_0(k,j,i), 3 ) ) 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 ) ) 1560 1561 dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL( a33 ) ) * 1.3_wp ) & 1562 * ( 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 ) ) 1663 1565 ENDDO 1664 1566 ENDDO … … 1666 1568 1667 1569 ! 1668 !-- Mass flux correction following Kim et al. (2013) 1669 !-- This correction factor insures that the mass flux is preserved at the 1670 !-- inflow boundary. First, calculate mean value of the imposed 1671 !-- perturbations at yz boundary. 1672 !-- Note, this needs to be done only after the last RK3-substep, else 1673 !-- the boundary values won't be accessed. 1570 !-- Mass flux correction following Kim et al. (2013). This correction factor ensures that the 1571 !-- mass flux is preserved at the inflow boundary. First, calculate mean value of the imposed 1572 !-- perturbations at yz boundary. Note, this needs to be done only after the last RK3-substep, 1573 !-- else the boundary values won't be accessed. 1674 1574 IF ( intermediate_timestep_count == intermediate_timestep_count_max ) THEN 1675 mc_factor_l 1676 mc_factor 1677 1575 mc_factor_l = 0.0_wp 1576 mc_factor = 0.0_wp 1577 1678 1578 IF ( myidy == id_stg_south ) j = nys 1679 1579 IF ( myidy == id_stg_north ) j = nyn+1 1680 1681 mc_factor_l(2) = SUM( dist_xz(nzb:nzt,nxl:nxr,2) & 1682 * MERGE( 1.0_wp, 0.0_wp, & 1580 1581 mc_factor_l(2) = SUM( dist_xz(nzb:nzt,nxl:nxr,2) * MERGE( 1.0_wp, 0.0_wp, & 1683 1582 BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 2 ) ) ) 1684 1685 IF ( myidy == id_stg_south ) j = nys-1 1686 IF ( myidy == id_stg_north ) j = nyn+1 1687 1688 mc_factor_l(1) = SUM( dist_xz(nzb:nzt,nxlu:nxr,1) & 1689 * MERGE( 1.0_wp, 0.0_wp, & 1583 1584 IF ( myidy == id_stg_south ) j = nys-1 1585 IF ( myidy == id_stg_north ) j = nyn+1 1586 1587 mc_factor_l(1) = SUM( dist_xz(nzb:nzt,nxlu:nxr,1) * MERGE( 1.0_wp, 0.0_wp, & 1690 1588 BTEST( wall_flags_total_0(nzb:nzt,j,nxlu:nxr), 1 ) ) ) 1691 mc_factor_l(3) = SUM( dist_xz(nzb:nzt,nxl:nxr,3) & 1692 * MERGE( 1.0_wp, 0.0_wp, & 1589 mc_factor_l(3) = SUM( dist_xz(nzb:nzt,nxl:nxr,3) * MERGE( 1.0_wp, 0.0_wp, & 1693 1590 BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 3 ) ) ) 1694 1591 1695 1592 #if defined( __parallel ) 1696 CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, & 1697 3, MPI_REAL, MPI_SUM, comm1dx, ierr ) 1698 #else 1593 CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 3, MPI_REAL, MPI_SUM, comm1dx, ierr ) 1594 #else 1699 1595 mc_factor = mc_factor_l 1700 1596 #endif 1701 1597 1702 1598 mc_factor = mc_factor / REAL( nr_non_topo_xz, KIND = wp ) 1703 1704 IF ( myidy == id_stg_south ) j = nys 1705 IF ( myidy == id_stg_north ) j = nyn+1 1706 1707 dist_xz(:,nxl:nxr,2) = ( dist_xz(:,nxl:nxr,2) - mc_factor(2) ) & 1708 * MERGE( 1.0_wp, 0.0_wp, & 1709 BTEST( wall_flags_total_0(:,j,nxl:nxr), 2 ) ) 1710 1711 1712 IF ( myidy == id_stg_south ) j = nys-1 1713 IF ( myidy == id_stg_north ) j = nyn+1 1714 1715 dist_xz(:,nxl:nxr,1) = ( dist_xz(:,nxl:nxr,1) - mc_factor(1) ) & 1716 * MERGE( 1.0_wp, 0.0_wp, & 1717 BTEST( wall_flags_total_0(:,j,nxl:nxr), 1 ) ) 1718 1719 dist_xz(:,nxl:nxr,3) = ( dist_xz(:,nxl:nxr,3) - mc_factor(3) ) & 1720 * MERGE( 1.0_wp, 0.0_wp, & 1721 BTEST( wall_flags_total_0(:,j,nxl:nxr), 3 ) ) 1722 ! 1599 1600 IF ( myidy == id_stg_south ) j = nys 1601 IF ( myidy == id_stg_north ) j = nyn+1 1602 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 ) ) 1605 1606 1607 IF ( myidy == id_stg_south ) j = nys-1 1608 IF ( myidy == id_stg_north ) j = nyn+1 1609 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 ) ) 1615 ! 1723 1616 !-- Add disturbances 1724 IF ( myidy == id_stg_south 1617 IF ( myidy == id_stg_south ) THEN 1725 1618 DO i = nxl, nxr 1726 1619 DO k = nzb+1, nzt 1727 u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) ) & 1728 * MERGE( 1.0_wp, 0.0_wp, & 1729 BTEST( wall_flags_total_0(k,-1,i), 1 ) ) 1730 v(k,0,i) = ( v(k,0,i) + dist_xz(k,i,2) ) & 1731 * MERGE( 1.0_wp, 0.0_wp, & 1732 BTEST( wall_flags_total_0(k,0,i), 2 ) ) 1620 u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) ) & 1621 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,-1,i), 1 ) ) 1622 v(k,0,i) = ( v(k,0,i) + dist_xz(k,i,2) ) & 1623 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,0,i), 2 ) ) 1733 1624 v(k,-1,i) = v(k,0,i) 1734 w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3) ) & 1735 * MERGE( 1.0_wp, 0.0_wp, & 1736 BTEST( wall_flags_total_0(k,-1,i), 3 ) ) 1625 w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3) ) & 1626 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,-1,i), 3 ) ) 1737 1627 ENDDO 1738 1628 ENDDO 1739 1629 ENDIF 1740 1630 IF ( myidy == id_stg_north ) THEN 1741 1631 1742 1632 DO i = nxl, nxr 1743 1633 DO k = nzb+1, nzt 1744 u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) ) & 1745 * MERGE( 1.0_wp, 0.0_wp, & 1746 BTEST( wall_flags_total_0(k,nyn+1,i), 1 ) ) 1747 v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) ) & 1748 * MERGE( 1.0_wp, 0.0_wp, & 1749 BTEST( wall_flags_total_0(k,nyn+1,i), 2 ) ) 1750 w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) ) & 1751 * MERGE( 1.0_wp, 0.0_wp, & 1752 BTEST( wall_flags_total_0(k,nyn+1,i), 3 ) ) 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 ) ) 1753 1640 ENDDO 1754 1641 ENDDO … … 1769 1656 END SUBROUTINE stg_main 1770 1657 1771 !------------------------------------------------------------------------------ !1658 !--------------------------------------------------------------------------------------------------! 1772 1659 ! Description: 1773 1660 ! ------------ 1774 !> Generate a set of random number rand_it wich is equal on each PE 1775 !> and calculate the velocity seed f_n. 1776 !> f_n is splitted in vertical direction by the number of PEs in x-direction and 1777 !> and each PE calculates a vertical subsection of f_n. At the the end, all 1778 !> parts are collected to form the full array. 1779 !------------------------------------------------------------------------------! 1661 !> Generate a set of random number rand_it wich is equal on each PE and calculate the velocity seed 1662 !> f_n. f_n is splitted in vertical direction by the number of PEs in x-direction and each PE 1663 !> calculates a vertical subsection of f_n. At the the end, all parts are collected to form the full 1664 !> array. 1665 !--------------------------------------------------------------------------------------------------! 1780 1666 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id_left, id_right ) 1781 1667 1782 INTEGER(iwp) :: id_left !< core ids at respective boundaries 1783 INTEGER(iwp), OPTIONAL :: id_right !< core ids at respective boundaries 1784 INTEGER(iwp) :: j !< loop index in y-direction 1785 INTEGER(iwp) :: jj !< loop index in y-direction 1786 INTEGER(iwp) :: k !< loop index in z-direction 1787 INTEGER(iwp) :: kk !< loop index in z-direction 1788 INTEGER(iwp) :: send_count !< send count for MPI_GATHERV 1789 1790 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y !< length scale in y-direction 1791 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z !< length scale in z-direction 1792 1793 REAL(wp) :: nyz_inv !< inverse of number of grid points in yz-slice 1794 REAL(wp) :: rand_av !< average of random number 1795 REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number 1796 1797 REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) :: b_y !< filter function in y-direction 1798 REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) :: b_z !< filter function in z-direction 1799 1800 REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nys:nyn) :: f_n_l !< local velocity seed 1801 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn) :: f_n !< velocity seed 1802 1803 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rand_it !< global array of random numbers 1804 ! 1805 !-- Generate random numbers using the parallel random generator. 1806 !-- The set of random numbers are modified to have an average of 0 and 1807 !-- unit variance. Note, at the end the random number array must be defined 1808 !-- globally in order to compute the correlation matrices. However, 1809 !-- the calculation and normalization of random numbers is done locally, 1810 !-- while the result is later distributed to all processes. Further, 1811 !-- please note, a set of random numbers is only calculated for the 1812 !-- left boundary, while the right boundary uses the same random numbers 1813 !-- and thus also computes the same correlation matrix. 1668 INTEGER(iwp) :: id_left !< core ids at respective boundaries 1669 INTEGER(iwp), OPTIONAL :: id_right !< core ids at respective boundaries 1670 INTEGER(iwp) :: j !< loop index in y-direction 1671 INTEGER(iwp) :: jj !< loop index in y-direction 1672 INTEGER(iwp) :: k !< loop index in z-direction 1673 INTEGER(iwp) :: kk !< loop index in z-direction 1674 INTEGER(iwp) :: send_count !< send count for MPI_GATHERV 1675 1676 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y !< length scale in y-direction 1677 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z !< length scale in z-direction 1678 1679 REAL(wp) :: nyz_inv !< inverse of number of grid points in yz-slice 1680 REAL(wp) :: rand_av !< average of random number 1681 REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number 1682 1683 REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) :: b_y !< filter function in y-direction 1684 REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) :: b_z !< filter function in z-direction 1685 1686 REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nys:nyn) :: f_n_l !< local velocity seed 1687 REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn) :: f_n !< velocity seed 1688 1689 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rand_it !< global array of random numbers 1690 ! 1691 !-- Generate random numbers using the parallel random generator. The set of random numbers are 1692 !-- modified to have an average of 0 and unit variance. Note, at the end the random number array 1693 !-- must be defined globally in order to compute the correlation matrices. However, the calculation 1694 !-- and normalization of random numbers is done locally, while the result is later distributed to 1695 !-- all processes. Further, please note, a set of random numbers is only calculated for the left 1696 !-- boundary, while the right boundary uses the same random numbers and thus also computes the same 1697 !-- correlation matrix. 1814 1698 ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp+nys:nyn+mergp) ) 1815 1699 rand_it = 0.0_wp … … 1817 1701 rand_av = 0.0_wp 1818 1702 rand_sigma_inv = 0.0_wp 1819 nyz_inv = 1.0_wp / REAL( ( nzt + 1 + mergp - ( nzb - mergp ) + 1 ) & 1820 * ( ny + mergp - ( 0 - mergp ) + 1 ), & 1821 KIND=wp ) 1703 nyz_inv = 1.0_wp / REAL( ( nzt + 1 + mergp - ( nzb - mergp ) + 1 ) & 1704 * ( ny + mergp - ( 0 - mergp ) + 1 ), KIND = wp ) 1822 1705 ! 1823 1706 !-- Compute and normalize random numbers. … … 1835 1718 ENDDO 1836 1719 ! 1837 !-- For normalization to zero mean, sum-up the global random numers. 1838 !-- To normalize the global set of random numbers, 1839 !-- the inner ghost layers mergp must not be summed-up, else 1840 !-- the random numbers on the ghost layers will be stronger weighted as they 1841 !-- also occur on the inner subdomains. 1842 DO j = MERGE( nys, nys - mergp, nys /= 0 ), & 1843 MERGE( nyn, nyn + mergp, nyn /= ny ) 1720 !-- For normalization to zero mean, sum-up the global random numers. To normalize the global set of 1721 !-- random numbers, the inner ghost layers mergp must not be summed-up, else the random numbers on 1722 !-- the ghost layers will be stronger weighted as they also occur on the inner subdomains. 1723 DO j = MERGE( nys, nys - mergp, nys /= 0 ), MERGE( nyn, nyn + mergp, nyn /= ny ) 1844 1724 DO k = nzb - mergp, nzt + 1 + mergp 1845 1725 rand_av = rand_av + rand_it(k,j) 1846 1726 ENDDO 1847 1727 ENDDO 1848 1728 1849 1729 #if defined( __parallel ) 1850 1730 ! 1851 1731 !-- Sum-up the local averages of the random numbers 1852 CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL, & 1853 MPI_SUM, comm1dy, ierr ) 1732 CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL, MPI_SUM, comm1dy, ierr ) 1854 1733 #endif 1855 1734 rand_av = rand_av * nyz_inv … … 1859 1738 ! 1860 1739 !-- Now, compute the variance 1861 DO j = MERGE( nys, nys - mergp, nys /= 0 ), & 1862 MERGE( nyn, nyn + mergp, nyn /= ny ) 1740 DO j = MERGE( nys, nys - mergp, nys /= 0 ), MERGE( nyn, nyn + mergp, nyn /= ny ) 1863 1741 DO k = nzb - mergp, nzt + 1 + mergp 1864 1742 rand_sigma_inv = rand_sigma_inv + rand_it(k,j)**2 … … 1869 1747 ! 1870 1748 !-- Sum-up the local quadratic averages of the random numbers 1871 CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL, & 1872 MPI_SUM, comm1dy, ierr ) 1749 CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL, MPI_SUM, comm1dy, ierr ) 1873 1750 #endif 1874 1751 ! … … 1885 1762 CALL cpu_log( log_point_s(31), 'STG f_n factors', 'start' ) 1886 1763 ! 1887 !-- Generate velocity seed following Eq.6 of Xie and Castro (2008). There 1888 !-- are two options. In the first one, the computation of the seeds is 1889 !-- distributed to all processes along the communicator comm1dy and 1890 !-- gathered on the leftmost and, if necessary, on the rightmost process. 1891 !-- For huge length scales the computational effort can become quite huge 1892 !-- (it scales with the turbulent length scales), so that gain by parallelization 1893 !-- exceeds the costs by the subsequent communication. 1894 !-- In the second option, which performs better when the turbulent length scales 1895 !-- are parametrized and thus the loops are smaller, the seeds are computed 1896 !-- locally and no communication is necessary. 1764 !-- Generate velocity seed following Eq.6 of Xie and Castro (2008). There are two options. In the 1765 !-- first one, the computation of the seeds is distributed to all processes along the communicator 1766 !-- comm1dy and gathered on the leftmost and, if necessary, on the rightmost process. For huge 1767 !-- length scales the computational effort can become quite huge (it scales with the turbulent 1768 !-- length scales), so that gain by parallelization exceeds the costs by the subsequent 1769 !-- communication. In the second option, which performs better when the turbulent length scales 1770 !-- are parametrized and thus the loops are smaller, the seeds are computed locally and no 1771 !-- communication is necessary. 1897 1772 IF ( compute_velocity_seeds_local ) THEN 1898 1773 … … 1928 1803 ! 1929 1804 !-- Gather the velocity seed matrix on left boundary mpi ranks. 1930 CALL MPI_GATHERV( f_n_l(nzb_x_stg,nys), send_count, stg_type_yz_small, & 1931 f_n(nzb+1,nys), recv_count_yz, displs_yz, stg_type_yz,& 1932 id_left, comm1dx, ierr ) 1933 ! 1934 !-- If required, gather the same velocity seed matrix on right boundary 1935 !-- mpi ranks (in offline nesting for example). 1805 CALL MPI_GATHERV( f_n_l(nzb_x_stg,nys), send_count, stg_type_yz_small, f_n(nzb+1,nys), & 1806 recv_count_yz, displs_yz, stg_type_yz, id_left, comm1dx, ierr ) 1807 ! 1808 !-- If required, gather the same velocity seed matrix on right boundary mpi ranks (in offline 1809 !-- nesting for example). 1936 1810 IF ( PRESENT( id_right ) ) THEN 1937 CALL MPI_GATHERV( f_n_l(nzb_x_stg,nys), send_count, stg_type_yz_small, & 1938 f_n(nzb+1,nys), recv_count_yz, displs_yz, stg_type_yz,& 1939 id_right, comm1dx, ierr ) 1811 CALL MPI_GATHERV( f_n_l(nzb_x_stg,nys), send_count, stg_type_yz_small, f_n(nzb+1,nys), & 1812 recv_count_yz, displs_yz, stg_type_yz, id_right, comm1dx, ierr ) 1940 1813 ENDIF 1941 1814 #else … … 1955 1828 1956 1829 1957 !------------------------------------------------------------------------------ !1830 !--------------------------------------------------------------------------------------------------! 1958 1831 ! Description: 1959 1832 ! ------------ 1960 !> Generate a set of random number rand_it wich is equal on each PE 1961 !> and calculate the velocity seedf_n.1962 !> f_n is splitted in vertical direction by the number of PEs in y-direction and 1963 !> and each PE calculates a vertical subsection of f_n. At the the end, all1964 !> parts are collected to form thefull array.1965 !------------------------------------------------------------------------------ !1833 !> Generate a set of random number rand_it wich is equal on each PE and calculate the velocity seed 1834 !> f_n. 1835 !> f_n is splitted in vertical direction by the number of PEs in y-direction and and each PE 1836 !> calculates a vertical subsection of f_n. At the the end, all parts are collected to form the 1837 !> full array. 1838 !--------------------------------------------------------------------------------------------------! 1966 1839 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id_south, id_north ) 1967 1840 1968 INTEGER(iwp) :: i !< loop index in x-direction 1969 INTEGER(iwp) :: id_north !< core ids at respective boundaries 1970 INTEGER(iwp) :: id_south !< core ids at respective boundaries 1971 INTEGER(iwp) :: ii !< loop index in x-direction 1972 INTEGER(iwp) :: k !< loop index in z-direction 1973 INTEGER(iwp) :: kk !< loop index in z-direction 1974 INTEGER(iwp) :: send_count !< send count for MPI_GATHERV 1975 1976 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x !< length scale in x-direction 1977 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z !< length scale in z-direction 1978 1979 REAL(wp) :: nxz_inv !< inverse of number of grid points in xz-slice 1980 REAL(wp) :: rand_av !< average of random number 1981 REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number 1982 1983 REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) :: b_x !< filter function in x-direction 1984 REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) :: b_z !< filter function in z-direction 1985 1986 REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxl:nxr) :: f_n_l !< local velocity seed 1987 REAL(wp), DIMENSION(nzb:nzt+1,nxl:nxr) :: f_n !< velocity seed 1988 1989 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rand_it !< global array of random numbers 1990 1991 ! 1992 !-- Generate random numbers using the parallel random generator. 1993 !-- The set of random numbers are modified to have an average of 0 and 1994 !-- unit variance. Note, at the end the random number array must be defined 1995 !-- globally in order to compute the correlation matrices. However, 1996 !-- the calculation and normalization of random numbers is done locally, 1997 !-- while the result is later distributed to all processes. Further, 1998 !-- please note, a set of random numbers is only calculated for the 1999 !-- left boundary, while the right boundary uses the same random numbers 2000 !-- and thus also computes the same correlation matrix. 1841 INTEGER(iwp) :: i !< loop index in x-direction 1842 INTEGER(iwp) :: id_north !< core ids at respective boundaries 1843 INTEGER(iwp) :: id_south !< core ids at respective boundaries 1844 INTEGER(iwp) :: ii !< loop index in x-direction 1845 INTEGER(iwp) :: k !< loop index in z-direction 1846 INTEGER(iwp) :: kk !< loop index in z-direction 1847 INTEGER(iwp) :: send_count !< send count for MPI_GATHERV 1848 1849 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x !< length scale in x-direction 1850 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z !< length scale in z-direction 1851 1852 REAL(wp) :: nxz_inv !< inverse of number of grid points in xz-slice 1853 REAL(wp) :: rand_av !< average of random number 1854 REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number 1855 1856 REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) :: b_x !< filter function in x-direction 1857 REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) :: b_z !< filter function in z-direction 1858 1859 REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxl:nxr) :: f_n_l !< local velocity seed 1860 REAL(wp), DIMENSION(nzb:nzt+1,nxl:nxr) :: f_n !< velocity seed 1861 1862 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rand_it !< global array of random numbers 1863 1864 ! 1865 !-- Generate random numbers using the parallel random generator. The set of random numbers are 1866 !-- modified to have an average of 0 and unit variance. Note, at the end the random number array 1867 !-- must be defined globally in order to compute the correlation matrices. However, the calculation 1868 !-- and normalization of random numbers is done locally, while the result is later distributed to 1869 !-- all processes. Further, please note, a set of random numbers is only calculated for the left 1870 !-- boundary, while the right boundary uses the same random numbers and thus also computes the same 1871 !-- correlation matrix. 2001 1872 ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp+nxl:nxr+mergp) ) 2002 1873 rand_it = 0.0_wp … … 2004 1875 rand_av = 0.0_wp 2005 1876 rand_sigma_inv = 0.0_wp 2006 nxz_inv = 1.0_wp / REAL( ( nzt + 1 + mergp - ( nzb - mergp ) + 1 ) & 2007 * ( nx + mergp - ( 0 - mergp ) +1 ), & 2008 KIND=wp ) 1877 nxz_inv = 1.0_wp / REAL( ( nzt + 1 + mergp - ( nzb - mergp ) + 1 ) & 1878 * ( nx + mergp - ( 0 - mergp ) +1 ), KIND = wp ) 2009 1879 ! 2010 1880 !-- Compute and normalize random numbers. … … 2023 1893 ! 2024 1894 !-- For normalization to zero mean, sum-up the global random numers. 2025 !-- To normalize the global set of random numbers, 2026 !-- the inner ghost layers mergp must not be summed-up, else 2027 !-- the random numbers on the ghost layers will be stronger weighted as they 1895 !-- To normalize the global set of random numbers, the inner ghost layers mergp must not be 1896 !-- summed-up, else the random numbers on the ghost layers will be stronger weighted as they 2028 1897 !-- also occur on the inner subdomains. 2029 DO i = MERGE( nxl, nxl - mergp, nxl /= 0 ), & 2030 MERGE( nxr, nxr + mergp, nxr /= nx ) 1898 DO i = MERGE( nxl, nxl - mergp, nxl /= 0 ), MERGE( nxr, nxr + mergp, nxr /= nx ) 2031 1899 DO k = nzb - mergp, nzt + 1 + mergp 2032 1900 rand_av = rand_av + rand_it(k,i) 2033 1901 ENDDO 2034 1902 ENDDO 2035 1903 2036 1904 #if defined( __parallel ) 2037 1905 ! 2038 1906 !-- Sum-up the local averages of the random numbers 2039 CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL, & 2040 MPI_SUM, comm1dx, ierr ) 1907 CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL, MPI_SUM, comm1dx, ierr ) 2041 1908 #endif 2042 1909 rand_av = rand_av * nxz_inv … … 2046 1913 ! 2047 1914 !-- Now, compute the variance 2048 DO i = MERGE( nxl, nxl - mergp, nxl /= 0 ), & 2049 MERGE( nxr, nxr + mergp, nxr /= nx ) 1915 DO i = MERGE( nxl, nxl - mergp, nxl /= 0 ), MERGE( nxr, nxr + mergp, nxr /= nx ) 2050 1916 DO k = nzb - mergp, nzt + 1 + mergp 2051 1917 rand_sigma_inv = rand_sigma_inv + rand_it(k,i)**2 … … 2056 1922 ! 2057 1923 !-- Sum-up the local quadratic averages of the random numbers 2058 CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL, & 2059 MPI_SUM, comm1dx, ierr ) 1924 CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL, MPI_SUM, comm1dx, ierr ) 2060 1925 #endif 2061 1926 ! … … 2072 1937 CALL cpu_log( log_point_s(31), 'STG f_n factors', 'start' ) 2073 1938 ! 2074 !-- Generate velocity seed following Eq.6 of Xie and Castro (2008). There 2075 !-- are two options. In the first one, the computation of the seeds is 2076 !-- distributed to all processes along the communicator comm1dx and 2077 !-- gathered on the southmost and, if necessary, on the northmost process. 2078 !-- For huge length scales the computational effort can become quite huge 2079 !-- (it scales with the turbulent length scales), so that gain by parallelization 2080 !-- exceeds the costs by the subsequent communication. 2081 !-- In the second option, which performs better when the turbulent length scales 2082 !-- are parametrized and thus the loops are smaller, the seeds are computed 2083 !-- locally and no communication is necessary. 1939 !-- Generate velocity seed following Eq.6 of Xie and Castro (2008). There are two options. In the 1940 !-- first one, the computation of the seeds is distributed to all processes along the communicator 1941 !-- comm1dx and gathered on the southmost and, if necessary, on the northmost process. For huge 1942 !-- length scales the computational effort can become quite huge (it scales with the turbulent 1943 !-- length scales), so that gain by parallelization exceeds the costs by the subsequent communication. 1944 !-- In the second option, which performs better when the turbulent length scales are parametrized 1945 !-- and thus the loops are smaller, the seeds are computed locally and no communication is necessary. 2084 1946 IF ( compute_velocity_seeds_local ) THEN 2085 1947 … … 2115 1977 ! 2116 1978 !-- Gather the processed velocity seed on south boundary mpi ranks. 2117 CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxl), send_count, stg_type_xz_small, & 2118 f_n(nzb+1,nxl), recv_count_xz, displs_xz, stg_type_xz, & 2119 id_south, comm1dy, ierr ) 1979 CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxl), send_count, stg_type_xz_small, f_n(nzb+1,nxl), & 1980 recv_count_xz, displs_xz, stg_type_xz, id_south, comm1dy, ierr ) 2120 1981 ! 2121 1982 !-- Gather the processed velocity seed on north boundary mpi ranks. 2122 CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxl), send_count, stg_type_xz_small, & 2123 f_n(nzb+1,nxl), recv_count_xz, displs_xz, stg_type_xz, & 2124 id_north, comm1dy, ierr ) 1983 CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxl), send_count, stg_type_xz_small, f_n(nzb+1,nxl), & 1984 recv_count_xz, displs_xz, stg_type_xz, id_north, comm1dy, ierr ) 2125 1985 #else 2126 1986 f_n(nzb+1:nzt+1,nxl:nxr) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxl:nxr) … … 2138 1998 END SUBROUTINE stg_generate_seed_xz 2139 1999 2140 !------------------------------------------------------------------------------ !2000 !--------------------------------------------------------------------------------------------------! 2141 2001 ! Description: 2142 2002 ! ------------ 2143 !> Parametrization of the Reynolds stress tensor, following the parametrization 2144 !> described in Rotach et al. (1996), which is applied in state-of-the-art 2145 !> dispserion modelling. Please note, the parametrization does not distinguish 2146 !> between along-wind and cross-wind turbulence. 2147 !------------------------------------------------------------------------------! 2003 !> Parametrization of the Reynolds stress tensor, following the parametrization described in Rotach 2004 !> et al. (1996), which is applied in state-of-the-art dispserion modelling. Please note, the 2005 !> parametrization does not distinguish between along-wind and cross-wind turbulence. 2006 !--------------------------------------------------------------------------------------------------! 2148 2007 SUBROUTINE parametrize_reynolds_stress 2149 2008 2150 INTEGER(iwp) :: k!< loop index in z-direction2151 2152 REAL(wp) :: zzi 2153 2009 INTEGER(iwp) :: k !< loop index in z-direction 2010 2011 REAL(wp) :: zzi !< ratio of z/zi 2012 2154 2013 DO k = nzb+1, nzt+1 2155 2014 ! 2156 2015 !-- Calculate function to gradually decrease Reynolds stress above ABL top. 2157 !-- The function decreases to 1/10 after one length scale above the 2158 !-- ABL top. 2016 !-- The function decreases to 1/10 after one length scale above the ABL top. 2159 2017 blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) ) 2160 2018 ! … … 2162 2020 zzi = zu(k) / zi_ribulk 2163 2021 ! 2164 !-- u'u' and v'v'. Assume isotropy. Note, add a small negative number 2165 !-- to the denominator, else the mergpe-function can crash if scale_l is 2166 !-- zero. 2167 r11(k) = scale_us**2 * ( & 2168 MERGE( 0.35_wp * ABS( & 2169 - zi_ribulk / ( kappa * scale_l - 10E-4_wp ) & 2170 )**( 2.0_wp / 3.0_wp ), & 2171 0.0_wp, & 2172 scale_l < 0.0_wp ) & 2173 + 5.0_wp - 4.0_wp * zzi & 2174 ) * blend 2175 2176 r22(k) = r11(k) 2177 ! 2178 !-- w'w' 2179 r33(k) = scale_wm**2 * ( & 2180 1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi ) & 2181 + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2 & 2182 ) * blend 2183 ! 2184 !-- u'w' and v'w'. Assume isotropy. 2022 !-- u'u' and v'v'. Assume isotropy. Note, add a small negative number to the denominator, else 2023 !-- 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 2028 r22(k) = r11(k) 2029 ! 2030 !-- w'w' 2031 r33(k) = scale_wm**2 * ( 1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi ) & 2032 + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2 ) * blend 2033 ! 2034 !-- u'w' and v'w'. Assume isotropy. 2185 2035 r31(k) = - scale_us**2 * ( 1.01_wp - MIN( zzi, 1.0_wp ) ) * blend 2186 2036 r32(k) = r31(k) 2187 2037 ! 2188 !-- For u'v' no parametrization exist so far - ?. For simplicity assume 2189 !-- a similar profile as for u'w'.2038 !-- For u'v' no parametrization exist so far - ?. For simplicity assume a similar profile as for 2039 !-- u'w'. 2190 2040 r21(k) = r31(k) 2191 2041 ENDDO 2192 2042 2193 2043 ! 2194 !-- Set bottom boundary condition 2044 !-- Set bottom boundary condition 2195 2045 r11(nzb) = r11(nzb+1) 2196 2046 r22(nzb) = r22(nzb+1) … … 2200 2050 r31(nzb) = r31(nzb+1) 2201 2051 r32(nzb) = r32(nzb+1) 2202 2203 2204 END SUBROUTINE parametrize_reynolds_stress 2205 2206 !------------------------------------------------------------------------------ !2052 2053 2054 END SUBROUTINE parametrize_reynolds_stress 2055 2056 !--------------------------------------------------------------------------------------------------! 2207 2057 ! Description: 2208 2058 ! ------------ 2209 !> Calculate the coefficient matrix from the Lund rotation. 2210 !------------------------------------------------------------------------------ !2059 !> Calculate the coefficient matrix from the Lund rotation. 2060 !--------------------------------------------------------------------------------------------------! 2211 2061 SUBROUTINE calc_coeff_matrix 2212 2062 2213 INTEGER(iwp) :: k!< loop index in z-direction2214 2063 INTEGER(iwp) :: k !< loop index in z-direction 2064 2215 2065 ! 2216 2066 !-- Calculate coefficient matrix. Split loops to allow for loop vectorization. 2217 2067 DO k = nzb+1, nzt+1 2218 2068 IF ( r11(k) > 10E-6_wp ) THEN 2219 a11(k) = SQRT( r11(k) ) 2069 a11(k) = SQRT( r11(k) ) 2220 2070 a21(k) = r21(k) / a11(k) 2221 2071 a31(k) = r31(k) / a11(k) … … 2227 2077 ENDDO 2228 2078 DO k = nzb+1, nzt+1 2229 a22(k) = r22(k) - a21(k)**2 2079 a22(k) = r22(k) - a21(k)**2 2230 2080 IF ( a22(k) > 10E-6_wp ) THEN 2231 2081 a22(k) = SQRT( a22(k) ) 2232 2082 a32(k) = r32(k) - a21(k) * a31(k) / a22(k) 2233 ELSE 2083 ELSE 2234 2084 a22(k) = 10E-8_wp 2235 2085 a32(k) = 10E-8_wp … … 2243 2093 a33(k) = 10E-8_wp 2244 2094 ENDIF 2245 ENDDO 2095 ENDDO 2246 2096 ! 2247 2097 !-- Set bottom boundary condition … … 2249 2099 a22(nzb) = a22(nzb+1) 2250 2100 a21(nzb) = a21(nzb+1) 2251 a33(nzb) = a33(nzb+1) 2101 a33(nzb) = a33(nzb+1) 2252 2102 a31(nzb) = a31(nzb+1) 2253 a32(nzb) = a32(nzb+1) 2103 a32(nzb) = a32(nzb+1) 2254 2104 2255 2105 END SUBROUTINE calc_coeff_matrix 2256 2257 !------------------------------------------------------------------------------ !2106 2107 !--------------------------------------------------------------------------------------------------! 2258 2108 ! Description: 2259 2109 ! ------------ 2260 !> This routine controls the re-adjustment of the turbulence statistics used 2261 !> for generating turbulence at the lateral boundaries.2262 !------------------------------------------------------------------------------ !2110 !> This routine controls the re-adjustment of the turbulence statistics used for generating 2111 !> turbulence at the lateral boundaries. 2112 !--------------------------------------------------------------------------------------------------! 2263 2113 SUBROUTINE stg_adjust 2264 2114 2265 2115 IF ( debug_output_timestep ) CALL debug_message( 'stg_adjust', 'start' ) 2266 2116 ! 2267 !-- In case of dirichlet inflow boundary conditions only at one lateral 2268 !-- boundary, i.e. in the case no offline or self nesting is applied but 2269 !-- synthetic turbulence shall be parametrized nevertheless, the 2270 !-- boundary-layer depth need to determined first. 2271 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) & 2272 CALL nesting_offl_calc_zi 2273 ! 2274 !-- Compute scaling parameters (domain-averaged), such as friction velocity 2275 !-- are calculated. 2117 !-- In case of dirichlet inflow boundary conditions only at one lateral boundary, i.e. in the case 2118 !-- no offline or self nesting is applied but synthetic turbulence shall be parametrized 2119 !-- nevertheless, the boundary-layer depth need to determined first. 2120 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) & 2121 CALL nesting_offl_calc_zi 2122 ! 2123 !-- Compute scaling parameters (domain-averaged), such as friction velocity are calculated. 2276 2124 CALL calc_scaling_variables 2277 2125 ! … … 2279 2127 CALL calc_length_and_time_scale 2280 2128 ! 2281 !-- Parametrize Reynolds-stress tensor, diagonal elements as well 2282 !-- as r21 (v'u'), r31 (w'u'), r32 (w'v'). Parametrization follows 2283 !-- Rotach et al. (1996) and is based on boundary-layer depth, 2284 !-- friction velocity and velocity scale. 2129 !-- Parametrize Reynolds-stress tensor, diagonal elements as well as r21 (v'u'), r31 (w'u'), 2130 !-- r32 (w'v'). Parametrization follows Rotach et al. (1996) and is based on boundary-layer depth, 2131 !-- friction velocity and velocity scale. 2285 2132 CALL parametrize_reynolds_stress 2286 2133 ! 2287 !-- Calculate coefficient matrix from Reynolds stress tensor 2288 !-- (Lund rotation) 2134 !-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation) 2289 2135 CALL calc_coeff_matrix 2290 2136 ! … … 2304 2150 2305 2151 IF ( debug_output_timestep ) CALL debug_message( 'stg_adjust', 'end' ) 2306 2307 END SUBROUTINE stg_adjust 2308 2309 2310 !------------------------------------------------------------------------------ !2152 2153 END SUBROUTINE stg_adjust 2154 2155 2156 !--------------------------------------------------------------------------------------------------! 2311 2157 ! Description: 2312 2158 ! ------------ 2313 !> Calculates turbuluent length and time scales if these are not available 2314 !> from measurements. 2315 !------------------------------------------------------------------------------! 2159 !> Calculates turbuluent length and time scales if these are not available from measurements. 2160 !--------------------------------------------------------------------------------------------------! 2316 2161 SUBROUTINE calc_length_and_time_scale 2317 2162 2318 INTEGER(iwp) :: k !< loop index in z-direction2319 2320 REAL(wp) :: length_scale_dum 2321 2322 ! 2323 !-- In initial call the boundary-layer depth can be zero. This case, set 2324 !-- minimum value forboundary-layer depth, to setup length scales correctly.2163 INTEGER(iwp) :: k !< loop index in z-direction 2164 2165 REAL(wp) :: length_scale_dum !< effectively used length scale 2166 2167 ! 2168 !-- In initial call the boundary-layer depth can be zero. This case, set minimum value for 2169 !-- boundary-layer depth, to setup length scales correctly. 2325 2170 zi_ribulk = MAX( zi_ribulk, zw(nzb+2) ) 2326 2171 ! 2327 !-- Set-up default turbulent length scales. From the numerical point of 2328 !-- view the imposed perturbations should not be immediately dissipated 2329 !-- by the numerics. The numerical dissipation, however, acts on scales 2330 !-- up to 8 x the grid spacing. For this reason, set the turbulence 2331 !-- length scale to 8 time the grid spacing. Further, above the boundary 2332 !-- layer height, set turbulence lenght scales to zero (equivalent to not 2333 !-- imposing any perturbations) in order to save computational costs. 2334 !-- Typical time scales are derived by assuming Taylors's hypothesis, 2335 !-- using the length scales and the mean profiles of u- and v-component. 2172 !-- Set-up default turbulent length scales. From the numerical point of view the imposed 2173 !-- perturbations should not be immediately dissipated by the numerics. The numerical dissipation, 2174 !-- however, acts on scales up to 8 x the grid spacing. For this reason, set the turbulence length 2175 !-- scale to 8 time the grid spacing. Further, above the boundary layer height, set turbulence 2176 !-- lenght scales to zero (equivalent to not imposing any perturbations) in order to save 2177 !-- computational costs. Typical time scales are derived by assuming Taylors's hypothesis, using 2178 !-- the length scales and the mean profiles of u- and v-component. 2336 2179 DO k = nzb+1, nzt+1 2337 2180 ! 2338 !-- Determine blending parameter. Within the boundary layer length scales 2339 !-- are constant, while above lengths scales approach gradully zero, 2340 !-- i.e. imposed turbulence is not correlated in space and time, 2341 !-- just white noise, which saves computations power as the loops for the 2342 !-- computation of the filter functions depend on the length scales. 2343 !-- The value decreases to 1/10 after one length scale above the 2344 !-- ABL top. 2181 !-- Determine blending parameter. Within the boundary layer length scales are constant, while 2182 !-- above lengths scales approach gradully zero, i.e. imposed turbulence is not correlated in 2183 !-- space and time, just white noise, which saves computations power as the loops for the 2184 !-- computation of the filter functions depend on the length scales. The value decreases to 2185 !-- 1/10 after one length scale above the ABL top. 2345 2186 blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) ) 2346 2187 ! 2347 2188 !-- Assume isotropic turbulence length scales 2348 nux(k) = MAX( INT( length_scale * ddx ), 1 ) * blend 2349 nuy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend 2350 nvx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend 2351 nvy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend 2352 nwx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend 2353 nwy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend 2354 ! 2355 !-- Along the vertical direction limit the length scale further by the 2356 !-- boundary-layer depth to assure that no length scales larger than 2357 !-- the boundary-layer depth are used 2189 nux(k) = MAX( INT( length_scale * ddx ), 1 ) * blend 2190 nuy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend 2191 nvx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend 2192 nvy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend 2193 nwx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend 2194 nwy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend 2195 ! 2196 !-- Along the vertical direction limit the length scale further by the boundary-layer depth to 2197 !-- assure that no length scales larger than the boundary-layer depth are used 2358 2198 length_scale_dum = MIN( length_scale, zi_ribulk ) 2359 2199 2360 2200 nuz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend 2361 2201 nvz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend 2362 2202 nwz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend 2363 2203 ! 2364 !-- Limit time scales, else they become very larger for low wind speed, 2365 !-- imposing long-living inflow perturbations which in turn propagate 2366 !-- further into the model domain. Use u_init and v_init to calculate 2367 !-- the time scales, which will be equal to the inflow profiles, both, 2368 !-- in offline nesting mode or in dirichlet/radiation mode. 2369 tu(k) = MIN( dt_stg_adjust, length_scale / & 2370 ( ABS( u_init(k) ) + 0.1_wp ) ) * blend 2371 tv(k) = MIN( dt_stg_adjust, length_scale / & 2372 ( ABS( v_init(k) ) + 0.1_wp ) ) * blend 2373 ! 2374 !-- Time scale of w-component is a mixture from u- and v-component. 2204 !-- Limit time scales, else they become very larger for low wind speed, imposing long-living 2205 !-- inflow perturbations which in turn propagate further into the model domain. Use u_init and 2206 !-- v_init to calculate the time scales, which will be equal to the inflow profiles, both, in 2207 !-- offline nesting mode or in dirichlet/radiation mode. 2208 tu(k) = MIN( dt_stg_adjust, length_scale / ( ABS( u_init(k) ) + 0.1_wp ) ) * blend 2209 tv(k) = MIN( dt_stg_adjust, length_scale / ( ABS( v_init(k) ) + 0.1_wp ) ) * blend 2210 ! 2211 !-- Time scale of w-component is a mixture from u- and v-component. 2375 2212 tw(k) = SQRT( tu(k)**2 + tv(k)**2 ) * blend 2376 2213 2377 2214 ENDDO 2378 2215 ! 2379 !-- Set bottom boundary condition for the length and time scales 2216 !-- Set bottom boundary condition for the length and time scales 2380 2217 nux(nzb) = nux(nzb+1) 2381 2218 nuy(nzb) = nuy(nzb+1) … … 2393 2230 2394 2231 2395 END SUBROUTINE calc_length_and_time_scale 2396 2397 !------------------------------------------------------------------------------ !2232 END SUBROUTINE calc_length_and_time_scale 2233 2234 !--------------------------------------------------------------------------------------------------! 2398 2235 ! Description: 2399 2236 ! ------------ 2400 !> Calculate scaling variables which are used for turbulence parametrization 2401 !> according to Rotach et al. (1996). Scaling variables are: friction velocity,2402 !> boundary-layer depth, momentum velocity scale, and Obukhov length.2403 !------------------------------------------------------------------------------ !2237 !> Calculate scaling variables which are used for turbulence parametrization according to Rotach 2238 !> et al. (1996). Scaling variables are: friction velocity, boundary-layer depth, momentum velocity 2239 !> scale, and Obukhov length. 2240 !--------------------------------------------------------------------------------------------------! 2404 2241 SUBROUTINE calc_scaling_variables 2405 2242 2406 INTEGER(iwp) :: i!< loop index in x-direction2407 INTEGER(iwp) :: j!< loop index in y-direction2408 INTEGER(iwp) :: k!< loop index in z-direction2409 INTEGER(iwp) :: m!< surface element index2410 2411 REAL(wp) :: friction_vel_l 2412 REAL(wp) :: pt_surf_mean 2413 REAL(wp) :: pt_surf_mean_l 2414 REAL(wp) :: scale_l_l 2415 REAL(wp) :: shf_mean 2416 REAL(wp) :: shf_mean_l 2417 REAL(wp) :: w_convective 2418 2419 ! 2420 !-- Calculate mean friction velocity, velocity scale, heat flux and 2421 !-- near-surface temperature in the model domain.2243 INTEGER(iwp) :: i !< loop index in x-direction 2244 INTEGER(iwp) :: j !< loop index in y-direction 2245 INTEGER(iwp) :: k !< loop index in z-direction 2246 INTEGER(iwp) :: m !< surface element index 2247 2248 REAL(wp) :: friction_vel_l !< mean friction veloctiy on subdomain 2249 REAL(wp) :: pt_surf_mean !< mean near surface temperature (at 1st grid point) 2250 REAL(wp) :: pt_surf_mean_l !< mean near surface temperature (at 1st grid point) on subdomain 2251 REAL(wp) :: scale_l_l !< mean Obukhov lenght on subdomain 2252 REAL(wp) :: shf_mean !< mean surface sensible heat flux 2253 REAL(wp) :: shf_mean_l !< mean surface sensible heat flux on subdomain 2254 REAL(wp) :: w_convective !< convective velocity scale 2255 2256 ! 2257 !-- Calculate mean friction velocity, velocity scale, heat flux and near-surface temperature in the 2258 !-- model domain. 2422 2259 pt_surf_mean_l = 0.0_wp 2423 2260 shf_mean_l = 0.0_wp … … 2432 2269 scale_l_l = scale_l_l + surf_def_h(0)%ol(m) 2433 2270 pt_surf_mean_l = pt_surf_mean_l + pt(k,j,i) 2434 ENDDO 2271 ENDDO 2435 2272 DO m = 1, surf_lsm_h%ns 2436 2273 i = surf_lsm_h%i(m) … … 2451 2288 pt_surf_mean_l = pt_surf_mean_l + pt(k,j,i) 2452 2289 ENDDO 2453 2290 2454 2291 #if defined( __parallel ) 2455 CALL MPI_ALLREDUCE( friction_vel_l, scale_us, 1, MPI_REAL, MPI_SUM, & 2456 comm2d, ierr ) 2457 CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM, & 2458 comm2d, ierr ) 2459 CALL MPI_ALLREDUCE( scale_l_l, scale_l, 1, MPI_REAL, MPI_SUM, & 2460 comm2d, ierr ) 2461 CALL MPI_ALLREDUCE( pt_surf_mean_l, pt_surf_mean, 1, MPI_REAL, MPI_SUM, & 2462 comm2d, ierr ) 2292 CALL MPI_ALLREDUCE( friction_vel_l, scale_us, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2293 CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2294 CALL MPI_ALLREDUCE( scale_l_l, scale_l, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2295 CALL MPI_ALLREDUCE( pt_surf_mean_l, pt_surf_mean, 1, MPI_REAL, MPI_SUM, comm2d, ierr ) 2463 2296 #else 2464 2297 scale_us = friction_vel_l … … 2471 2304 shf_mean = shf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2472 2305 scale_l = scale_l / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2473 pt_surf_mean = pt_surf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2474 ! 2475 !-- Compute mean convective velocity scale. Note, in case the mean heat flux 2476 !-- is negative, setconvective velocity scale to zero.2306 pt_surf_mean = pt_surf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2307 ! 2308 !-- Compute mean convective velocity scale. Note, in case the mean heat flux is negative, set 2309 !-- convective velocity scale to zero. 2477 2310 IF ( shf_mean > 0.0_wp ) THEN 2478 2311 w_convective = ( g * shf_mean * zi_ribulk / pt_surf_mean )**( 1.0_wp / 3.0_wp ) … … 2481 2314 ENDIF 2482 2315 ! 2483 !-- Finally, in order to consider also neutral or stable stratification, 2484 !-- compute momentum velocity scale from u* and convective velocity scale, 2485 !-- according to Rotach et al. (1996). 2316 !-- Finally, in order to consider also neutral or stable stratification, compute momentum velocity 2317 !-- scale from u* and convective velocity scale, according to Rotach et al. (1996). 2486 2318 scale_wm = ( scale_us**3 + 0.6_wp * w_convective**3 )**( 1.0_wp / 3.0_wp ) 2487 2319 2488 2320 END SUBROUTINE calc_scaling_variables 2489 2321
Note: See TracChangeset
for help on using the changeset viewer.
