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

Last change on this file since 4495 was 4495, checked in by raasch, 5 years ago

restart data handling with MPI-IO added, first part

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