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

Last change on this file since 4603 was 4603, checked in by suehring, 4 years ago

Synthetic turbulence generator: Bugfix in initialization from ASCII file - x-length scales at the bottom boundary were not initialized properly

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