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

Last change on this file since 4403 was 4346, checked in by motisi, 4 years ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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