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

Last change on this file since 4028 was 4022, checked in by suehring, 5 years ago

Synthetic turbulence generator: Revise bias correction of imposed perturbations (correction via volume flow can create instabilities in case the mean volume flow is close to zero); Introduce lower limits in calculation of coefficient matrix, else the calculation may become numerically unstable; Impose perturbations every timestep, even though no new set of perturbations is generated in case dt_stg_call /= dt_3d; Implement a gradual decrease of Reynolds stress and length scales above ABL height (within 1 length scale above ABL depth to 1/10) rather than an abrupt decline; Bugfix in non-nested case: use ABL height for parametrized turbulence; Offline nesting: Rename subroutine for ABL height calculation and make it public; Change top boundary condition for pressure back to Neumann

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