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

Last change on this file since 3891 was 3891, checked in by suehring, 5 years ago

Bugfix in initialization of turbulence generator in case of restart runs and offline nesting

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