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

Last change on this file since 4071 was 4071, checked in by suehring, 2 years ago

Bugfix in synthetic turbulence generator in non-nested cases (offline or self nesting) when no cyclic fill is used: in order to get correct initial inflow profiles at the left boundary the array mean_inflow_profiles need to be initialized accordingly.

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