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

Last change on this file since 4629 was 4629, checked in by raasch, 4 years ago

support for MPI Fortran77 interface (mpif.h) removed

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