source: palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90 @ 4649

Last change on this file since 4649 was 4647, checked in by suehring, 4 years ago

Change default value of synthetic turbulence adjustment as well as compute_velocity_seeds_local. By default, the random-seed computation is now distributed among several cores. Especially for large length scales this is significantly faster.

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