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

Last change on this file since 3044 was 3044, checked in by gronemeier, 3 years ago

add missing variable descriptions

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