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

Last change on this file since 4329 was 4329, checked in by motisi, 18 months ago

Renamed wall_flags_0 to wall_flags_static_0

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