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

Last change on this file since 4598 was 4566, checked in by suehring, 10 months ago

synthetic turbulence generator: revise parametrization for reynolds-stress components, turbulent length- and time scales; revise computation of velocity disturbances to be consistent to Xie and Castro (2008); change default value of time interval to adjust turbulence parametrization; bugfix in computation of amplitude-tensor (vertical flux of horizontal momentum)

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