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

Last change on this file since 4148 was 4148, checked in by suehring, 23 months ago

indoor_model_mod: Revision of indoor model; urban_surface_mod: parameters for waste heat from cooling and heating are introduced to building data base; initialization of building data base moved to an earlier point of time before indoor model will be initialized; radiation_model_mod: minor improvement in some comments; synthetic_turbulence_generator: unused variable removed

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