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

Last change on this file since 4335 was 4335, checked in by suehring, 16 months ago

Plant canopy: implement fix for LAD on building edges also for the vectorl-optimized branch

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