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

Last change on this file since 4441 was 4441, checked in by suehring, 19 months ago

Change order of dimension in surface arrays %frac, %emissivity and %albedo to allow for better vectorization in the radiation interactions; Set back turbulent length scale to 8 x grid spacing in the parametrized mode for the synthetic turbulence generator (was accidentally changed in last commit)

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