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

Last change on this file since 4332 was 4332, checked in by suehring, 4 years ago

Synthetic turbulence: Limit initial velocity seeds in restart runs, if not the seed calculation may become unstable.

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