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

Last change on this file since 4463 was 4457, checked in by raasch, 4 years ago

ghost point exchange modularized, bugfix for wrong 2d-exchange

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