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

Last change on this file since 4442 was 4442, checked in by suehring, 3 years ago

last commit documented

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