source: palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

Last change on this file was 4851, checked in by gronemeier, 3 years ago

bugfix: deactivated header output (dynamics_mod); change: formatting clean-up (synthetic_turbulence_generator_mod)

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