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

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

Synthetic turbulence generator: Computation of velocity seeds optimized. This implies that random numbers are computed now using the parallel random number generator. Random number are now only computed and normalized locally, while distributed over all mpi ranks afterwards, instead of computing random numbers on a global array. urther, the number of calls for the time-consuming velocity-seed generation is reduced - now the left and right, as well as the north and south boundary share the same velocity-seed matrices.

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