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

Last change on this file since 4444 was 4444, checked in by raasch, 3 years ago

bugfix: cpp-directives for serial mode added

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