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

Last change on this file since 3049 was 3049, checked in by Giersch, 3 years ago

Revision history corrected

  • Property svn:keywords set to Id
File size: 62.1 KB
Line 
1!> @synthetic_turbulence_generator_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: synthetic_turbulence_generator_mod.f90 3049 2018-05-29 13:52:36Z Giersch $
27! Error messages revised
28!
29! 3045 2018-05-28 07:55:41Z Giersch
30! Error message revised
31!
32! 3044 2018-05-25 10:59:41Z gronemeier
33! Add missing variable descriptions
34!
35! 3038 2018-05-24 10:54:00Z gronemeier
36! updated variable descriptions
37!
38! 2967 2018-04-13 11:22:08Z raasch
39! bugfix: missing parallel cpp-directives added
40!
41! 2946 2018-04-04 17:01:23Z suehring
42! Remove unused module load
43!
44! 2945 2018-04-04 16:27:14Z suehring
45! - Bugfix in parallelization of synthetic turbulence generator in case the
46!   z-dimension is not an integral divisor of the number of processors along
47!   the x- and y-dimension
48! - Revision in control mimic in case of RAN-LES nesting
49!
50! 2938 2018-03-27 15:52:42Z suehring
51! Apply turbulence generator at all non-cyclic lateral boundaries in case of
52! realistic Inifor large-scale forcing or RANS-LES nesting
53!
54! 2936 2018-03-27 14:49:27Z suehring
55! variable named found has been introduced for checking if restart data was found,
56! reading of restart strings has been moved completely to read_restart_data_mod,
57! redundant skipping function has been removed, stg_read/write_restart_data
58! have been renamed to stg_r/wrd_global, stg_rrd_global is called in
59! read_restart_data_mod now, flag syn_turb_gen_prerun and marker *** end stg
60! *** have been removed (Giersch), strings and their respective lengths are
61! written out and read now in case of restart runs to get rid of prescribed
62! character lengths (Giersch), CASE DEFAULT was added if restart data is read
63!
64! 2841 2018-02-27 15:02:57Z suehring
65! Bugfix: wrong placement of include 'mpif.h' corrected
66!
67! 2836 2018-02-26 13:40:05Z Giersch
68! The variables synthetic_turbulence_generator and
69! use_synthetic_turbulence_generator have been abbreviated + syn_turb_gen_prerun
70! flag is used to define if module related parameters were outputted as restart
71! data
72!
73! 2716 2017-12-29 16:35:59Z kanani
74! Corrected "Former revisions" section
75!
76! 2696 2017-12-14 17:12:51Z kanani
77! Change in file header (GPL part)
78!
79! 2669 2017-12-06 16:03:27Z raasch
80! unit number for file containing turbulence generator data changed to 90
81! bugfix: preprocessor directives added for MPI specific code
82!
83! 2576 2017-10-24 13:49:46Z Giersch
84! Definition of a new function called stg_skip_global to skip module
85! parameters during reading restart data
86!
87! 2563 2017-10-19 15:36:10Z Giersch
88! stg_read_restart_data is called in stg_parin in the case of a restart run
89!
90! 2259 2017-06-08 09:09:11Z gronemeier
91! Initial revision
92!
93!
94!
95! Authors:
96! --------
97! @author Tobias Gronemeier, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
98!
99!
100! Description:
101! ------------
102!> The module generates turbulence at the inflow boundary based on a method by
103!> Xie and Castro (2008) utilizing a Lund rotation (Lund, 1998) and a mass-flux
104!> correction by Kim et al. (2013).
105!> The turbulence is correlated based on length scales in y- and z-direction and
106!> a time scale for each velocity component. The profiles of length and time
107!> scales, mean u, v, w, e and pt, and all components of the Reynolds stress
108!> tensor are read from file STG_PROFILES.
109!>
110!> @todo test restart
111!>       enable cyclic_fill
112!>       implement turbulence generation for e and pt
113!> @todo Input of height-constant length scales via namelist
114!> @note <Enter notes on the module>
115!> @bug  Height information from input file is not used. Profiles from input
116!>       must match with current PALM grid.
117!>       Transformation of length scales to number of gridpoints does not
118!>       consider grid stretching.
119!>       In case of restart, velocity seeds differ from precursor run if a11,
120!>       a22, or a33 are zero.
121!------------------------------------------------------------------------------!
122 MODULE synthetic_turbulence_generator_mod
123
124
125    USE arrays_3d,                                                             &
126        ONLY:  mean_inflow_profiles, u, v, w
127
128    USE constants,                                                             &
129        ONLY:  pi
130
131    USE control_parameters,                                                    &
132        ONLY:  initializing_actions, message_string, syn_turb_gen
133
134    USE cpulog,                                                                &
135        ONLY:  cpu_log, log_point, log_point_s
136
137    USE indices,                                                               &
138        ONLY:  nbgp, nzb, nzt, nxl, nxlg, nxr, nxrg, nys, nyn, nyng, nysg
139
140    USE kinds
141
142#if defined( __parallel )  &&  !defined( __mpifh )
143    USE MPI
144#endif
145
146    USE pegrid,                                                                &
147        ONLY:  comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
148
149    USE transpose_indices,                                                     &
150        ONLY: nzb_x, nzt_x
151
152
153    IMPLICIT NONE
154
155#if defined( __parallel )  &&  defined( __mpifh )
156    INCLUDE "mpif.h"
157#endif
158
159    LOGICAL :: velocity_seed_initialized = .FALSE.  !< true after first call of stg_main
160    LOGICAL :: use_syn_turb_gen = .FALSE.           !< switch to use synthetic turbulence generator
161
162    INTEGER(iwp) ::  id_stg_left        !< left lateral boundary core id in case of turbulence generator
163    INTEGER(iwp) ::  id_stg_north       !< north lateral boundary core id in case of turbulence generator
164    INTEGER(iwp) ::  id_stg_right       !< right lateral boundary core id in case of turbulence generator
165    INTEGER(iwp) ::  id_stg_south       !< south lateral boundary core id in case of turbulence generator
166    INTEGER(iwp) ::  stg_type_xz        !< MPI type for full z range
167    INTEGER(iwp) ::  stg_type_xz_small  !< MPI type for small z range
168    INTEGER(iwp) ::  stg_type_yz        !< MPI type for full z range
169    INTEGER(iwp) ::  stg_type_yz_small  !< MPI type for small z range
170    INTEGER(iwp) ::  merg               !< maximum length scale (in gp)
171    INTEGER(iwp) ::  mergp              !< merg + nbgp
172    INTEGER(iwp) ::  nzb_x_stg          !< lower bound of z coordinate (required for transposing z on PEs along x)
173    INTEGER(iwp) ::  nzt_x_stg          !< upper bound of z coordinate (required for transposing z on PEs along x)
174    INTEGER(iwp) ::  nzb_y_stg          !< lower bound of z coordinate (required for transposing z on PEs along y)
175    INTEGER(iwp) ::  nzt_y_stg          !< upper bound of z coordinate (required for transposing z on PEs along y)
176
177    REAL(wp) :: mc_factor    !< mass flux correction factor
178
179    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_xz      !< displacement for MPI_GATHERV
180    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_xz  !< receive count for MPI_GATHERV
181    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_yz      !< displacement for MPI_GATHERV
182    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_yz  !< receive count for MPI_GATHERV
183    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nux            !< length scale of u in x direction (in gp)
184    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy            !< length scale of u in y direction (in gp)
185    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz            !< length scale of u in z direction (in gp)
186    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvx            !< length scale of v in x direction (in gp)
187    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy            !< length scale of v in y direction (in gp)
188    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz            !< length scale of v in z direction (in gp)
189    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwx            !< length scale of w in x direction (in gp)
190    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy            !< length scale of w in y direction (in gp)
191    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz            !< length scale of w in z direction (in gp)
192
193    INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed        !< seed of random number for rn-generator
194
195    REAL(wp), DIMENSION(:), ALLOCATABLE :: a11             !< coefficient for Lund rotation
196    REAL(wp), DIMENSION(:), ALLOCATABLE :: a21             !< coefficient for Lund rotation
197    REAL(wp), DIMENSION(:), ALLOCATABLE :: a22             !< coefficient for Lund rotation
198    REAL(wp), DIMENSION(:), ALLOCATABLE :: a31             !< coefficient for Lund rotation
199    REAL(wp), DIMENSION(:), ALLOCATABLE :: a32             !< coefficient for Lund rotation
200    REAL(wp), DIMENSION(:), ALLOCATABLE :: a33             !< coefficient for Lund rotation
201    REAL(wp), DIMENSION(:), ALLOCATABLE :: tu              !< Lagrangian time scale of u
202    REAL(wp), DIMENSION(:), ALLOCATABLE :: tv              !< Lagrangian time scale of v
203    REAL(wp), DIMENSION(:), ALLOCATABLE :: tw              !< Lagrangian time scale of w
204
205    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bux           !< filter function for u in x direction
206    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy           !< filter function for u in y direction
207    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz           !< filter function for u in z direction
208    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvx           !< filter function for v in x direction
209    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy           !< filter function for v in y direction
210    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz           !< filter function for v in z direction
211    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwx           !< filter function for w in y direction
212    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy           !< filter function for w in y direction
213    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz           !< filter function for w in z direction
214    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_xz         !< velocity seed for u at xz plane
215    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_xz        !< velocity seed for u at xz plane with new random number
216    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_yz         !< velocity seed for u at yz plane
217    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_yz        !< velocity seed for u at yz plane with new random number
218    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_xz         !< velocity seed for v at xz plane
219    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_xz        !< velocity seed for v at xz plane with new random number
220    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_yz         !< velocity seed for v at yz plane
221    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_yz        !< velocity seed for v at yz plane with new random number
222    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_xz         !< velocity seed for w at xz plane
223    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_xz        !< velocity seed for w at xz plane with new random number
224    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_yz         !< velocity seed for w at yz plane
225    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz        !< velocity seed for w at yz plane with new random number
226
227
228!
229!-- PALM interfaces:
230!-- Input parameter checks to be done in check_parameters
231    INTERFACE stg_check_parameters
232       MODULE PROCEDURE stg_check_parameters
233    END INTERFACE stg_check_parameters
234
235!
236!-- Calculate filter functions
237    INTERFACE stg_filter_func
238       MODULE PROCEDURE stg_filter_func
239    END INTERFACE stg_filter_func
240
241!
242!-- Generate velocity seeds at south and north domain boundary
243    INTERFACE stg_generate_seed_xz
244       MODULE PROCEDURE stg_generate_seed_xz
245    END INTERFACE stg_generate_seed_xz
246!
247!-- Generate velocity seeds at left and/or right domain boundary
248    INTERFACE stg_generate_seed_yz
249       MODULE PROCEDURE stg_generate_seed_yz
250    END INTERFACE stg_generate_seed_yz
251
252!
253!-- Output of information to the header file
254    INTERFACE stg_header
255       MODULE PROCEDURE stg_header
256    END INTERFACE stg_header
257
258!
259!-- Initialization actions
260    INTERFACE stg_init
261       MODULE PROCEDURE stg_init
262    END INTERFACE stg_init
263
264!
265!-- Main procedure of synth. turb. gen.
266    INTERFACE stg_main
267       MODULE PROCEDURE stg_main
268    END INTERFACE stg_main
269
270!
271!-- Reading of NAMELIST parameters
272    INTERFACE stg_parin
273       MODULE PROCEDURE stg_parin
274    END INTERFACE stg_parin
275
276!
277!-- Reading of parameters for restart runs
278    INTERFACE stg_rrd_global
279       MODULE PROCEDURE stg_rrd_global
280    END INTERFACE stg_rrd_global
281
282!
283!-- Writing of binary output for restart runs
284    INTERFACE stg_wrd_global
285       MODULE PROCEDURE stg_wrd_global
286    END INTERFACE stg_wrd_global
287
288    SAVE
289
290    PRIVATE
291
292!
293!-- Public interfaces
294    PUBLIC  stg_check_parameters, stg_header, stg_init, stg_main, stg_parin,   &
295            stg_wrd_global, stg_rrd_global
296
297!
298!-- Public variables
299    PUBLIC  id_stg_left, id_stg_north, id_stg_right, id_stg_south,             &
300            use_syn_turb_gen
301
302
303 CONTAINS
304
305
306!------------------------------------------------------------------------------!
307! Description:
308! ------------
309!> Check parameters routine for synthetic turbulence generator
310!------------------------------------------------------------------------------!
311 SUBROUTINE stg_check_parameters
312
313
314    USE control_parameters,                                                    &
315        ONLY:  bc_lr, bc_ns, forcing, nest_domain, rans_mode, turbulent_inflow
316
317    USE pmc_interface,                                                         &
318        ONLY : rans_mode_parent
319
320
321    IMPLICIT NONE
322
323    IF ( .NOT. use_syn_turb_gen  .AND.  .NOT. rans_mode  .AND.  forcing )  THEN
324       message_string = 'Synthetic turbulence generator has to be applied ' // &
325                        'when forcing is used and model operates in LES mode.'
326       CALL message( 'stg_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
327    ENDIF
328
329    IF ( .NOT. use_syn_turb_gen  .AND.  nest_domain                            &
330         .AND. rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
331       message_string = 'Synthetic turbulence generator has to be applied ' // &
332                        'when nesting is applied and parent operates in '  //  &
333                        'RANS-mode but current child in LES mode.'
334       CALL message( 'stg_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
335    ENDIF
336
337    IF ( use_syn_turb_gen )  THEN
338
339       IF ( .NOT. forcing  .AND.  .NOT. nest_domain )  THEN
340
341          IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0     &
342        .AND.  INDEX( initializing_actions, 'read_restart_data' ) == 0 )  THEN
343             message_string = 'Using synthetic turbulence generator ' //       &
344                              'requires %initializing_actions = '         //   &
345                              '"set_constant_profiles" or "read_restart_data"'
346             CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
347          ENDIF
348
349          IF ( bc_lr /= 'dirichlet/radiation' )  THEN
350             message_string = 'Using synthetic turbulence generator ' //       &
351                              'requires &bc_lr = "dirichlet/radiation"'
352             CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
353          ENDIF
354          IF ( bc_ns /= 'cyclic' )  THEN
355             message_string = 'Using synthetic turbulence generator ' //       &
356                              'requires &bc_ns = "cyclic"'
357             CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
358          ENDIF
359
360       ENDIF
361
362       IF ( turbulent_inflow )  THEN
363          message_string = 'Using synthetic turbulence generator ' //          &
364                           'in combination &with turbulent_inflow = .T. '//    &
365                              'is not allowed'
366          CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
367       ENDIF
368
369    ENDIF
370
371 END SUBROUTINE stg_check_parameters
372
373
374!------------------------------------------------------------------------------!
375! Description:
376! ------------
377!> Header output for synthetic turbulence generator
378!------------------------------------------------------------------------------!
379 SUBROUTINE stg_header ( io )
380
381
382    IMPLICIT NONE
383
384    INTEGER(iwp), INTENT(IN) ::  io   !< Unit of the output file
385
386!
387!-- Write synthetic turbulence generator Header
388    WRITE( io, 1 )
389    IF ( use_syn_turb_gen )  THEN
390       WRITE( io, 2 )
391    ELSE
392       WRITE( io, 3 )
393    ENDIF
394
3951   FORMAT (//' Synthetic turbulence generator information:'/                  &
396              ' ------------------------------------------'/)
3972   FORMAT ('    synthetic turbulence generator switched on')
3983   FORMAT ('    synthetic turbulence generator switched off')
399
400 END SUBROUTINE stg_header
401
402
403!------------------------------------------------------------------------------!
404! Description:
405! ------------
406!> Initialization of the synthetic turbulence generator
407!------------------------------------------------------------------------------!
408 SUBROUTINE stg_init
409
410
411    USE arrays_3d,                                                             &
412        ONLY:  ddzw, u_init, v_init, zu
413
414    USE control_parameters,                                                    &
415        ONLY:  coupling_char, dz, e_init, forcing, nest_domain, rans_mode
416
417    USE grid_variables,                                                        &
418        ONLY:  ddx, ddy
419
420    USE indices,                                                               &
421        ONLY:  nz
422
423    USE pmc_interface,                                                         &
424        ONLY : rans_mode_parent
425
426
427    IMPLICIT NONE
428
429    LOGICAL ::  file_stg_exist = .FALSE. !< flag indication whether parameter file for Reynolds stress and length scales exist
430
431#if defined( __parallel )
432    INTEGER(KIND=MPI_ADDRESS_KIND) :: extent !< extent of new MPI type
433    INTEGER(KIND=MPI_ADDRESS_KIND) :: tob=0  !< dummy variable
434#endif
435
436    INTEGER(iwp) :: i                        !> grid index in x-direction
437    INTEGER(iwp) :: j                        !> loop index
438    INTEGER(iwp) :: k                        !< index
439    INTEGER(iwp) :: newtype                  !< dummy MPI type
440    INTEGER(iwp) :: realsize                 !< size of REAL variables
441    INTEGER(iwp) :: nseed                    !< dimension of random number seed
442    INTEGER(iwp) :: startseed = 1234567890   !< start seed for random number generator
443
444!
445!-- Dummy variables used for reading profiles
446    REAL(wp) :: d1      !< u profile
447    REAL(wp) :: d2      !< v profile
448    REAL(wp) :: d3      !< w profile
449    REAL(wp) :: d5      !< e profile
450    REAL(wp) :: d11     !< vertical interpolation for a11
451    REAL(wp) :: d21     !< vertical interpolation for a21
452    REAL(wp) :: d22     !< vertical interpolation for a22
453    REAL(wp) :: dum_exp !< dummy variable used for exponential vertical decrease of turbulent length scales
454    REAL(wp) :: luy     !< length scale for u in y direction
455    REAL(wp) :: luz     !< length scale for u in z direction
456    REAL(wp) :: lvy     !< length scale for v in y direction
457    REAL(wp) :: lvz     !< length scale for v in z direction
458    REAL(wp) :: lwy     !< length scale for w in y direction
459    REAL(wp) :: lwz     !< length scale for w in z direction
460    REAL(wp) :: nnz     !< increment used to determine processor decomposition of z-axis along x and y direction
461    REAL(wp) :: zz      !< height
462
463    REAL(wp) :: length_scale_surface !< typical length scale
464    REAL(wp) :: r_ii_0               !< correction factor
465    REAL(wp) :: time_scale           !< typical time scale
466    REAL(wp) :: length_scale_z       !< typical length scale
467
468    REAL(wp),DIMENSION(nzb:nzt+1) :: r11  !< Reynolds parameter
469    REAL(wp),DIMENSION(nzb:nzt+1) :: r21  !< Reynolds parameter
470    REAL(wp),DIMENSION(nzb:nzt+1) :: r22  !< Reynolds parameter
471    REAL(wp),DIMENSION(nzb:nzt+1) :: r31  !< Reynolds parameter
472    REAL(wp),DIMENSION(nzb:nzt+1) :: r32  !< Reynolds parameter
473    REAL(wp),DIMENSION(nzb:nzt+1) :: r33  !< Reynolds parameter
474
475
476#if defined( __parallel )
477    CALL MPI_BARRIER( comm2d, ierr )
478#endif
479
480    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
481
482    IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )                             &
483       ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )
484
485    ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                 &
486               a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1),                 &
487               nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1),                 &
488               nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1),                 &
489               nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
490               tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
491
492#if defined( __parallel )
493!
494!-- Determine processor decomposition of z-axis along x- and y-direction
495    nnz = nz / pdims(1)
496    nzb_x_stg = 1 + myidx * INT( nnz )
497    nzt_x_stg = ( myidx + 1 ) * INT( nnz )
498
499    IF ( MOD( nz , pdims(1) ) /= 0  .AND.  myidx == id_stg_right )             &
500       nzt_x_stg = nzt_x_stg + myidx * ( nnz - INT( nnz ) )
501!        nzt_x_stg = myidx * nnz + MOD( nz , pdims(1) )
502
503    IF ( forcing  .OR.  ( nest_domain .AND.  rans_mode_parent  .AND.           &
504                   .NOT.  rans_mode ) )  THEN
505       nnz = nz / pdims(2)
506       nzb_y_stg = 1 + myidy * INT( nnz )
507       nzt_y_stg = ( myidy + 1 ) * INT( nnz )
508
509       IF ( MOD( nz , pdims(2) ) /= 0  .AND.  myidy == id_stg_north )          &
510          nzt_y_stg = nzt_y_stg + myidy * ( nnz - INT( nnz ) )
511!           nzt_y_stg = myidy * nnz + MOD( nz , pdims(2) )
512    ENDIF
513
514!
515!-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted
516!-- velocity seeds
517    CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
518    extent = 1 * realsize
519!
520!-- Set-up MPI datatyp to involve all cores for turbulence generation at yz
521!-- layer
522!-- stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
523    CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyng-nysg+1],                 &
524            [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
525    CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz, ierr )
526    CALL MPI_TYPE_COMMIT( stg_type_yz, ierr )
527    CALL MPI_TYPE_FREE( newtype, ierr )
528
529    ! stg_type_yz_small: yz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
530    CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyng-nysg+1],     &
531            [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
532    CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr )
533    CALL MPI_TYPE_COMMIT( stg_type_yz_small, ierr )
534    CALL MPI_TYPE_FREE( newtype, ierr )
535
536    ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
537    ALLOCATE( recv_count_yz(pdims(1)), displs_yz(pdims(1)) )
538
539    recv_count_yz           = nzt_x_stg-nzb_x_stg + 1
540    recv_count_yz(pdims(1)) = recv_count_yz(pdims(1)) + 1
541
542    DO  j = 1, pdims(1)
543       displs_yz(j) = 0 + (nzt_x_stg-nzb_x_stg+1) * (j-1)
544    ENDDO
545!
546!-- Set-up MPI datatyp to involve all cores for turbulence generation at xz
547!-- layer
548!-- stg_type_xz: xz-slice with vertical bounds nzb:nzt+1
549    IF ( forcing  .OR.  ( nest_domain .AND.  rans_mode_parent  .AND.           &
550                   .NOT.  rans_mode ) )  THEN
551       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nxrg-nxlg+1],              &
552               [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
553       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz, ierr )
554       CALL MPI_TYPE_COMMIT( stg_type_xz, ierr )
555       CALL MPI_TYPE_FREE( newtype, ierr )
556
557       ! stg_type_yz_small: xz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
558       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_y_stg-nzb_y_stg+2,nxrg-nxlg+1],  &
559               [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
560       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz_small, ierr )
561       CALL MPI_TYPE_COMMIT( stg_type_xz_small, ierr )
562       CALL MPI_TYPE_FREE( newtype, ierr )
563
564       ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
565       ALLOCATE( recv_count_xz(pdims(2)), displs_xz(pdims(2)) )
566
567       recv_count_xz           = nzt_y_stg-nzb_y_stg + 1
568       recv_count_xz(pdims(2)) = recv_count_xz(pdims(2)) + 1
569
570       DO  j = 1, pdims(2)
571          displs_xz(j) = 0 + (nzt_y_stg-nzb_y_stg+1) * (j-1)
572       ENDDO
573
574    ENDIF
575
576#endif
577!
578!-- Define seed of random number
579    CALL RANDOM_SEED()
580    CALL RANDOM_SEED( size=nseed )
581    ALLOCATE( seed(1:nseed) )
582    DO  j = 1, nseed
583       seed(j) = startseed + j
584    ENDDO
585    CALL RANDOM_SEED(put=seed)
586
587!-- Read inflow profile
588!-- Height levels of profiles in input profile is as follows:
589!-- zu: luy, luz, tu, lvy, lvz, tv, r11, r21, r22, d1, d2, d5
590!-- zw: lwy, lwz, tw, r31, r32, r33, d3
591!-- WARNING: zz is not used at the moment
592    INQUIRE( FILE = 'STG_PROFILES' // TRIM( coupling_char ),                   &
593             EXIST = file_stg_exist  )
594
595    IF ( file_stg_exist )  THEN
596
597       OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD',     &
598                      FORM='FORMATTED')
599!
600!--    Skip header
601       READ( 90, * )
602
603       DO  k = nzb, nzt+1
604          READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k), &
605                        r11(k), r21(k), r22(k), r31(k), r32(k), r33(k),        &
606                        d1, d2, d3, d5
607
608!
609!--       Convert length scales from meter to number of grid points
610          nuy(k) = INT( luy * ddy )
611          nuz(k) = INT( luz / dz  )
612          nvy(k) = INT( lvy * ddy )
613          nvz(k) = INT( lvz / dz  )
614          nwy(k) = INT( lwy * ddy )
615          nwz(k) = INT( lwz / dz  )
616!
617!--       Workaround, assume isotropic turbulence
618          nwx(k) = nwy(k)
619          nvx(k) = nvy(k)
620          nux(k) = nuy(k)
621!
622!--       Save Mean inflow profiles
623          IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
624             mean_inflow_profiles(k,1) = d1
625             mean_inflow_profiles(k,2) = d2
626            !  mean_inflow_profiles(k,4) = d4
627             mean_inflow_profiles(k,5) = d5
628          ENDIF
629       ENDDO
630
631       CLOSE( 90 )
632
633    ELSE
634!
635!--    Set-up defaul length scales. Assume exponentially decreasing length
636!--    scales and isotropic turbulence.
637!--    Typical length (time) scales of 100 m (s) should be a good compromise
638!--    between all stratrifications. Near-surface variances are fixed to
639!--    0.1 m2/s2, vertical fluxes are one order of magnitude smaller.
640!--    Vertical fluxes
641       length_scale_surface = 100.0_wp
642       r_ii_0               = 0.1_wp
643       time_scale           = 100.0_wp
644       DO  k = nzb+1, nzt+1
645          dum_exp        = MERGE( -zu(k) / ( 0.3* zu(nzt) ),                   &
646                                  0.0_wp,                                      &
647                                  zu(k) > 0.3 * zu(nzt)                        &
648                                )
649          length_scale_z = length_scale_surface * EXP( dum_exp )
650
651          nux(k) = MAX( INT( length_scale_z * ddx     ), 1 )
652          nuy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
653          nuz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
654          nvx(k) = MAX( INT( length_scale_z * ddx     ), 1 )
655          nvy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
656          nvz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
657          nwx(k) = MAX( INT( length_scale_z * ddx     ), 1 )
658          nwy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
659          nwz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
660
661          r11(k) = r_ii_0 * EXP( dum_exp )
662          r22(k) = r_ii_0 * EXP( dum_exp )
663          r33(k) = r_ii_0 * EXP( dum_exp )
664
665          r21(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
666          r31(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
667          r32(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
668
669          tu(k)  = time_scale
670          tv(k)  = time_scale
671          tw(k)  = time_scale
672
673       ENDDO
674       nux(nzb) = nux(nzb+1)
675       nuy(nzb) = nuy(nzb+1)
676       nuz(nzb) = nuz(nzb+1)
677       nvx(nzb) = nvx(nzb+1)
678       nvy(nzb) = nvy(nzb+1)
679       nvz(nzb) = nvz(nzb+1)
680       nwx(nzb) = nwx(nzb+1)
681       nwy(nzb) = nwy(nzb+1)
682       nwz(nzb) = nwz(nzb+1)
683
684       r11(nzb) = r11(nzb+1)
685       r22(nzb) = r22(nzb+1)
686       r33(nzb) = r33(nzb+1)
687
688       r21(nzb) = r11(nzb+1)
689       r31(nzb) = r31(nzb+1)
690       r32(nzb) = r32(nzb+1)
691
692       tu(nzb)  = time_scale
693       tv(nzb)  = time_scale
694       tw(nzb)  = time_scale
695
696    ENDIF
697
698!
699!-- Assign initial profiles
700    IF ( .NOT. forcing  .AND.  .NOT.  nest_domain )  THEN
701       u_init = mean_inflow_profiles(:,1)
702       v_init = mean_inflow_profiles(:,2)
703      !pt_init = mean_inflow_profiles(:,4)
704       e_init = MAXVAL( mean_inflow_profiles(:,5) )
705    ENDIF
706!
707!-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
708    DO  k = nzb, nzt+1
709       IF ( r11(k) > 0.0_wp )  THEN
710          a11(k) = SQRT( r11(k) )
711          a21(k) = r21(k) / a11(k)
712       ELSE
713          a11(k) = 0.0_wp
714          a21(k) = 0.0_wp
715       ENDIF
716
717       a22(k) = r22(k) - a21(k)**2
718       IF ( a22(k) > 0.0_wp )  THEN
719          a22(k) = SQRT( a22(k) )
720          a32(k) = ( r32(k) - a21(k) * a31(k) ) / a22(k)
721       ELSE
722          a22(k) = 0.0_wp
723          a32(k) = 0.0_wp
724       ENDIF
725
726!
727!--    a31, a32, a33 must be calculated with interpolated a11, a21, a22 (d11,
728!--    d21, d22) because of different vertical grid
729       IF ( k .le. nzt )  THEN
730          d11 = 0.5_wp * ( r11(k) + r11(k+1) )
731          IF ( d11 > 0.0_wp )  THEN
732             d11 = SQRT( d11 )
733             d21 = ( 0.5_wp * ( r21(k) + r21(k+1) ) ) / d11
734             a31(k) = r31(k) / d11
735          ELSE
736             d21 = 0.0_wp
737             a31(k) = 0.0_wp
738          ENDIF
739
740          d22 = 0.5_wp * ( r22(k) + r22(k+1) ) - d21 ** 2
741          IF ( d22 > 0.0_wp )  THEN
742             a32(k) = ( r32(k) - d21 * a31(k) ) / SQRT( d22 )
743          ELSE
744             a32(k) = 0.0_wp
745          ENDIF
746
747          a33(k) = r33(k) - a31(k) ** 2 - a32(k) ** 2
748          IF ( a33(k) > 0.0_wp )  THEN
749             a33(k) = SQRT( a33(k) )
750          ELSE
751             a33(k) = 0.0_wp
752          ENDIF
753       ELSE
754          a31(k) = a31(k-1)
755          a32(k) = a32(k-1)
756          a33(k) = a33(k-1)
757       ENDIF
758
759    ENDDO
760!
761!-- Define the size of the filter functions and allocate them.
762    merg = 0
763
764    ! arrays must be large enough to cover the largest length scale
765    DO  k = nzb, nzt+1
766       j = MAX( ABS(nux(k)), ABS(nuy(k)), ABS(nuz(k)), &
767                ABS(nvx(k)), ABS(nvy(k)), ABS(nvz(k)), &
768                ABS(nwx(k)), ABS(nwy(k)), ABS(nwz(k))  )
769       IF ( j > merg )  merg = j
770    ENDDO
771
772    merg  = 2 * merg
773    mergp = merg + nbgp
774
775    ALLOCATE ( bux(-merg:merg,nzb:nzt+1),                                      &
776               buy(-merg:merg,nzb:nzt+1),                                      &
777               buz(-merg:merg,nzb:nzt+1),                                      &
778               bvx(-merg:merg,nzb:nzt+1),                                      &
779               bvy(-merg:merg,nzb:nzt+1),                                      &
780               bvz(-merg:merg,nzb:nzt+1),                                      &
781               bwx(-merg:merg,nzb:nzt+1),                                      &
782               bwy(-merg:merg,nzb:nzt+1),                                      &
783               bwz(-merg:merg,nzb:nzt+1)  )
784
785!
786!-- Allocate velocity seeds for turbulence at xz-layer
787    ALLOCATE ( fu_xz( nzb:nzt+1,nxlg:nxrg), fuo_xz(nzb:nzt+1,nxlg:nxrg),       &
788               fv_xz( nzb:nzt+1,nxlg:nxrg), fvo_xz(nzb:nzt+1,nxlg:nxrg),       &
789               fw_xz( nzb:nzt+1,nxlg:nxrg), fwo_xz(nzb:nzt+1,nxlg:nxrg)  )
790
791!
792!-- Allocate velocity seeds for turbulence at yz-layer
793    ALLOCATE ( fu_yz( nzb:nzt+1,nysg:nyng), fuo_yz(nzb:nzt+1,nysg:nyng),       &
794               fv_yz( nzb:nzt+1,nysg:nyng), fvo_yz(nzb:nzt+1,nysg:nyng),       &
795               fw_yz( nzb:nzt+1,nysg:nyng), fwo_yz(nzb:nzt+1,nysg:nyng)  )
796
797    fu_xz  = 0.0_wp
798    fuo_xz = 0.0_wp
799    fv_xz  = 0.0_wp
800    fvo_xz = 0.0_wp
801    fw_xz  = 0.0_wp
802    fwo_xz = 0.0_wp
803
804    fu_yz  = 0.0_wp
805    fuo_yz = 0.0_wp
806    fv_yz  = 0.0_wp
807    fvo_yz = 0.0_wp
808    fw_yz  = 0.0_wp
809    fwo_yz = 0.0_wp
810
811!
812!-- Create filter functions
813    CALL stg_filter_func( nux, bux ) !filter ux
814    CALL stg_filter_func( nuy, buy ) !filter uy
815    CALL stg_filter_func( nuz, buz ) !filter uz
816    CALL stg_filter_func( nvx, bvx ) !filter vx
817    CALL stg_filter_func( nvy, bvy ) !filter vy
818    CALL stg_filter_func( nvz, bvz ) !filter vz
819    CALL stg_filter_func( nwx, bwx ) !filter wx
820    CALL stg_filter_func( nwy, bwy ) !filter wy
821    CALL stg_filter_func( nwz, bwz ) !filter wz
822
823#if defined( __parallel )
824    CALL MPI_BARRIER( comm2d, ierr )
825#endif
826
827!
828!-- In case of restart, calculate velocity seeds fu, fv, fw from former
829!   time step.
830!   Bug: fu, fv, fw are different in those heights where a11, a22, a33
831!        are 0 compared to the prerun. This is mostly for k=nzt+1.
832    IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
833       IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
834
835          IF ( myidx == id_stg_left  )  i = -1
836          IF ( myidx == id_stg_right )  i = nxr+1
837
838          DO  j = nysg, nyng
839             DO  k = nzb, nzt+1
840
841                IF  ( a11(k) .NE. 0._wp ) THEN
842                   fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
843                ELSE
844                   fu_yz(k,j) = 0._wp
845                ENDIF
846
847                IF  ( a22(k) .NE. 0._wp ) THEN
848                   fv_yz(k,j) = ( v(k,j,i) / mc_factor - a21(k) * fu_yz(k,j) - &
849                               v_init(k) ) / a22(k)
850                ELSE
851                   fv_yz(k,j) = 0._wp
852                ENDIF
853
854                IF  ( a33(k) .NE. 0._wp ) THEN
855                   fw_yz(k,j) = ( w(k,j,i) / mc_factor - a31(k) * fu_yz(k,j) - &
856                               a32(k) * fv_yz(k,j) ) / a33(k)
857                ELSE
858                   fw_yz = 0._wp
859                ENDIF
860
861             ENDDO
862          ENDDO
863       ENDIF
864    ENDIF
865
866    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
867
868 END SUBROUTINE stg_init
869
870
871!------------------------------------------------------------------------------!
872! Description:
873! ------------
874!> Calculate filter function bxx from length scale nxx following Eg.9 and 10
875!> (Xie and Castro, 2008)
876!------------------------------------------------------------------------------!
877 SUBROUTINE stg_filter_func( nxx, bxx )
878
879
880    IMPLICIT NONE
881
882    INTEGER(iwp) :: k         !< loop index
883    INTEGER(iwp) :: n_k       !< length scale nXX in height k
884    INTEGER(iwp) :: n_k2      !< n_k * 2
885    INTEGER(iwp) :: nf        !< index for length scales
886
887    REAL(wp) :: bdenom        !< denominator for filter functions bXX
888    REAL(wp) :: qsi = 1.0_wp  !< minimization factor
889
890    INTEGER(iwp), DIMENSION(:) :: nxx(nzb:nzt+1)           !< length scale (in gp)
891
892    REAL(wp), DIMENSION(:,:) :: bxx(-merg:merg,nzb:nzt+1)  !< filter function
893
894
895    bxx = 0.0_wp
896
897    DO  k = nzb, nzt+1
898       bdenom = 0.0_wp
899       n_k    = nxx(k)
900       IF ( n_k /= 0 )  THEN
901          n_k2 = n_k * 2
902
903!
904!--       ( Eq.10 )^2
905          DO  nf = -n_k2, n_k2
906             bdenom = bdenom + EXP( -qsi * pi * ABS(nf) / n_k )**2
907          ENDDO
908
909!
910!--       ( Eq.9 )
911          bdenom = SQRT( bdenom )
912          DO  nf = -n_k2, n_k2
913             bxx(nf,k) = EXP( -qsi * pi * ABS(nf) / n_k ) / bdenom
914          ENDDO
915       ENDIF
916    ENDDO
917
918END SUBROUTINE stg_filter_func
919
920
921!------------------------------------------------------------------------------!
922! Description:
923! ------------
924!> Parin for &stg_par for synthetic turbulence generator
925!------------------------------------------------------------------------------!
926 SUBROUTINE stg_parin
927
928
929    IMPLICIT NONE
930
931    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
932
933
934    NAMELIST /stg_par/   use_syn_turb_gen
935
936    line = ' '
937
938!
939!-- Try to find stg package
940    REWIND ( 11 )
941    line = ' '
942    DO WHILE ( INDEX( line, '&stg_par' ) == 0 )
943       READ ( 11, '(A)', END=10 )  line
944    ENDDO
945    BACKSPACE ( 11 )
946
947!
948!-- Read namelist
949    READ ( 11, stg_par )
950
951!
952!-- Set flag that indicates that the synthetic turbulence generator is switched
953!-- on
954    syn_turb_gen = .TRUE.
955
956
957 10 CONTINUE
958
959 END SUBROUTINE stg_parin
960
961
962!------------------------------------------------------------------------------!
963! Description:
964! ------------
965!> This routine reads the respective restart data.
966!------------------------------------------------------------------------------!
967 SUBROUTINE stg_rrd_global( found )
968
969
970    USE control_parameters,                                                    &
971        ONLY: length, restart_string
972
973
974    IMPLICIT NONE
975
976    LOGICAL, INTENT(OUT)  ::  found !< flag indicating if variable was found
977
978
979    found = .TRUE.
980
981
982    SELECT CASE ( restart_string(1:length) )
983
984       CASE ( 'mc_factor' )
985          READ ( 13 )  mc_factor
986       CASE ( 'use_syn_turb_gen' )
987          READ ( 13 )  use_syn_turb_gen
988
989       CASE DEFAULT
990
991          found = .FALSE.
992
993    END SELECT
994
995
996 END SUBROUTINE stg_rrd_global
997
998
999!------------------------------------------------------------------------------!
1000! Description:
1001! ------------
1002!> This routine writes the respective restart data.
1003!------------------------------------------------------------------------------!
1004 SUBROUTINE stg_wrd_global
1005
1006
1007    IMPLICIT NONE
1008
1009    CALL wrd_write_string( 'mc_factor' )
1010    WRITE ( 14 )  mc_factor
1011
1012    CALL wrd_write_string( 'use_syn_turb_gen' )
1013    WRITE ( 14 )  use_syn_turb_gen
1014
1015
1016 END SUBROUTINE stg_wrd_global
1017
1018
1019!------------------------------------------------------------------------------!
1020! Description:
1021! ------------
1022!> Create turbulent inflow fields for u, v, w with prescribed length scales and
1023!> Reynolds stress tensor after a method of Xie and Castro (2008), modified
1024!> following suggestions of Kim et al. (2013), and using a Lund rotation
1025!> (Lund, 1998).
1026!------------------------------------------------------------------------------!
1027 SUBROUTINE stg_main
1028
1029
1030    USE arrays_3d,                                                             &
1031        ONLY:  dzw
1032
1033    USE control_parameters,                                                    &
1034        ONLY:  dt_3d, forcing, intermediate_timestep_count,  nest_domain,      &
1035               rans_mode, simulated_time, volume_flow_initial
1036
1037    USE grid_variables,                                                        &
1038        ONLY:  dx, dy
1039
1040    USE indices,                                                               &
1041        ONLY:  wall_flags_0
1042
1043    USE statistics,                                                            &
1044        ONLY:  weight_substep
1045
1046    USE pmc_interface,                                                         &
1047        ONLY : rans_mode_parent
1048
1049
1050    IMPLICIT NONE
1051
1052    INTEGER(iwp) :: i           !< grid index in x-direction
1053    INTEGER(iwp) :: j           !< loop index in y-direction
1054    INTEGER(iwp) :: k           !< loop index in z-direction
1055
1056    REAL(wp) :: dt_stg          !< wheighted subtimestep
1057    REAL(wp) :: mc_factor_l     !< local mass flux correction factor
1058    REAL(wp) :: volume_flow     !< mass flux through lateral boundary
1059    REAL(wp) :: volume_flow_l   !< local mass flux through lateral boundary
1060
1061    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg,5) :: dist_xz !< imposed disturbances at north/south boundary
1062    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5) :: dist_yz !< imposed disturbances at left/right boundary
1063
1064
1065    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1066
1067!
1068!-- Calculate time step which is needed for filter functions
1069    dt_stg = dt_3d * weight_substep(intermediate_timestep_count)
1070
1071!
1072!-- Initial value of fu, fv, fw
1073    IF ( simulated_time == 0.0_wp .AND. .NOT. velocity_seed_initialized )  THEN
1074       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
1075       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
1076       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
1077
1078       IF ( forcing  .OR.  ( nest_domain .AND.  rans_mode_parent  .AND.        &
1079                      .NOT.  rans_mode ) )  THEN
1080!
1081!--       Generate turbulence at right boundary
1082          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_right )
1083          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_right )
1084          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_right )
1085!
1086!--       Generate turbulence at north boundary
1087          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_north )
1088          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_north )
1089          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_north )
1090!
1091!--       Generate turbulence at south boundary
1092          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_south )
1093          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_south )
1094          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_south )
1095       ENDIF
1096       velocity_seed_initialized = .TRUE.
1097    ENDIF
1098!
1099!-- New set of fu, fv, fw
1100    CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
1101    CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
1102    CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
1103
1104    IF ( forcing  .OR.  ( nest_domain .AND.  rans_mode_parent  .AND.           &
1105                   .NOT.  rans_mode ) )  THEN
1106!
1107!--       Generate turbulence at right boundary
1108          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
1109          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
1110          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
1111!
1112!--       Generate turbulence at north boundary
1113          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
1114          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
1115          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
1116!
1117!--       Generate turbulence at south boundary
1118          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
1119          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
1120          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
1121    ENDIF
1122!
1123!-- Turbulence generation at left and or right boundary
1124    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
1125
1126       DO  j = nysg, nyng
1127          DO  k = nzb, nzt + 1
1128!
1129!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1130             IF ( tu(k) == 0.0_wp )  THEN
1131                fu_yz(k,j) = fuo_yz(k,j)
1132             ELSE
1133                fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
1134                         fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
1135             ENDIF
1136
1137             IF ( tv(k) == 0.0_wp )  THEN
1138                fv_yz(k,j) = fvo_yz(k,j)
1139             ELSE
1140                fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
1141                         fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
1142             ENDIF
1143
1144             IF ( tw(k) == 0.0_wp )  THEN
1145                fw_yz(k,j) = fwo_yz(k,j)
1146             ELSE
1147                fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
1148                         fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
1149             ENDIF
1150!
1151!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
1152!--          Additional factors are added to improve the variance of v and w
1153             IF( k == 0 )  THEN
1154                dist_yz(k,j,1) = 0.0_wp
1155                dist_yz(k,j,2) = 0.0_wp
1156                dist_yz(k,j,3) = 0.0_wp
1157!                 dist_yz(k,j,4) = 0.0_wp
1158!                 dist_yz(k,j,5) = 0.0_wp
1159             ELSE
1160                dist_yz(k,j,1) = a11(k) * fu_yz(k,j)
1161                !experimental test of 1.2
1162                dist_yz(k,j,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
1163                                         * 1.2_wp )                            &
1164                                       * (   a21(k) * fu_yz(k,j)               &
1165                                           + a22(k) * fv_yz(k,j) )
1166                dist_yz(k,j,3) = ( SQRT(a33(k) / MAXVAL(a33) )                 &
1167                                         * 1.3_wp )                            &
1168                                       * (   a31(k) * fu_yz(k,j)               &
1169                                           + a32(k) * fv_yz(k,j)               &
1170                                           + a33(k) * fw_yz(k,j) )
1171                ! Calculation for pt and e not yet implemented
1172!                 dist_yz(k,j,4) = 0.0_wp
1173!                 dist_yz(k,j,5) = 0.0_wp
1174             ENDIF
1175
1176          ENDDO
1177       ENDDO
1178
1179!
1180!--    Mass flux correction following Kim et al. (2013)
1181!--    This correction factor insures that the mass flux is preserved at the
1182!--    inflow boundary
1183       IF ( .NOT. forcing  .AND.  .NOT. nest_domain )  THEN
1184          mc_factor_l = 0.0_wp
1185          mc_factor   = 0.0_wp
1186          DO  j = nys, nyn
1187             DO  k = nzb+1, nzt
1188                mc_factor_l = mc_factor_l + dzw(k)  *                          &
1189                              ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) )
1190             ENDDO
1191          ENDDO
1192
1193#if defined( __parallel )
1194          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,  &
1195                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1196#else
1197          mc_factor = mc_factor_l
1198#endif
1199
1200          mc_factor = volume_flow_initial(1) / mc_factor
1201
1202!
1203!--       Add disturbance at the inflow
1204          DO  j = nysg, nyng
1205             DO  k = nzb, nzt+1
1206                 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +              &
1207                                      dist_yz(k,j,1)             ) * mc_factor
1208                 v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +              &
1209                                      dist_yz(k,j,2)             ) * mc_factor
1210                 w(k,j,-nbgp:-1)  =   dist_yz(k,j,3)               * mc_factor
1211             ENDDO
1212          ENDDO
1213
1214       ELSE
1215!
1216!--       First, calculate volume flow at yz boundary
1217          IF ( myidx == id_stg_left  )  i = nxl
1218          IF ( myidx == id_stg_right )  i = nxr+1
1219
1220          volume_flow_l = 0.0_wp
1221          volume_flow   = 0.0_wp
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                volume_flow_l = volume_flow_l + u(k,j,i) * dzw(k) * dy         &
1227                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1228                                              BTEST( wall_flags_0(k,j,i), 1 ) )
1229
1230                mc_factor_l = mc_factor_l     + ( u(k,j,i) + dist_yz(k,j,1) )  &
1231                                                         * dzw(k) * dy         &
1232                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1233                                              BTEST( wall_flags_0(k,j,i), 1 ) )
1234             ENDDO
1235          ENDDO
1236#if defined( __parallel )
1237          CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                      &
1238                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1239          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                          &
1240                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1241#else
1242          volume_flow = volume_flow_l
1243          mc_factor   = mc_factor_l
1244#endif
1245
1246          mc_factor = volume_flow / mc_factor
1247
1248!
1249!--       Add disturbances
1250          IF ( myidx == id_stg_left  )  THEN
1251
1252             DO  j = nysg, nyng
1253                DO  k = nzb, nzt+1
1254                   u(k,j,-nbgp+1:0) = ( u(k,j,-nbgp+1:0) + dist_yz(k,j,1) )    &
1255                                        * mc_factor
1256                   v(k,j,-nbgp:-1)  = ( v(k,j,-nbgp:-1)  + dist_yz(k,j,2)  )   &
1257                                        * mc_factor
1258                   w(k,j,-nbgp:-1)  = ( w(k,j,-nbgp:-1)  + dist_yz(k,j,3)  )   &
1259                                        * mc_factor
1260                ENDDO
1261             ENDDO
1262          ENDIF
1263          IF ( myidx == id_stg_right  )  THEN
1264
1265             DO  j = nysg, nyng
1266                DO  k = nzb, nzt+1
1267                   u(k,j,nxr+1:nxr+nbgp) = ( u(k,j,nxr+1:nxr+nbgp) +           &
1268                                             dist_yz(k,j,1) ) * mc_factor
1269                   v(k,j,nxr+1:nxr+nbgp) = ( v(k,j,nxr+1:nxr+nbgp) +           &
1270                                             dist_yz(k,j,2) ) * mc_factor
1271                   w(k,j,nxr+1:nxr+nbgp) = ( w(k,j,nxr+1:nxr+nbgp) +           &
1272                                             dist_yz(k,j,3) ) * mc_factor
1273                ENDDO
1274             ENDDO
1275          ENDIF
1276       ENDIF
1277
1278    ENDIF
1279!
1280!-- Turbulence generation at north and south boundary
1281    IF ( myidy == id_stg_north  .OR.  myidy == id_stg_south )  THEN
1282
1283       DO  i = nxlg, nxrg
1284          DO  k = nzb, nzt + 1
1285!
1286!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1287             IF ( tu(k) == 0.0_wp )  THEN
1288                fu_xz(k,i) = fuo_xz(k,i)
1289             ELSE
1290                fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
1291                         fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
1292             ENDIF
1293
1294             IF ( tv(k) == 0.0_wp )  THEN
1295                fv_xz(k,i) = fvo_xz(k,i)
1296             ELSE
1297                fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
1298                         fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
1299             ENDIF
1300
1301             IF ( tw(k) == 0.0_wp )  THEN
1302                fw_xz(k,i) = fwo_xz(k,i)
1303             ELSE
1304                fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
1305                         fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
1306             ENDIF
1307!
1308!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
1309!--          Additional factors are added to improve the variance of v and w
1310             IF( k == 0 )  THEN
1311                dist_xz(k,i,1) = 0.0_wp
1312                dist_xz(k,i,2) = 0.0_wp
1313                dist_xz(k,i,3) = 0.0_wp
1314
1315             ELSE
1316                dist_xz(k,i,1) = a11(k) * fu_xz(k,i)
1317                !experimental test of 1.2
1318                dist_xz(k,i,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
1319                                         * 1.2_wp )                            &
1320                                       * (   a21(k) * fu_xz(k,i)               &
1321                                           + a22(k) * fv_xz(k,i) )
1322                dist_xz(k,i,3) = ( SQRT(a33(k) / MAXVAL(a33) )                 &
1323                                         * 1.3_wp )                            &
1324                                       * (   a31(k) * fu_xz(k,i)               &
1325                                           + a32(k) * fv_xz(k,i)               &
1326                                           + a33(k) * fw_xz(k,i) )
1327             ENDIF
1328
1329          ENDDO
1330       ENDDO
1331!
1332!--    Mass flux correction following Kim et al. (2013)
1333!--    This correction factor insures that the mass flux is preserved at the
1334!--    inflow boundary.
1335!--    First, calculate volume flow at xz boundary
1336       IF ( myidy == id_stg_south  ) j = nys
1337       IF ( myidy == id_stg_north )  j = nyn+1
1338
1339       volume_flow_l = 0.0_wp
1340       volume_flow   = 0.0_wp
1341       mc_factor_l   = 0.0_wp
1342       mc_factor     = 0.0_wp
1343       DO  i = nxl, nxr
1344          DO  k = nzb+1, nzt
1345             volume_flow_l = volume_flow_l + v(k,j,i) * dzw(k) * dx            &
1346                                  * MERGE( 1.0_wp, 0.0_wp,                     &
1347                                           BTEST( wall_flags_0(k,j,i), 2 ) )
1348
1349             mc_factor_l = mc_factor_l     + ( v(k,j,i) + dist_xz(k,i,2) )     &
1350                                                      * dzw(k) * dx            &
1351                                  * MERGE( 1.0_wp, 0.0_wp,                     &
1352                                           BTEST( wall_flags_0(k,j,i), 2 ) )
1353          ENDDO
1354       ENDDO
1355#if defined( __parallel )
1356       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                         &
1357                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1358       CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                             &
1359                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1360#else
1361       volume_flow = volume_flow_l
1362       mc_factor   = mc_factor_l
1363#endif
1364
1365       mc_factor = volume_flow / mc_factor
1366
1367!
1368!--    Add disturbances
1369       IF ( myidy == id_stg_south  )  THEN
1370
1371          DO  i = nxlg, nxrg
1372             DO  k = nzb, nzt+1
1373                u(k,-nbgp:-1,i) = ( u(k,-nbgp:-1,i) + dist_xz(k,i,1) )         &
1374                                     * mc_factor
1375                v(k,-nbgp:0,i)  = ( v(k,-nbgp:0,i)  + dist_xz(k,i,2)  )        &
1376                                     * mc_factor
1377                w(k,-nbgp:-1,i) = ( w(k,-nbgp:-1,i) + dist_xz(k,i,3)  )        &
1378                                     * mc_factor
1379             ENDDO
1380          ENDDO
1381       ENDIF
1382       IF ( myidy == id_stg_north  )  THEN
1383
1384          DO  i = nxlg, nxrg
1385             DO  k = nzb, nzt+1
1386                u(k,nyn+1:nyn+nbgp,i) = ( u(k,nyn+1:nyn+nbgp,i) +              &
1387                                          dist_xz(k,i,1) ) * mc_factor
1388                v(k,nyn+1:nyn+nbgp,i) = ( v(k,nyn+1:nyn+nbgp,i) +              &
1389                                          dist_xz(k,i,2) ) * mc_factor
1390                w(k,nyn+1:nyn+nbgp,i) = ( w(k,nyn+1:nyn+nbgp,i) +              &
1391                                          dist_xz(k,i,3) ) * mc_factor
1392             ENDDO
1393          ENDDO
1394       ENDIF
1395    ENDIF
1396
1397    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1398
1399 END SUBROUTINE stg_main
1400
1401!------------------------------------------------------------------------------!
1402! Description:
1403! ------------
1404!> Generate a set of random number rand_it wich is equal on each PE
1405!> and calculate the velocity seed f_n.
1406!> f_n is splitted in vertical direction by the number of PEs in x-direction and
1407!> and each PE calculates a vertical subsection of f_n. At the the end, all
1408!> parts are collected to form the full array.
1409!------------------------------------------------------------------------------!
1410 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
1411
1412
1413    USE indices,                                                               &
1414        ONLY: ny
1415
1416    IMPLICIT NONE
1417
1418    INTEGER(iwp) :: id          !< core ids at respective boundaries
1419    INTEGER(iwp) :: j           !< loop index in y-direction
1420    INTEGER(iwp) :: jj          !< loop index in y-direction
1421    INTEGER(iwp) :: k           !< loop index in z-direction
1422    INTEGER(iwp) :: kk          !< loop index in z-direction
1423    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
1424
1425    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y    !< length scale in y-direction
1426    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
1427    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y2   !< n_y*2
1428    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
1429
1430    REAL(wp) :: nyz_inv         !< inverse of number of grid points in yz-slice
1431    REAL(wp) :: rand_av         !< average of random number
1432    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
1433
1434    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter func in y-dir
1435    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
1436    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
1437    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)     :: f_n     !<  velocity seed
1438    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
1439
1440
1441!
1442!-- Generate random numbers using a seed generated in stg_init.
1443!-- The set of random numbers are modified to have an average of 0 and
1444!-- unit variance.
1445    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:ny+mergp) )
1446
1447    rand_av        = 0.0_wp
1448    rand_sigma_inv = 0.0_wp
1449    nyz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( ny+1 ), KIND=wp )
1450
1451    DO  j = 0, ny
1452       DO  k = nzb, nzt+1
1453          CALL RANDOM_NUMBER( rand_it(k,j) )
1454          rand_av = rand_av + rand_it(k,j)
1455       ENDDO
1456    ENDDO
1457
1458    rand_av = rand_av * nyz_inv
1459
1460    DO  j = 0, ny
1461       DO  k = nzb, nzt+1
1462          rand_it(k,j)   = rand_it(k,j) - rand_av
1463          rand_sigma_inv = rand_sigma_inv + rand_it(k,j) ** 2
1464       ENDDO
1465    ENDDO
1466
1467    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nyz_inv)
1468
1469    DO  j = 0, ny
1470       DO  k = nzb, nzt+1
1471          rand_it(k,j) = rand_it(k,j) * rand_sigma_inv
1472       ENDDO
1473    ENDDO
1474
1475!
1476!-- Periodic fill of random number in space
1477    DO  j = 0, ny
1478       DO  k = 1, mergp
1479          rand_it(nzb  -k,j) = rand_it(nzt+2-k,j)    ! bottom margin
1480          rand_it(nzt+1+k,j) = rand_it(nzb+k-1,j)    ! top margin
1481       ENDDO
1482    ENDDO
1483    DO  j = 1, mergp
1484       DO  k = nzb-mergp, nzt+1+mergp
1485          rand_it(k,  -j) = rand_it(k,ny-j+1)        ! south margin
1486          rand_it(k,ny+j) = rand_it(k,   j-1)        ! north margin
1487       ENDDO
1488    ENDDO
1489
1490!
1491!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1492    n_y2 = n_y * 2
1493    n_z2 = n_z * 2
1494    f_n_l  = 0.0_wp
1495
1496    DO  j = nysg, nyng
1497       DO  k = nzb_x_stg, nzt_x_stg+1
1498          DO  jj = -n_y2(k), n_y2(k)
1499             DO  kk = -n_z2(k), n_z2(k)
1500                f_n_l(k,j) = f_n_l(k,j)                                        &
1501                           + b_y(jj,k) * b_z(kk,k) * rand_it(k+kk,j+jj)
1502             ENDDO
1503          ENDDO
1504       ENDDO
1505    ENDDO
1506
1507    DEALLOCATE( rand_it )
1508!
1509!-- Gather velocity seeds of full subdomain
1510    send_count = nzt_x_stg - nzb_x_stg + 1
1511    IF ( nzt_x_stg == nzt )  send_count = send_count + 1
1512
1513#if defined( __parallel )
1514    CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,    &
1515                      f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,  &
1516                      id, comm1dx, ierr )
1517#else
1518    f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x_stg:nzt_x_stg+1,nysg:nyng)
1519#endif
1520
1521
1522 END SUBROUTINE stg_generate_seed_yz
1523
1524
1525
1526
1527!------------------------------------------------------------------------------!
1528! Description:
1529! ------------
1530!> Generate a set of random number rand_it wich is equal on each PE
1531!> and calculate the velocity seed f_n.
1532!> f_n is splitted in vertical direction by the number of PEs in y-direction and
1533!> and each PE calculates a vertical subsection of f_n. At the the end, all
1534!> parts are collected to form the full array.
1535!------------------------------------------------------------------------------!
1536 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
1537
1538
1539    USE indices,                                                               &
1540        ONLY: nx
1541
1542
1543    IMPLICIT NONE
1544
1545    INTEGER(iwp) :: id          !< core ids at respective boundaries
1546    INTEGER(iwp) :: i           !< loop index in x-direction
1547    INTEGER(iwp) :: ii          !< loop index in x-direction
1548    INTEGER(iwp) :: k           !< loop index in z-direction
1549    INTEGER(iwp) :: kk          !< loop index in z-direction
1550    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
1551
1552    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x    !< length scale in x-direction
1553    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
1554    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_y*2
1555    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
1556
1557    REAL(wp) :: nxz_inv         !< inverse of number of grid points in xz-slice
1558    REAL(wp) :: rand_av         !< average of random number
1559    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
1560
1561    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter func in y-dir
1562    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
1563    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
1564    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)     :: f_n     !<  velocity seed
1565    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
1566
1567
1568!
1569!-- Generate random numbers using a seed generated in stg_init.
1570!-- The set of random numbers are modified to have an average of 0 and
1571!-- unit variance.
1572    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:nx+mergp) )
1573
1574    rand_av        = 0.0_wp
1575    rand_sigma_inv = 0.0_wp
1576    nxz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( nx+1 ), KIND=wp )
1577
1578    DO  i = 0, nx
1579       DO  k = nzb, nzt+1
1580          CALL RANDOM_NUMBER( rand_it(k,i) )
1581          rand_av = rand_av + rand_it(k,i)
1582       ENDDO
1583    ENDDO
1584
1585    rand_av = rand_av * nxz_inv
1586
1587    DO  i = 0, nx
1588       DO  k = nzb, nzt+1
1589          rand_it(k,i)   = rand_it(k,i) - rand_av
1590          rand_sigma_inv = rand_sigma_inv + rand_it(k,i) ** 2
1591       ENDDO
1592    ENDDO
1593
1594    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nxz_inv)
1595
1596    DO  i = 0, nx
1597       DO  k = nzb, nzt+1
1598          rand_it(k,i) = rand_it(k,i) * rand_sigma_inv
1599       ENDDO
1600    ENDDO
1601
1602!
1603!-- Periodic fill of random number in space
1604    DO  i = 0, nx
1605       DO  k = 1, mergp
1606          rand_it(nzb-k,i)   = rand_it(nzt+2-k,i)    ! bottom margin
1607          rand_it(nzt+1+k,i) = rand_it(nzb+k-1,i)    ! top margin
1608       ENDDO
1609    ENDDO
1610    DO  i = 1, mergp
1611       DO  k = nzb-mergp, nzt+1+mergp
1612          rand_it(k,-i)   = rand_it(k,nx-i+1)        ! left margin
1613          rand_it(k,nx+i) = rand_it(k,i-1)           ! right margin
1614       ENDDO
1615    ENDDO
1616
1617!
1618!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1619    n_x2 = n_x * 2
1620    n_z2 = n_z * 2
1621    f_n_l  = 0.0_wp
1622
1623    DO  i = nxlg, nxrg
1624       DO  k = nzb_y_stg, nzt_y_stg+1
1625          DO  ii = -n_x2(k), n_x2(k)
1626             DO  kk = -n_z2(k), n_z2(k)
1627                f_n_l(k,i) = f_n_l(k,i)                                        &
1628                           + b_x(ii,k) * b_z(kk,k) * rand_it(k+kk,i+ii)
1629             ENDDO
1630          ENDDO
1631       ENDDO
1632    ENDDO
1633
1634    DEALLOCATE( rand_it )
1635
1636!
1637!-- Gather velocity seeds of full subdomain
1638    send_count = nzt_y_stg - nzb_y_stg + 1
1639    IF ( nzt_y_stg == nzt )  send_count = send_count + 1
1640
1641
1642#if defined( __parallel )
1643    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,    &
1644                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz,  &
1645                      id, comm1dy, ierr )
1646#else
1647    f_n(nzb+1:nzt+1,nxlg:nxrg) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg)
1648#endif
1649
1650
1651 END SUBROUTINE stg_generate_seed_xz
1652
1653 END MODULE synthetic_turbulence_generator_mod
Note: See TracBrowser for help on using the repository browser.