source: palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

Last change on this file was 4851, checked in by gronemeier, 6 months 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
Line 
1!> @synthetic_turbulence_generator_mod.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: synthetic_turbulence_generator_mod.f90 4851 2021-01-22 09:25:05Z banzhafs $
27! formatting clean-up
28!
29! 4848 2021-01-21 15:51:51Z gronemeier
30! replaced use_syn_turb_gen by syn_turb_gen
31!
32! 4843 2021-01-15 15:22:11Z raasch
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
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
41! Implementation of downward facing USM and LSM surfaces
42!
43! 4647 2020-08-24 16:36:18Z suehring
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
46! large length scales this is significantly faster.
47!
48! 4640 2020-08-11 16:28:32Z suehring
49! - to avoid that the correction term in r11/r22 computation becomes unrealistically high, limit
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
55! support for MPI Fortran77 interface (mpif.h) removed
56!
57! 4603 2020-07-14 16:08:30Z suehring
58! Bugfix in initialization from ASCII file - x-length scales at the bottom boundary were not
59! initialized properly
60!
61! 4566 2020-06-16 10:11:51Z suehring
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)
66!
67! 4562 2020-06-12 08:38:47Z raasch
68! Parts of r4559 re-formatted
69!
70! 4559 2020-06-11 08:51:48Z raasch
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!
76! 4495 2020-04-13 20:11:20Z raasch
77! Restart data handling with MPI-IO added
78!
79! 4481 2020-03-31 18:55:54Z maronga
80! Bugfix: cpp-directives for serial mode added, dummy statements to prevent compile errors added
81!
82! 4442 2020-03-04 19:21:13Z suehring
83! Set back turbulent length scale to 8 x grid spacing in the parametrized mode
84! (was accidantly changed).
85!
86! 4441 2020-03-04 19:20:35Z suehring
87! Correct misplaced preprocessor directive
88!
89! 4438 2020-03-03 20:49:28Z suehring
90! Performance optimizations in velocity-seed calculation:
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!
98! 4346 2019-12-18 11:55:56Z motisi
99! Introduction of wall_flags_total_0, which currently sets bits based on static topography
100! information used in wall_flags_static_0
101!
102! 4335 2019-12-12 16:39:05Z suehring
103! Commentation of last commit
104!
105! 4332 2019-12-10 19:44:12Z suehring
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!
109! 4329 2019-12-10 15:46:36Z motisi
110! Renamed wall_flags_0 to wall_flags_static_0
111!
112! 4309 2019-11-26 18:49:59Z suehring
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!
120! 4182 2019-08-22 15:20:23Z scharf
121! Corrected "Former revisions" section
122!
123! 4148 2019-08-08 11:26:00Z suehring
124! Remove unused variable
125!
126! 4144 2019-08-06 09:11:47Z raasch
127! Relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
128!
129! 4071 2019-07-03 20:02:00Z suehring
130! Bugfix, initialize mean_inflow_profiles in case turbulence and inflow information is not read from
131! file.
132!
133! 4022 2019-06-12 11:52:39Z suehring
134! Several bugfixes and improvements
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
143! - Bugfix in non-nested case: use ABL height for parametrized turbulence
144!
145! 3987 2019-05-22 09:52:13Z kanani
146! Introduce alternative switch for debug output during timestepping
147!
148! 3938 2019-04-29 16:06:25Z suehring
149! Remove unused variables
150!
151! 3937 2019-04-29 15:09:07Z suehring
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!
156! 3909 2019-04-17 09:13:25Z suehring
157! Minor bugfix for last commit
158!
159! 3900 2019-04-16 15:17:43Z suehring
160! Missing re-calculation of perturbation seeds in case of restarts
161!
162! 3891 2019-04-12 17:52:01Z suehring
163! Bugfix in initialization in case of restart runs.
164!
165! 3885 2019-04-11 11:29:34Z kanani
166! Changes related to global restructuring of location messages and introduction of additional debug
167! messages
168!
169!
170! Removed unused variables
171!
172! 3719 2019-02-06 13:10:18Z kanani
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!
177! 2259 2017-06-08 09:09:11Z gronemeier
178! Initial revision
179!
180! Authors:
181! --------
182! @author Tobias Gronemeier, Matthias Suehring, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
183!
184!
185! Description:
186! ------------
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.
193!>
194!> @todo Enable cyclic_fill
195!>       Implement turbulence generation for e and pt
196!> @note <Enter notes on the module>
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!--------------------------------------------------------------------------------------------------!
201 MODULE synthetic_turbulence_generator_mod
202
203
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,                                                                                 &
219               zw
220
221    USE basic_constants_and_equations_mod,                                                         &
222        ONLY:  g,                                                                                  &
223               kappa,                                                                              &
224               pi
225
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,                                                         &
249               turbulent_inflow
250
251    USE cpulog,                                                                                    &
252        ONLY:  cpu_log,                                                                            &
253               log_point_s
254
255    USE grid_variables,                                                                            &
256        ONLY:  ddx,                                                                                &
257               ddy,                                                                                &
258               dx,                                                                                 &
259               dy
260
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,                                                                                &
274               wall_flags_total_0
275
276    USE kinds
277
278#if defined( __parallel )
279    USE MPI
280#endif
281
282    USE nesting_offl_mod,                                                                          &
283        ONLY:  nesting_offl_calc_zi,                                                               &
284               zi_ribulk
285
286    USE pegrid,                                                                                    &
287        ONLY:  comm1dx,                                                                            &
288               comm1dy,                                                                            &
289               comm2d,                                                                             &
290               ierr,                                                                               &
291               myidx,                                                                              &
292               myidy,                                                                              &
293               pdims
294
295    USE pmc_interface,                                                                             &
296        ONLY : rans_mode_parent
297
298    USE random_generator_parallel,                                                                 &
299        ONLY:  init_parallel_random_generator,                                                     &
300               random_dummy,                                                                       &
301               random_number_parallel,                                                             &
302               random_seed_parallel
303
304    USE restart_data_mpi_io_mod,                                                                   &
305        ONLY:  rrd_mpi_io,                                                                         &
306               wrd_mpi_io
307
308    USE transpose_indices,                                                                         &
309        ONLY:  nzb_x,                                                                              &
310               nzt_x
311
312    USE surface_mod,                                                                               &
313        ONLY:  surf_def_h,                                                                         &
314               surf_lsm_h,                                                                         &
315               surf_usm_h
316
317    IMPLICIT NONE
318
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
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)
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)
333#if defined( __parallel )
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
338#endif
339
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
345#if defined( __parallel )
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
350#endif
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)
360
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
365
366
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
369    LOGICAL ::  compute_velocity_seeds_local  = .FALSE.   !< switch to decide whether velocity seeds are computed locally
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
374
375
376    REAL(wp) ::  blend                     !< value to create gradually and smooth blending of Reynolds stress and length
377                                           !< scales above the boundary layer
378    REAL(wp) ::  blend_coeff = -9.3_wp     !< coefficient used to ensure that blending functions decreases to 1/10 after
379                                           !< one length scale above ABL top
380    REAL(wp) ::  d_l                       !< blend_coeff/length_scale
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
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
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
434
435!
436!-- PALM interfaces:
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!
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!
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
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
490    INTERFACE stg_rrd_global
491       MODULE PROCEDURE stg_rrd_global_ftn
492       MODULE PROCEDURE stg_rrd_global_mpi
493    END INTERFACE stg_rrd_global
494
495!
496!-- Writing of binary output for restart runs
497    INTERFACE stg_wrd_global
498       MODULE PROCEDURE stg_wrd_global
499    END INTERFACE stg_wrd_global
500
501    SAVE
502
503    PRIVATE
504
505!
506!-- Public interfaces
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
515
516!
517!-- Public variables
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,                                                                       &
526            time_stg_call
527
528
529 CONTAINS
530
531
532!--------------------------------------------------------------------------------------------------!
533! Description:
534! ------------
535!> Check parameters routine for synthetic turbulence generator
536!--------------------------------------------------------------------------------------------------!
537 SUBROUTINE stg_check_parameters
538
539    IF ( .NOT. syn_turb_gen  .AND.  .NOT. rans_mode  .AND.                                         &
540          nesting_offline )  THEN
541       message_string = 'Synthetic turbulence generator is required ' //                           &
542                        'if offline nesting is applied and PALM operates in LES mode.'
543       CALL message( 'stg_check_parameters', 'PA0520', 0, 0, 0, 6, 0 )
544    ENDIF
545
546    IF ( .NOT. syn_turb_gen  .AND.  child_domain                                                   &
547         .AND. rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
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.'
550       CALL message( 'stg_check_parameters', 'PA0524', 1, 2, 0, 6, 0 )
551    ENDIF
552
553    IF ( syn_turb_gen )  THEN
554
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.'
557          CALL message( 'stg_check_parameters', 'PA0620', 1, 2, 0, 6, 0 )
558
559       ENDIF
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.'
563          CALL message( 'stg_check_parameters', 'PA0621', 1, 2, 0, 6, 0 )
564
565       ENDIF
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.'
574             CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
575          ENDIF
576
577          IF ( bc_lr /= 'dirichlet/radiation' )  THEN
578             message_string = 'Using synthetic turbulence generator requires &bc_lr = ' //         &
579                              ' "dirichlet/radiation", if not offline nesting is applied.'
580             CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
581          ENDIF
582          IF ( bc_ns /= 'cyclic' )  THEN
583             message_string = 'Using synthetic turbulence generator requires &bc_ns = ' //         &
584                              ' "cyclic", if not offline nesting is applied.'
585             CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
586          ENDIF
587
588       ENDIF
589
590       IF ( turbulent_inflow )  THEN
591          message_string = 'Using synthetic turbulence generator in combination ' //               &
592                           '&with turbulent_inflow = .T. is not allowed'
593          CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
594       ENDIF
595!
596!--    Synthetic turbulence generator requires the parallel random generator
597       IF ( random_generator /= 'random-parallel' )  THEN
598          message_string = 'Using synthetic turbulence generator requires random_generator = ' //  &
599                           'random-parallel.'
600          CALL message( 'stg_check_parameters', 'PA0421', 1, 2, 0, 6, 0 )
601       ENDIF
602
603    ENDIF
604
605 END SUBROUTINE stg_check_parameters
606
607
608!--------------------------------------------------------------------------------------------------!
609! Description:
610! ------------
611!> Header output for synthetic turbulence generator
612!--------------------------------------------------------------------------------------------------!
613 SUBROUTINE stg_header ( io )
614
615    INTEGER(iwp), INTENT(IN) ::  io  !< Unit of the output file
616
617!
618!-- Write synthetic turbulence generator Header
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
626    ENDIF
627
6281   FORMAT (//' Synthetic turbulence generator information:'/                                      &
629              ' ------------------------------------------'/)
6304   FORMAT ('    imposed turbulence statistics are parametrized and ajdusted to boundary-layer development each ', F8.2, ' s' )
6315   FORMAT ('    imposed turbulence is read from file' )
632
633 END SUBROUTINE stg_header
634
635
636!--------------------------------------------------------------------------------------------------!
637! Description:
638! ------------
639!> Initialization of the synthetic turbulence generator
640!--------------------------------------------------------------------------------------------------!
641 SUBROUTINE stg_init
642
643#if defined( __parallel )
644    INTEGER(KIND=MPI_ADDRESS_KIND) :: extent   !< extent of new MPI type
645    INTEGER(KIND=MPI_ADDRESS_KIND) :: tob = 0  !< dummy variable
646#endif
647
648    INTEGER(iwp) :: i         !> grid index in x-direction
649    INTEGER(iwp) :: j         !> loop index
650    INTEGER(iwp) :: k         !< index
651#if defined( __parallel )
652    INTEGER(iwp) :: newtype   !< dummy MPI type
653    INTEGER(iwp) :: realsize  !< size of REAL variables
654#endif
655
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
658
659
660    LOGICAL ::  file_stg_exist = .FALSE.  !< flag indicating whether parameter file for Reynolds stress and length scales exist
661
662!
663!-- Dummy variables used for reading profiles
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
674#if defined( __parallel )
675    REAL(wp) :: nnz  !< increment used to determine processor decomposition of z-axis along x and y direction
676#endif
677    REAL(wp) :: zz   !< height
678
679
680#if defined( __parallel )
681    CALL MPI_BARRIER( comm2d, ierr )
682#endif
683!
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).
689    IF ( .NOT. compute_velocity_seeds_local )  THEN
690#if defined( __parallel )
691!
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 )
696
697       IF ( MOD( nz , pdims(1) ) /= 0  .AND.  myidx == id_stg_right )                              &
698          nzt_x_stg = nzt_x_stg + myidx * ( nnz - INT( nnz ) )
699
700       IF ( nesting_offline   .OR.  ( child_domain  .AND.  rans_mode_parent                        &
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 )
705
706          IF ( MOD( nz , pdims(2) ) /= 0  .AND.  myidy == id_stg_north )                           &
707             nzt_y_stg = nzt_y_stg + myidy * ( nnz - INT( nnz ) )
708       ENDIF
709
710!
711!--    Define MPI type used in stg_generate_seed_yz to gather vertical splitted velocity seeds
712       CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
713       extent = 1 * realsize
714!
715!--    Set-up MPI datatyp to involve all cores for turbulence generation at yz layer
716!--    stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
717       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyn-nys+1],                                    &
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 )
721       CALL MPI_TYPE_FREE( newtype, ierr )
722
723       ! stg_type_yz_small: yz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
724       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyn-nys+1],                        &
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 )
728       CALL MPI_TYPE_FREE( newtype, ierr )
729
730       ! Receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
731       ALLOCATE( recv_count_yz(pdims(1)), displs_yz(pdims(1)) )
732
733       recv_count_yz           = nzt_x_stg-nzb_x_stg + 1
734       recv_count_yz(pdims(1)) = recv_count_yz(pdims(1)) + 1
735
736       DO  j = 1, pdims(1)
737          displs_yz(j) = 0 + (nzt_x_stg-nzb_x_stg+1) * (j-1)
738       ENDDO
739!
740!--    Set-up MPI datatyp to involve all cores for turbulence generation at xz layer
741!--    stg_type_xz: xz-slice with vertical bounds nzb:nzt+1
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],                                 &
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 )
749
750          ! stg_type_yz_small: xz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
751          CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_y_stg-nzb_y_stg+2,nxr-nxl+1],                     &
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 )
756
757          ! Receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
758          ALLOCATE( recv_count_xz(pdims(2)), displs_xz(pdims(2)) )
759
760          recv_count_xz           = nzt_y_stg-nzb_y_stg + 1
761          recv_count_xz(pdims(2)) = recv_count_xz(pdims(2)) + 1
762
763          DO  j = 1, pdims(2)
764             displs_xz(j) = 0 + (nzt_y_stg-nzb_y_stg+1) * (j-1)
765          ENDDO
766
767       ENDIF
768#endif
769    ENDIF
770!
771!-- Allocate required arrays.
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
777       IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )  THEN
778          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
779          mean_inflow_profiles = 0.0_wp
780          mean_inflow_profiles(:,1) = u_init
781          mean_inflow_profiles(:,2) = v_init
782!
783!--       Even though potential temperature and humidity are not perturbed, they need to be
784!--       initialized appropriately.
785          IF ( .NOT. neutral )                                                                     &
786             mean_inflow_profiles(:,4) = pt_init
787          IF ( humidity )                                                                          &
788             mean_inflow_profiles(:,6) = q_init
789       ENDIF
790    ENDIF
791!
792!-- Assign initial profiles. Note, this is only required if turbulent inflow from the left is
793!-- desired, not in case of any of the nesting (offline or self nesting) approaches.
794    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
795       u_init = mean_inflow_profiles(:,1)
796       v_init = mean_inflow_profiles(:,2)
797       e_init = MAXVAL( mean_inflow_profiles(:,5) )
798    ENDIF
799
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
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
819    ALLOCATE ( dist_xz(nzb:nzt+1,nxl:nxr,3) )
820    ALLOCATE ( dist_yz(nzb:nzt+1,nys:nyn,3) )
821    dist_xz = 0.0_wp
822    dist_yz = 0.0_wp
823!
824!-- Read inflow profile
825!-- Height levels of profiles in input profile is as follows:
826!-- zu: luy, luz, tu, lvy, lvz, tv, r11, r21, r22, d1, d2, d5 zw: lwy, lwz, tw, r31, r32, r33, d3
827!-- WARNING: zz is not used at the moment
828    INQUIRE( FILE = 'STG_PROFILES' // TRIM( coupling_char ), EXIST = file_stg_exist  )
829
830    IF ( file_stg_exist )  THEN
831
832       OPEN( 90, FILE = 'STG_PROFILES' // TRIM( coupling_char ), STATUS = 'OLD', FORM = 'FORMATTED' )
833!
834!--    Skip header
835       READ( 90, * )
836
837       DO  k = nzb, nzt
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
840
841!
842!--       Convert length scales from meter to number of grid points.
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
853          nuy(k) = INT( luy * ddy )
854          nvy(k) = INT( lvy * ddy )
855          nwy(k) = INT( lwy * ddy )
856!
857!--       Set length scales in x-direction. As a workaround assume isotropic turbulence in x- and
858!--       y-direction.
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
871!
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)
877       nuy(nzb) = nuy(nzb+1)
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)
883
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
894       CLOSE( 90 )
895!
896!--    Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
897       CALL calc_coeff_matrix
898!
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.
904    ELSE
905!
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 )
908!
909!--    Set flag indicating that turbulence is parametrized
910       parametrize_inflow_turbulence = .TRUE.
911!
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 )                                    &
916          CALL nesting_offl_calc_zi
917!
918!--    Determine boundary-layer depth, which is used to initialize lenght scales
919       CALL calc_scaling_variables
920!
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.
924       CALL parametrize_turbulence
925!
926!--    Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
927       CALL calc_coeff_matrix
928
929    ENDIF
930!
931!-- Initial filter functions
932    CALL stg_setup_filter_function
933!
934!-- Allocate velocity seeds for turbulence at xz-layer
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) )
938
939!
940!-- Allocate velocity seeds for turbulence at yz-layer
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) )
944
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
959#if defined( __parallel )
960    CALL MPI_BARRIER( comm2d, ierr )
961#endif
962
963!
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.
967    IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
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
973          DO  j = nys, nyn
974             DO  k = nzb, nzt+1
975                IF  ( a11(k) > 10E-8_wp )  THEN
976                   fu_yz(k,j) = ( u(k,j,i) - u_init(k) ) / a11(k)
977                ELSE
978                   fu_yz(k,j) = 10E-8_wp
979                ENDIF
980
981                IF  ( a22(k) > 10E-8_wp )  THEN
982                   fv_yz(k,j) = ( v(k,j,i) - a21(k) * fu_yz(k,j) - v_init(k) ) / a22(k)
983                ELSE
984                   fv_yz(k,j) = 10E-8_wp
985                ENDIF
986
987                IF  ( a33(k) > 10E-8_wp )  THEN
988                   fw_yz(k,j) = ( w(k,j,i) - a31(k) * fu_yz(k,j) - a32(k) * fv_yz(k,j) ) / a33(k)
989                ELSE
990                   fw_yz(k,j) = 10E-8_wp
991                ENDIF
992             ENDDO
993          ENDDO
994       ENDIF
995
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
1001          DO  i = nxl, nxr
1002             DO  k = nzb, nzt+1
1003!
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.
1007                IF  ( a11(k) > 10E-8_wp )  THEN
1008                   fu_xz(k,i) = ( u(k,j,i) - u_init(k) ) / a11(k)
1009                ELSE
1010                   fu_xz(k,i) = 10E-8_wp
1011                ENDIF
1012
1013                IF  ( a22(k) > 10E-8_wp )  THEN
1014                   fv_xz(k,i) = ( v(k,j,i) - a21(k) * fu_xz(k,i) - v_init(k) ) / a22(k)
1015                ELSE
1016                   fv_xz(k,i) = 10E-8_wp
1017                ENDIF
1018
1019                IF  ( a33(k) > 10E-8_wp )  THEN
1020                   fw_xz(k,i) = ( w(k,j,i) - a31(k) * fu_xz(k,i) - a32(k) * fv_xz(k,i) ) / a33(k)
1021                ELSE
1022                   fw_xz(k,i) = 10E-8_wp
1023                ENDIF
1024
1025             ENDDO
1026          ENDDO
1027       ENDIF
1028    ENDIF
1029!
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.
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
1039
1040       nr_non_topo_yz_l(1) = SUM( MERGE( 1, 0, BTEST( wall_flags_total_0(nzb:nzt,nys:nyn,i), 1 ) ) )
1041!
1042!--    Number of grid points where perturbations are imposed on v and w
1043       IF ( myidx == id_stg_left  )  i = nxl-1
1044       IF ( myidx == id_stg_right )  i = nxr+1
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
1049#if defined( __parallel )
1050       CALL MPI_ALLREDUCE( nr_non_topo_yz_l, nr_non_topo_yz, 3, MPI_INTEGER, MPI_SUM, comm1dy, ierr )
1051#else
1052       nr_non_topo_yz = nr_non_topo_yz_l
1053#endif
1054    ENDIF
1055
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
1061
1062       nr_non_topo_xz_l(2) = SUM( MERGE( 1, 0, BTEST( wall_flags_total_0(nzb:nzt,j,nxl:nxr), 2 ) ) )
1063!
1064!--    Number of grid points where perturbations are imposed on u and w
1065       IF ( myidy == id_stg_south )  j = nys-1
1066       IF ( myidy == id_stg_north )  j = nyn+1
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
1071#if defined( __parallel )
1072       CALL MPI_ALLREDUCE( nr_non_topo_xz_l, nr_non_topo_xz, 3, MPI_INTEGER, MPI_SUM, comm1dx, ierr )
1073#else
1074       nr_non_topo_xz = nr_non_topo_xz_l
1075#endif
1076    ENDIF
1077!
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.
1081    ALLOCATE ( id_rand_yz(-mergp_limit+nys:nyn+mergp_limit) )
1082    ALLOCATE ( seq_rand_yz(5,-mergp_limit+nys:nyn+mergp_limit) )
1083    id_rand_yz  = 0
1084    seq_rand_yz = 0
1085
1086    CALL init_parallel_random_generator( ny, -mergp_limit+nys, nyn+mergp_limit, id_rand_yz, seq_rand_yz )
1087
1088    IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent                             &
1089         .AND.  .NOT.  rans_mode ) )  THEN
1090       ALLOCATE ( id_rand_xz(-mergp_limit+nxl:nxr+mergp_limit) )
1091       ALLOCATE ( seq_rand_xz(5,-mergp_limit+nxl:nxr+mergp_limit) )
1092       id_rand_xz  = 0
1093       seq_rand_xz = 0
1094
1095       CALL init_parallel_random_generator( nx, -mergp_limit+nxl, nxr+mergp_limit, id_rand_xz, seq_rand_xz )
1096    ENDIF
1097
1098 END SUBROUTINE stg_init
1099
1100
1101!--------------------------------------------------------------------------------------------------!
1102! Description:
1103! ------------
1104!> Calculate filter function bxx from length scale nxx following Eg.9 and 10 (Xie and Castro, 2008)
1105!--------------------------------------------------------------------------------------------------!
1106 SUBROUTINE stg_filter_func( nxx, bxx, mergp )
1107
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
1112
1113    INTEGER(iwp), DIMENSION(nzb:nzt+1) ::  nxx  !< length scale (in gp)
1114
1115    REAL(wp) :: bdenom        !< denominator for filter functions bXX
1116    REAL(wp) :: qsi = 1.0_wp  !< minimization factor
1117
1118    REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) ::  bxx  !< filter function
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
1130          DO  nf = -n_k, n_k
1131             bdenom = bdenom + EXP( -qsi * pi * ABS( nf ) / n_k )**2
1132          ENDDO
1133
1134!
1135!--       ( Eq.9 )
1136          bdenom = SQRT( bdenom )
1137          DO  nf = -n_k, n_k
1138             bxx(nf,k) = EXP( -qsi * pi * ABS( nf ) / n_k ) / bdenom
1139          ENDDO
1140       ENDIF
1141    ENDDO
1142
1143 END SUBROUTINE stg_filter_func
1144
1145
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
1210!--------------------------------------------------------------------------------------------------!
1211! Description:
1212! ------------
1213!> Parin for &stg_par for synthetic turbulence generator
1214!--------------------------------------------------------------------------------------------------!
1215 SUBROUTINE stg_parin
1216
1217    CHARACTER(LEN=100) ::  line  !< dummy string that contains the current line of the parameter file
1218
1219    INTEGER(iwp) ::  io_status   !< status after reading the namelist file
1220
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,                                                             &
1227                        dt_stg_call,                                                               &
1228                        switch_off_module
1229
1230
1231!
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 )
1235
1236!
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.
1242       IF ( .NOT. switch_off_module )  syn_turb_gen = .TRUE.
1243
1244    ELSEIF ( io_status > 0 )  THEN
1245!
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 )
1251
1252    ENDIF
1253
1254 END SUBROUTINE stg_parin
1255
1256
1257!--------------------------------------------------------------------------------------------------!
1258! Description:
1259! ------------
1260!> Read module-specific global restart data (Fortran binary format).
1261!--------------------------------------------------------------------------------------------------!
1262 SUBROUTINE stg_rrd_global_ftn( found )
1263
1264    LOGICAL, INTENT(OUT) ::  found !< flag indicating if variable was found
1265
1266    found = .TRUE.
1267
1268
1269    SELECT CASE ( restart_string(1:length) )
1270
1271       CASE ( 'time_stg_adjust' )
1272          READ ( 13 )  time_stg_adjust
1273
1274       CASE ( 'time_stg_call' )
1275          READ ( 13 )  time_stg_call
1276
1277       CASE DEFAULT
1278
1279          found = .FALSE.
1280
1281    END SELECT
1282
1283
1284 END SUBROUTINE stg_rrd_global_ftn
1285
1286
1287!--------------------------------------------------------------------------------------------------!
1288! Description:
1289! ------------
1290!> Read module-specific global restart data (MPI-IO).
1291!--------------------------------------------------------------------------------------------------!
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
1300!--------------------------------------------------------------------------------------------------!
1301! Description:
1302! ------------
1303!> This routine writes the respective restart data.
1304!--------------------------------------------------------------------------------------------------!
1305 SUBROUTINE stg_wrd_global
1306
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
1311
1312       CALL wrd_write_string( 'time_stg_call' )
1313       WRITE ( 14 )  time_stg_call
1314
1315    ELSEIF ( restart_data_format_output(1:3) == 'mpi' )  THEN
1316
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
1322 END SUBROUTINE stg_wrd_global
1323
1324
1325!--------------------------------------------------------------------------------------------------!
1326! Description:
1327! ------------
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!--------------------------------------------------------------------------------------------------!
1332 SUBROUTINE stg_main
1333
1334    USE exchange_horiz_mod,                                                                        &
1335        ONLY:  exchange_horiz
1336
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
1340
1341    LOGICAL  ::  stg_call = .FALSE.  !< control flag indicating whether turbulence was updated or only restored from last call
1342
1343    REAL(wp) ::  dt_stg  !< time interval the STG is called
1344
1345    REAL(wp), DIMENSION(3) ::  mc_factor_l  !< local mass flux correction factor
1346
1347    IF ( debug_output_timestep )  CALL debug_message( 'stg_main', 'start' )
1348!
1349!-- Calculate time step which is needed for filter functions
1350    dt_stg = MAX( dt_3d, dt_stg_call )
1351!
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.                                                       &
1355         intermediate_timestep_count == intermediate_timestep_count_max  )  THEN
1356       stg_call = .TRUE.
1357    ELSE
1358       stg_call = .FALSE.
1359    ENDIF
1360!
1361!-- Initial value of fu, fv, fw
1362    IF ( time_since_reference_point == 0.0_wp .AND. .NOT. velocity_seed_initialized )  THEN
1363!
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
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
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 )
1380!
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 )
1387       ENDIF
1388       velocity_seed_initialized = .TRUE.
1389    ENDIF
1390!
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.
1394    IF ( stg_call )  THEN
1395       IF ( .NOT. ( nesting_offline  .OR.                                                          &
1396                  ( child_domain .AND.  rans_mode_parent                                           &
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
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 )
1406
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 )
1410       ENDIF
1411    ENDIF
1412
1413!
1414!-- Turbulence generation at left and/or right boundary
1415    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
1416!
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
1419       IF ( stg_call )  THEN
1420
1421          DO  j = nys, nyn
1422             DO  k = nzb, nzt + 1
1423!
1424!--             Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1425                IF ( tu(k) == 0.0_wp  .OR.  adjustment_step )  THEN
1426                   fu_yz(k,j) = fuo_yz(k,j)
1427                ELSE
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) ) )
1430                ENDIF
1431
1432                IF ( tv(k) == 0.0_wp  .OR.  adjustment_step )  THEN
1433                   fv_yz(k,j) = fvo_yz(k,j)
1434                ELSE
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) ) )
1437                ENDIF
1438
1439                IF ( tw(k) == 0.0_wp  .OR.  adjustment_step )  THEN
1440                   fw_yz(k,j) = fwo_yz(k,j)
1441                ELSE
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) ) )
1444                ENDIF
1445             ENDDO
1446          ENDDO
1447
1448          dist_yz(nzb,:,1) = 0.0_wp
1449          dist_yz(nzb,:,2) = 0.0_wp
1450          dist_yz(nzb,:,3) = 0.0_wp
1451
1452          IF ( myidx == id_stg_left  )  i = nxl
1453          IF ( myidx == id_stg_right )  i = nxr+1
1454          DO  j = nys, nyn
1455             DO  k = nzb+1, nzt + 1
1456!
1457!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
1458                dist_yz(k,j,1) = a11(k) * fu_yz(k,j) *                                             &
1459                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
1460
1461             ENDDO
1462          ENDDO
1463
1464          IF ( myidx == id_stg_left  )  i = nxl-1
1465          IF ( myidx == id_stg_right )  i = nxr+1
1466          DO  j = nys, nyn
1467             DO  k = nzb+1, nzt + 1
1468!
1469!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
1470!--             Additional factors are added to improve the variance of v and w experimental test
1471!--             of 1.2
1472                dist_yz(k,j,2) = ( a21(k) * fu_yz(k,j) + a22(k) * fv_yz(k,j) ) *                   &
1473                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
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
1478             ENDDO
1479          ENDDO
1480       ENDIF
1481!
1482!--    Mass flux correction following Kim et al. (2013)
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.
1486       IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
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.
1492          IF ( myidx == id_stg_left  )  i = nxl
1493          IF ( myidx == id_stg_right )  i = nxr+1
1494
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 ) ) )
1497
1498          IF ( myidx == id_stg_left  )  i = nxl-1
1499          IF ( myidx == id_stg_right )  i = nxr+1
1500
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 ) ) )
1505
1506#if defined( __parallel )
1507          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 3, MPI_REAL, MPI_SUM, comm1dy, ierr )
1508#else
1509          mc_factor   = mc_factor_l
1510#endif
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
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 ) )
1520
1521
1522          IF ( myidx == id_stg_left  )  i = nxl-1
1523          IF ( myidx == id_stg_right )  i = nxr+1
1524
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 ) )
1527
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 ) )
1530!
1531!--       Add disturbances
1532          IF ( myidx == id_stg_left  )  THEN
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!
1539!--             Add disturbance at the inflow
1540                DO  j = nys, nyn
1541                   DO  k = nzb, nzt+1
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 ) )
1548                   ENDDO
1549                ENDDO
1550             ELSE
1551
1552                DO  j = nys, nyn
1553                   DO  k = nzb+1, nzt
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 ) )
1556                      u(k,j,-1)  = u(k,j,0)
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 ) )
1561                   ENDDO
1562                ENDDO
1563             ENDIF
1564          ENDIF
1565          IF ( myidx == id_stg_right  )  THEN
1566             DO  j = nys, nyn
1567                DO  k = nzb+1, nzt
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 ) )
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!
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
1585       IF ( stg_call )  THEN
1586          DO  i = nxl, nxr
1587             DO  k = nzb, nzt + 1
1588!
1589!--             Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1590                IF ( tu(k) == 0.0_wp  .OR.  adjustment_step )  THEN
1591                   fu_xz(k,i) = fuo_xz(k,i)
1592                ELSE
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) ) )
1595                ENDIF
1596
1597                IF ( tv(k) == 0.0_wp  .OR.  adjustment_step )  THEN
1598                   fv_xz(k,i) = fvo_xz(k,i)
1599                ELSE
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) ) )
1602                ENDIF
1603
1604                IF ( tw(k) == 0.0_wp  .OR.  adjustment_step )  THEN
1605                   fw_xz(k,i) = fwo_xz(k,i)
1606                ELSE
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) ) )
1609                ENDIF
1610             ENDDO
1611          ENDDO
1612
1613
1614          dist_xz(nzb,:,1) = 0.0_wp
1615          dist_xz(nzb,:,2) = 0.0_wp
1616          dist_xz(nzb,:,3) = 0.0_wp
1617
1618          IF ( myidy == id_stg_south  ) j = nys
1619          IF ( myidy == id_stg_north )  j = nyn+1
1620          DO  i = nxl, nxr
1621             DO  k = nzb+1, nzt + 1
1622!
1623!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
1624!--             Additional factors are added to improve the variance of v and w
1625                !experimental test of 1.2
1626                dist_xz(k,i,2) = ( a21(k) * fu_xz(k,i) + a22(k) * fv_xz(k,i) ) *                   &
1627                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
1628             ENDDO
1629          ENDDO
1630
1631          IF ( myidy == id_stg_south  ) j = nys-1
1632          IF ( myidy == id_stg_north )  j = nyn+1
1633          DO  i = nxl, nxr
1634             DO  k = nzb+1, nzt + 1
1635!
1636!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
1637!--             Additional factors are added to improve the variance of v and w
1638                dist_xz(k,i,1) = a11(k) * fu_xz(k,i) *                                             &
1639                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
1640
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) ) *                                         &
1643                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
1644             ENDDO
1645          ENDDO
1646       ENDIF
1647
1648!
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.
1653       IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
1654          mc_factor_l = 0.0_wp
1655          mc_factor   = 0.0_wp
1656
1657          IF ( myidy == id_stg_south  ) j = nys
1658          IF ( myidy == id_stg_north )  j = nyn+1
1659
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 ) ) )
1662
1663          IF ( myidy == id_stg_south ) j = nys-1
1664          IF ( myidy == id_stg_north ) j = nyn+1
1665
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 ) ) )
1670
1671#if defined( __parallel )
1672          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 3, MPI_REAL, MPI_SUM, comm1dx, ierr )
1673#else
1674          mc_factor   = mc_factor_l
1675#endif
1676
1677          mc_factor = mc_factor / REAL( nr_non_topo_xz, KIND = wp )
1678
1679          IF ( myidy == id_stg_south ) j = nys
1680          IF ( myidy == id_stg_north ) j = nyn+1
1681
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 ) )
1684
1685
1686          IF ( myidy == id_stg_south ) j = nys-1
1687          IF ( myidy == id_stg_north ) j = nyn+1
1688
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 ) )
1691
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 ) )
1694!
1695!--       Add disturbances
1696          IF ( myidy == id_stg_south )  THEN
1697             DO  i = nxl, nxr
1698                DO  k = nzb+1, nzt
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 ) )
1703                   v(k,-1,i) = v(k,0,i)
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 ) )
1706                ENDDO
1707             ENDDO
1708          ENDIF
1709          IF ( myidy == id_stg_north  )  THEN
1710
1711             DO  i = nxl, nxr
1712                DO  k = nzb+1, nzt
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 ) )
1719                ENDDO
1720             ENDDO
1721          ENDIF
1722       ENDIF
1723    ENDIF
1724!
1725!-- Exchange ghost points.
1726    CALL exchange_horiz( u, nbgp )
1727    CALL exchange_horiz( v, nbgp )
1728    CALL exchange_horiz( w, nbgp )
1729!
1730!-- Finally, set time counter for calling STG to zero
1731    IF ( stg_call )  time_stg_call = 0.0_wp
1732!
1733!-- Set adjustment step to False to indicate that time correlation can be
1734!-- switched-on again.
1735    adjustment_step = .FALSE.
1736
1737    IF ( debug_output_timestep )  CALL debug_message( 'stg_main', 'end' )
1738
1739 END SUBROUTINE stg_main
1740
1741!--------------------------------------------------------------------------------------------------!
1742! Description:
1743! ------------
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!--------------------------------------------------------------------------------------------------!
1749 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id_left, id_right )
1750
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
1758
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
1761
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
1765
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
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
1773!
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
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
1780!-- correlation matrix.
1781    ALLOCATE( rand_it(nzb-mergp_z:nzt+1+mergp_z,-mergp_y+nys:nyn+mergp_y) )
1782    rand_it = 0.0_wp
1783
1784    rand_av        = 0.0_wp
1785    rand_sigma_inv = 0.0_wp
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 )
1788!
1789!-- Compute and normalize random numbers.
1790    DO  j = nys - mergp_y, nyn + mergp_y
1791!
1792!--    Put the random seeds at grid point j
1793       CALL random_seed_parallel( put=seq_rand_yz(:,j) )
1794       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
1795          CALL random_number_parallel( random_dummy )
1796          rand_it(k,j) = random_dummy
1797       ENDDO
1798!
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
1802!
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.
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
1808          rand_av = rand_av + rand_it(k,j)
1809       ENDDO
1810    ENDDO
1811
1812#if defined( __parallel )
1813!
1814!-- Sum-up the local averages of the random numbers
1815    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1816#endif
1817    rand_av = rand_av * nyz_inv
1818!
1819!-- Obtain zero mean
1820    rand_it= rand_it - rand_av
1821!
1822!-- Now, compute the variance
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
1825          rand_sigma_inv = rand_sigma_inv + rand_it(k,j)**2
1826       ENDDO
1827    ENDDO
1828
1829#if defined( __parallel )
1830!
1831!-- Sum-up the local quadratic averages of the random numbers
1832    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1833#endif
1834!
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
1840    ENDIF
1841!
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' )
1846!
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.
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
1866       ENDDO
1867
1868    ELSE
1869
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
1877             ENDDO
1878          ENDDO
1879       ENDDO
1880!
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
1884
1885#if defined( __parallel )
1886!
1887!--    Gather the velocity seed matrix on left boundary mpi ranks.
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 )
1890!
1891!--    If required, gather the same velocity seed matrix on right boundary mpi ranks (in offline
1892!--    nesting for example).
1893       IF ( PRESENT( id_right ) )  THEN
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 )
1896       ENDIF
1897#else
1898       f_n(nzb+1:nzt+1,nys:nyn) = f_n_l(nzb_x_stg:nzt_x_stg+1,nys:nyn)
1899!
1900!--    Next line required to avoid compile errors because of unused dummy arguments
1901       IF ( id_left == 0 )  id_right = 0
1902#endif
1903
1904    ENDIF
1905
1906    DEALLOCATE( rand_it )
1907
1908    CALL cpu_log( log_point_s(31), 'STG f_n factors', 'stop' )
1909
1910 END SUBROUTINE stg_generate_seed_yz
1911
1912
1913!--------------------------------------------------------------------------------------------------!
1914! Description:
1915! ------------
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!--------------------------------------------------------------------------------------------------!
1922 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id_south, id_north )
1923
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
1931
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
1934
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
1938
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
1941
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
1944
1945    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it   !< global array of random numbers
1946
1947!
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
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
1954!-- correlation matrix.
1955    ALLOCATE( rand_it(nzb-mergp_z:nzt+1+mergp_z,-mergp_x+nxl:nxr+mergp_x) )
1956    rand_it = 0.0_wp
1957
1958    rand_av        = 0.0_wp
1959    rand_sigma_inv = 0.0_wp
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 )
1962!
1963!-- Compute and normalize random numbers.
1964    DO  i = nxl - mergp_x, nxr + mergp_x
1965!
1966!--    Put the random seeds at grid point ii
1967       CALL random_seed_parallel( put=seq_rand_xz(:,i) )
1968       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
1969          CALL random_number_parallel( random_dummy )
1970          rand_it(k,i) = random_dummy
1971       ENDDO
1972!
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
1976!
1977!-- For normalization to zero mean, sum-up the global random numers.
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
1980!-- also occur on the inner subdomains.
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
1983          rand_av = rand_av + rand_it(k,i)
1984       ENDDO
1985    ENDDO
1986
1987#if defined( __parallel )
1988!
1989!-- Sum-up the local averages of the random numbers
1990    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1991#endif
1992    rand_av = rand_av * nxz_inv
1993!
1994!-- Obtain zero mean
1995    rand_it= rand_it - rand_av
1996!
1997!-- Now, compute the variance
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
2000          rand_sigma_inv = rand_sigma_inv + rand_it(k,i)**2
2001       ENDDO
2002    ENDDO
2003
2004#if defined( __parallel )
2005!
2006!-- Sum-up the local quadratic averages of the random numbers
2007    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL, MPI_SUM, comm1dx, ierr )
2008#endif
2009!
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
2015    ENDIF
2016!
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' )
2021!
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.
2029    IF ( compute_velocity_seeds_local )  THEN
2030
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
2038             ENDDO
2039          ENDDO
2040       ENDDO
2041
2042    ELSE
2043
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
2054!
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
2058
2059#if defined( __parallel )
2060!
2061!--    Gather the processed velocity seed on south boundary mpi ranks.
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 )
2064!
2065!--    Gather the processed velocity seed on north boundary mpi ranks.
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 )
2068#else
2069       f_n(nzb+1:nzt+1,nxl:nxr) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxl:nxr)
2070!
2071!--    Next line required to avoid compile errors because of unused dummy arguments
2072       IF ( id_north == 0 )  id_south = 0
2073#endif
2074
2075    ENDIF
2076
2077    DEALLOCATE( rand_it )
2078
2079    CALL cpu_log( log_point_s(31), 'STG f_n factors', 'stop' )
2080
2081 END SUBROUTINE stg_generate_seed_xz
2082
2083!--------------------------------------------------------------------------------------------------!
2084! Description:
2085! ------------
2086!> Parametrization of the Reynolds-stress componentes, turbulent length- and time scales. The
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.
2089!--------------------------------------------------------------------------------------------------!
2090  SUBROUTINE parametrize_turbulence
2091
2092    INTEGER(iwp) :: k                 !< loop index in z-direction
2093
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
2113
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
2118!
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 )
2124!
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,                                                                  &
2148                          scale_l < -5.0_wp )
2149!
2150!--    Determine normalized height coordinate
2151       zzi = zu(k) * d_zi
2152!
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 )
2158
2159       r11(k) = ABS(   cos_phi(k) * rlon1 + sin_phi(k) * rlat1 )
2160       r22(k) = ABS( - sin_phi(k) * rlon1 + cos_phi(k) * rlat1 )
2161!
2162!--    w'w'
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                              )
2167!
2168!--    u'w' and v'w'. After calculation of the longitudinal and crosswind component
2169!--    these are projected along the x- and y-direction. Note, it is assumed that
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) ) )
2178!
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!
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
2206    ENDDO
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)
2213
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
2224!
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
2235       r21(k) = 0.5_wp * ( r31(k) + r32(k) )
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
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
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)
2280
2281    tu(nzb)  = tu(nzb+1)
2282    tv(nzb)  = tv(nzb+1)
2283    tw(nzb)  = tw(nzb+1)
2284
2285    adjustment_step = .TRUE.
2286
2287 END SUBROUTINE parametrize_turbulence
2288
2289!--------------------------------------------------------------------------------------------------!
2290! Description:
2291! ------------
2292!> Calculate the coefficient matrix from the Lund rotation.
2293!--------------------------------------------------------------------------------------------------!
2294 SUBROUTINE calc_coeff_matrix
2295
2296    INTEGER(iwp) ::  k  !< loop index in z-direction
2297
2298!
2299!-- Calculate coefficient matrix. Split loops to allow for loop vectorization.
2300    DO  k = nzb+1, nzt+1
2301       IF ( r11(k) > 10E-6_wp )  THEN
2302          a11(k) = SQRT( r11(k) )
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
2312       a22(k) = r22(k) - a21(k)**2
2313       IF ( a22(k) > 10E-6_wp )  THEN
2314          a22(k) = SQRT( a22(k) )
2315          a32(k) = ( r32(k) - a21(k) * a31(k) ) / a22(k)
2316       ELSE
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
2323       IF ( a33(k) > 10E-6_wp )  THEN
2324          a33(k) =  SQRT( a33(k) )
2325       ELSE
2326          a33(k) = 10E-8_wp
2327       ENDIF
2328    ENDDO
2329!
2330!-- Set bottom boundary condition
2331    a11(nzb) = a11(nzb+1)
2332    a22(nzb) = a22(nzb+1)
2333    a21(nzb) = a21(nzb+1)
2334    a33(nzb) = a33(nzb+1)
2335    a31(nzb) = a31(nzb+1)
2336    a32(nzb) = a32(nzb+1)
2337
2338 END SUBROUTINE calc_coeff_matrix
2339
2340!--------------------------------------------------------------------------------------------------!
2341! Description:
2342! ------------
2343!> This routine controls the re-adjustment of the turbulence statistics used for generating
2344!> turbulence at the lateral boundaries.
2345!--------------------------------------------------------------------------------------------------!
2346 SUBROUTINE stg_adjust
2347
2348    IF ( debug_output_timestep )  CALL debug_message( 'stg_adjust', 'start' )
2349!
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
2355!
2356!-- Compute scaling parameters (domain-averaged), such as friction velocity are calculated.
2357    CALL calc_scaling_variables
2358!
2359!-- Parametrize Reynolds-stress tensor. Parametrization follows Brost et al. (1982) with
2360!-- modifications described in Rotach et al. (1996) and is based on boundary-layer depth, friction
2361!-- velocity and velocity scale.
2362    CALL parametrize_turbulence
2363!
2364!-- Calculate coefficient matrix from Reynolds stress tensor
2365!-- (Lund rotation)
2366    CALL calc_coeff_matrix
2367!
2368!-- Setup filter functions according to updated length scales.
2369    CALL stg_setup_filter_function
2370!
2371!-- Reset time counter for controlling next adjustment to zero
2372    time_stg_adjust = 0.0_wp
2373
2374    IF ( debug_output_timestep )  CALL debug_message( 'stg_adjust', 'end' )
2375
2376 END SUBROUTINE stg_adjust
2377
2378
2379!--------------------------------------------------------------------------------------------------!
2380! Description:
2381! ------------
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!--------------------------------------------------------------------------------------------------!
2386 SUBROUTINE calc_scaling_variables
2387
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
2392
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
2401!
2402!-- Calculate mean friction velocity, velocity scale, heat flux and near-surface temperature in the
2403!-- model domain.
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
2408    DO  m = 1, surf_def_h(0)%ns
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)
2416    ENDDO
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)
2424       pt_surf_mean_l = pt_surf_mean_l  + pt(k,j,i)
2425    ENDDO
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)
2433       pt_surf_mean_l = pt_surf_mean_l  + pt(k,j,i)
2434    ENDDO
2435
2436#if defined( __parallel )
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 )
2441#else
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
2446#endif
2447
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
2452!
2453!-- Compute mean convective velocity scale. Note, in case the mean heat flux is negative, set
2454!-- convective velocity scale to zero.
2455    IF ( shf_mean > 0.0_wp )  THEN
2456       w_convective = ( g * shf_mean * zi_ribulk / pt_surf_mean )**( 1.0_wp / 3.0_wp )
2457    ELSE
2458       w_convective = 0.0_wp
2459    ENDIF
2460!
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!
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).
2468    scale_wm = ( scale_us**3 + 0.6_wp * w_convective**3 )**( 1.0_wp / 3.0_wp )
2469
2470 END SUBROUTINE calc_scaling_variables
2471
2472 END MODULE synthetic_turbulence_generator_mod
Note: See TracBrowser for help on using the repository browser.