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

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

synthetic turbulence generator: revise parametrization for reynolds-stress components, turbulent length- and time scales; revise computation of velocity disturbances to be consistent to Xie and Castro (2008); change default value of time interval to adjust turbulence parametrization; bugfix in computation of amplitude-tensor (vertical flux of horizontal momentum)

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