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

Last change on this file since 4440 was 4440, checked in by suehring, 3 years ago

Correct misplaced preprocessor directive

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