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

Last change on this file since 4698 was 4671, checked in by pavelkrc, 4 years ago

Radiative transfer model RTM version 4.1

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