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

Last change on this file since 4438 was 4438, checked in by suehring, 19 months ago

Synthetic turbulence: performance optimizations - random numbers only defined and computed locally, option to compute velocity seeds locally without need of global communication; paralell random number generator: new routine to initialize 1D random number arrays; virtual measurements: CPU-log points added

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