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

Last change on this file since 3051 was 3051, checked in by suehring, 3 years ago

Speed-up NetCDF input; Revise NetCDF-input routines and remove input via io-blocks; Temporarily revoke renaming of input variables in dynamic driver; More detailed error messages created; Bugfix in mapping 3D buildings; Bugfix in land-surface model at pavement surfaces; Bugfix in initialization with inifor

  • Property svn:keywords set to Id
File size: 62.0 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 3051 2018-05-30 17:43:55Z suehring $
27! Bugfix in calculation of initial Reynolds-stress tensor.
28!
29! 3045 2018-05-28 07:55:41Z Giersch
30! Error messages 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       ELSE
721          a22(k) = 0.0_wp
722       ENDIF
723
724!
725!--    a31, a32, a33 must be calculated with interpolated a11, a21, a22 (d11,
726!--    d21, d22) because of different vertical grid
727       IF ( k .le. nzt )  THEN
728          d11 = 0.5_wp * ( r11(k) + r11(k+1) )
729          IF ( d11 > 0.0_wp )  THEN
730             d11 = SQRT( d11 )
731             d21 = ( 0.5_wp * ( r21(k) + r21(k+1) ) ) / d11
732             a31(k) = r31(k) / d11
733          ELSE
734             d21 = 0.0_wp
735             a31(k) = 0.0_wp
736          ENDIF
737
738          d22 = 0.5_wp * ( r22(k) + r22(k+1) ) - d21 ** 2
739          IF ( d22 > 0.0_wp )  THEN
740             a32(k) = ( r32(k) - d21 * a31(k) ) / SQRT( d22 )
741          ELSE
742             a32(k) = 0.0_wp
743          ENDIF
744
745          a33(k) = r33(k) - a31(k) ** 2 - a32(k) ** 2
746          IF ( a33(k) > 0.0_wp )  THEN
747             a33(k) = SQRT( a33(k) )
748          ELSE
749             a33(k) = 0.0_wp
750          ENDIF
751       ELSE
752          a31(k) = a31(k-1)
753          a32(k) = a32(k-1)
754          a33(k) = a33(k-1)
755       ENDIF
756
757    ENDDO
758!
759!-- Define the size of the filter functions and allocate them.
760    merg = 0
761
762    ! arrays must be large enough to cover the largest length scale
763    DO  k = nzb, nzt+1
764       j = MAX( ABS(nux(k)), ABS(nuy(k)), ABS(nuz(k)), &
765                ABS(nvx(k)), ABS(nvy(k)), ABS(nvz(k)), &
766                ABS(nwx(k)), ABS(nwy(k)), ABS(nwz(k))  )
767       IF ( j > merg )  merg = j
768    ENDDO
769
770    merg  = 2 * merg
771    mergp = merg + nbgp
772
773    ALLOCATE ( bux(-merg:merg,nzb:nzt+1),                                      &
774               buy(-merg:merg,nzb:nzt+1),                                      &
775               buz(-merg:merg,nzb:nzt+1),                                      &
776               bvx(-merg:merg,nzb:nzt+1),                                      &
777               bvy(-merg:merg,nzb:nzt+1),                                      &
778               bvz(-merg:merg,nzb:nzt+1),                                      &
779               bwx(-merg:merg,nzb:nzt+1),                                      &
780               bwy(-merg:merg,nzb:nzt+1),                                      &
781               bwz(-merg:merg,nzb:nzt+1)  )
782
783!
784!-- Allocate velocity seeds for turbulence at xz-layer
785    ALLOCATE ( fu_xz( nzb:nzt+1,nxlg:nxrg), fuo_xz(nzb:nzt+1,nxlg:nxrg),       &
786               fv_xz( nzb:nzt+1,nxlg:nxrg), fvo_xz(nzb:nzt+1,nxlg:nxrg),       &
787               fw_xz( nzb:nzt+1,nxlg:nxrg), fwo_xz(nzb:nzt+1,nxlg:nxrg)  )
788
789!
790!-- Allocate velocity seeds for turbulence at yz-layer
791    ALLOCATE ( fu_yz( nzb:nzt+1,nysg:nyng), fuo_yz(nzb:nzt+1,nysg:nyng),       &
792               fv_yz( nzb:nzt+1,nysg:nyng), fvo_yz(nzb:nzt+1,nysg:nyng),       &
793               fw_yz( nzb:nzt+1,nysg:nyng), fwo_yz(nzb:nzt+1,nysg:nyng)  )
794
795    fu_xz  = 0.0_wp
796    fuo_xz = 0.0_wp
797    fv_xz  = 0.0_wp
798    fvo_xz = 0.0_wp
799    fw_xz  = 0.0_wp
800    fwo_xz = 0.0_wp
801
802    fu_yz  = 0.0_wp
803    fuo_yz = 0.0_wp
804    fv_yz  = 0.0_wp
805    fvo_yz = 0.0_wp
806    fw_yz  = 0.0_wp
807    fwo_yz = 0.0_wp
808
809!
810!-- Create filter functions
811    CALL stg_filter_func( nux, bux ) !filter ux
812    CALL stg_filter_func( nuy, buy ) !filter uy
813    CALL stg_filter_func( nuz, buz ) !filter uz
814    CALL stg_filter_func( nvx, bvx ) !filter vx
815    CALL stg_filter_func( nvy, bvy ) !filter vy
816    CALL stg_filter_func( nvz, bvz ) !filter vz
817    CALL stg_filter_func( nwx, bwx ) !filter wx
818    CALL stg_filter_func( nwy, bwy ) !filter wy
819    CALL stg_filter_func( nwz, bwz ) !filter wz
820
821#if defined( __parallel )
822    CALL MPI_BARRIER( comm2d, ierr )
823#endif
824
825!
826!-- In case of restart, calculate velocity seeds fu, fv, fw from former
827!   time step.
828!   Bug: fu, fv, fw are different in those heights where a11, a22, a33
829!        are 0 compared to the prerun. This is mostly for k=nzt+1.
830    IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
831       IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
832
833          IF ( myidx == id_stg_left  )  i = -1
834          IF ( myidx == id_stg_right )  i = nxr+1
835
836          DO  j = nysg, nyng
837             DO  k = nzb, nzt+1
838
839                IF  ( a11(k) .NE. 0._wp ) THEN
840                   fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
841                ELSE
842                   fu_yz(k,j) = 0._wp
843                ENDIF
844
845                IF  ( a22(k) .NE. 0._wp ) THEN
846                   fv_yz(k,j) = ( v(k,j,i) / mc_factor - a21(k) * fu_yz(k,j) - &
847                               v_init(k) ) / a22(k)
848                ELSE
849                   fv_yz(k,j) = 0._wp
850                ENDIF
851
852                IF  ( a33(k) .NE. 0._wp ) THEN
853                   fw_yz(k,j) = ( w(k,j,i) / mc_factor - a31(k) * fu_yz(k,j) - &
854                               a32(k) * fv_yz(k,j) ) / a33(k)
855                ELSE
856                   fw_yz = 0._wp
857                ENDIF
858
859             ENDDO
860          ENDDO
861       ENDIF
862    ENDIF
863
864    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
865
866 END SUBROUTINE stg_init
867
868
869!------------------------------------------------------------------------------!
870! Description:
871! ------------
872!> Calculate filter function bxx from length scale nxx following Eg.9 and 10
873!> (Xie and Castro, 2008)
874!------------------------------------------------------------------------------!
875 SUBROUTINE stg_filter_func( nxx, bxx )
876
877
878    IMPLICIT NONE
879
880    INTEGER(iwp) :: k         !< loop index
881    INTEGER(iwp) :: n_k       !< length scale nXX in height k
882    INTEGER(iwp) :: n_k2      !< n_k * 2
883    INTEGER(iwp) :: nf        !< index for length scales
884
885    REAL(wp) :: bdenom        !< denominator for filter functions bXX
886    REAL(wp) :: qsi = 1.0_wp  !< minimization factor
887
888    INTEGER(iwp), DIMENSION(:) :: nxx(nzb:nzt+1)           !< length scale (in gp)
889
890    REAL(wp), DIMENSION(:,:) :: bxx(-merg:merg,nzb:nzt+1)  !< filter function
891
892
893    bxx = 0.0_wp
894
895    DO  k = nzb, nzt+1
896       bdenom = 0.0_wp
897       n_k    = nxx(k)
898       IF ( n_k /= 0 )  THEN
899          n_k2 = n_k * 2
900
901!
902!--       ( Eq.10 )^2
903          DO  nf = -n_k2, n_k2
904             bdenom = bdenom + EXP( -qsi * pi * ABS(nf) / n_k )**2
905          ENDDO
906
907!
908!--       ( Eq.9 )
909          bdenom = SQRT( bdenom )
910          DO  nf = -n_k2, n_k2
911             bxx(nf,k) = EXP( -qsi * pi * ABS(nf) / n_k ) / bdenom
912          ENDDO
913       ENDIF
914    ENDDO
915
916END SUBROUTINE stg_filter_func
917
918
919!------------------------------------------------------------------------------!
920! Description:
921! ------------
922!> Parin for &stg_par for synthetic turbulence generator
923!------------------------------------------------------------------------------!
924 SUBROUTINE stg_parin
925
926
927    IMPLICIT NONE
928
929    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
930
931
932    NAMELIST /stg_par/   use_syn_turb_gen
933
934    line = ' '
935
936!
937!-- Try to find stg package
938    REWIND ( 11 )
939    line = ' '
940    DO WHILE ( INDEX( line, '&stg_par' ) == 0 )
941       READ ( 11, '(A)', END=10 )  line
942    ENDDO
943    BACKSPACE ( 11 )
944
945!
946!-- Read namelist
947    READ ( 11, stg_par )
948
949!
950!-- Set flag that indicates that the synthetic turbulence generator is switched
951!-- on
952    syn_turb_gen = .TRUE.
953
954
955 10 CONTINUE
956
957 END SUBROUTINE stg_parin
958
959
960!------------------------------------------------------------------------------!
961! Description:
962! ------------
963!> This routine reads the respective restart data.
964!------------------------------------------------------------------------------!
965 SUBROUTINE stg_rrd_global( found )
966
967
968    USE control_parameters,                                                    &
969        ONLY: length, restart_string
970
971
972    IMPLICIT NONE
973
974    LOGICAL, INTENT(OUT)  ::  found !< flag indicating if variable was found
975
976
977    found = .TRUE.
978
979
980    SELECT CASE ( restart_string(1:length) )
981
982       CASE ( 'mc_factor' )
983          READ ( 13 )  mc_factor
984       CASE ( 'use_syn_turb_gen' )
985          READ ( 13 )  use_syn_turb_gen
986
987       CASE DEFAULT
988
989          found = .FALSE.
990
991    END SELECT
992
993
994 END SUBROUTINE stg_rrd_global
995
996
997!------------------------------------------------------------------------------!
998! Description:
999! ------------
1000!> This routine writes the respective restart data.
1001!------------------------------------------------------------------------------!
1002 SUBROUTINE stg_wrd_global
1003
1004
1005    IMPLICIT NONE
1006
1007    CALL wrd_write_string( 'mc_factor' )
1008    WRITE ( 14 )  mc_factor
1009
1010    CALL wrd_write_string( 'use_syn_turb_gen' )
1011    WRITE ( 14 )  use_syn_turb_gen
1012
1013
1014 END SUBROUTINE stg_wrd_global
1015
1016
1017!------------------------------------------------------------------------------!
1018! Description:
1019! ------------
1020!> Create turbulent inflow fields for u, v, w with prescribed length scales and
1021!> Reynolds stress tensor after a method of Xie and Castro (2008), modified
1022!> following suggestions of Kim et al. (2013), and using a Lund rotation
1023!> (Lund, 1998).
1024!------------------------------------------------------------------------------!
1025 SUBROUTINE stg_main
1026
1027
1028    USE arrays_3d,                                                             &
1029        ONLY:  dzw
1030
1031    USE control_parameters,                                                    &
1032        ONLY:  dt_3d, forcing, intermediate_timestep_count,  nest_domain,      &
1033               rans_mode, simulated_time, volume_flow_initial
1034
1035    USE grid_variables,                                                        &
1036        ONLY:  dx, dy
1037
1038    USE indices,                                                               &
1039        ONLY:  wall_flags_0
1040
1041    USE statistics,                                                            &
1042        ONLY:  weight_substep
1043
1044    USE pmc_interface,                                                         &
1045        ONLY : rans_mode_parent
1046
1047
1048    IMPLICIT NONE
1049
1050    INTEGER(iwp) :: i           !< grid index in x-direction
1051    INTEGER(iwp) :: j           !< loop index in y-direction
1052    INTEGER(iwp) :: k           !< loop index in z-direction
1053
1054    REAL(wp) :: dt_stg          !< wheighted subtimestep
1055    REAL(wp) :: mc_factor_l     !< local mass flux correction factor
1056    REAL(wp) :: volume_flow     !< mass flux through lateral boundary
1057    REAL(wp) :: volume_flow_l   !< local mass flux through lateral boundary
1058
1059    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg,5) :: dist_xz !< imposed disturbances at north/south boundary
1060    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5) :: dist_yz !< imposed disturbances at left/right boundary
1061
1062
1063    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1064
1065!
1066!-- Calculate time step which is needed for filter functions
1067    dt_stg = dt_3d * weight_substep(intermediate_timestep_count)
1068
1069!
1070!-- Initial value of fu, fv, fw
1071    IF ( simulated_time == 0.0_wp .AND. .NOT. velocity_seed_initialized )  THEN
1072       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
1073       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
1074       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
1075
1076       IF ( forcing  .OR.  ( nest_domain .AND.  rans_mode_parent  .AND.        &
1077                      .NOT.  rans_mode ) )  THEN
1078!
1079!--       Generate turbulence at right boundary
1080          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_right )
1081          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_right )
1082          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_right )
1083!
1084!--       Generate turbulence at north boundary
1085          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_north )
1086          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_north )
1087          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_north )
1088!
1089!--       Generate turbulence at south boundary
1090          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_south )
1091          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_south )
1092          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_south )
1093       ENDIF
1094       velocity_seed_initialized = .TRUE.
1095    ENDIF
1096!
1097!-- New set of fu, fv, fw
1098    CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
1099    CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
1100    CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
1101
1102    IF ( forcing  .OR.  ( nest_domain .AND.  rans_mode_parent  .AND.           &
1103                   .NOT.  rans_mode ) )  THEN
1104!
1105!--       Generate turbulence at right boundary
1106          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
1107          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
1108          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
1109!
1110!--       Generate turbulence at north boundary
1111          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
1112          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
1113          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
1114!
1115!--       Generate turbulence at south boundary
1116          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
1117          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
1118          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
1119    ENDIF
1120!
1121!-- Turbulence generation at left and or right boundary
1122    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
1123
1124       DO  j = nysg, nyng
1125          DO  k = nzb, nzt + 1
1126!
1127!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1128             IF ( tu(k) == 0.0_wp )  THEN
1129                fu_yz(k,j) = fuo_yz(k,j)
1130             ELSE
1131                fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
1132                         fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
1133             ENDIF
1134
1135             IF ( tv(k) == 0.0_wp )  THEN
1136                fv_yz(k,j) = fvo_yz(k,j)
1137             ELSE
1138                fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
1139                         fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
1140             ENDIF
1141
1142             IF ( tw(k) == 0.0_wp )  THEN
1143                fw_yz(k,j) = fwo_yz(k,j)
1144             ELSE
1145                fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
1146                         fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
1147             ENDIF
1148!
1149!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
1150!--          Additional factors are added to improve the variance of v and w
1151             IF( k == 0 )  THEN
1152                dist_yz(k,j,1) = 0.0_wp
1153                dist_yz(k,j,2) = 0.0_wp
1154                dist_yz(k,j,3) = 0.0_wp
1155!                 dist_yz(k,j,4) = 0.0_wp
1156!                 dist_yz(k,j,5) = 0.0_wp
1157             ELSE
1158                dist_yz(k,j,1) = a11(k) * fu_yz(k,j)
1159                !experimental test of 1.2
1160                dist_yz(k,j,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
1161                                         * 1.2_wp )                            &
1162                                       * (   a21(k) * fu_yz(k,j)               &
1163                                           + a22(k) * fv_yz(k,j) )
1164                dist_yz(k,j,3) = ( SQRT(a33(k) / MAXVAL(a33) )                 &
1165                                         * 1.3_wp )                            &
1166                                       * (   a31(k) * fu_yz(k,j)               &
1167                                           + a32(k) * fv_yz(k,j)               &
1168                                           + a33(k) * fw_yz(k,j) )
1169                ! Calculation for pt and e not yet implemented
1170!                 dist_yz(k,j,4) = 0.0_wp
1171!                 dist_yz(k,j,5) = 0.0_wp
1172             ENDIF
1173
1174          ENDDO
1175       ENDDO
1176
1177!
1178!--    Mass flux correction following Kim et al. (2013)
1179!--    This correction factor insures that the mass flux is preserved at the
1180!--    inflow boundary
1181       IF ( .NOT. forcing  .AND.  .NOT. nest_domain )  THEN
1182          mc_factor_l = 0.0_wp
1183          mc_factor   = 0.0_wp
1184          DO  j = nys, nyn
1185             DO  k = nzb+1, nzt
1186                mc_factor_l = mc_factor_l + dzw(k)  *                          &
1187                              ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) )
1188             ENDDO
1189          ENDDO
1190
1191#if defined( __parallel )
1192          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,  &
1193                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1194#else
1195          mc_factor = mc_factor_l
1196#endif
1197
1198          mc_factor = volume_flow_initial(1) / mc_factor
1199
1200!
1201!--       Add disturbance at the inflow
1202          DO  j = nysg, nyng
1203             DO  k = nzb, nzt+1
1204                 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +              &
1205                                      dist_yz(k,j,1)             ) * mc_factor
1206                 v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +              &
1207                                      dist_yz(k,j,2)             ) * mc_factor
1208                 w(k,j,-nbgp:-1)  =   dist_yz(k,j,3)               * mc_factor
1209             ENDDO
1210          ENDDO
1211
1212       ELSE
1213!
1214!--       First, calculate volume flow at yz boundary
1215          IF ( myidx == id_stg_left  )  i = nxl
1216          IF ( myidx == id_stg_right )  i = nxr+1
1217
1218          volume_flow_l = 0.0_wp
1219          volume_flow   = 0.0_wp
1220          mc_factor_l   = 0.0_wp
1221          mc_factor     = 0.0_wp
1222          DO  j = nys, nyn
1223             DO  k = nzb+1, nzt
1224                volume_flow_l = volume_flow_l + u(k,j,i) * dzw(k) * dy         &
1225                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1226                                              BTEST( wall_flags_0(k,j,i), 1 ) )
1227
1228                mc_factor_l = mc_factor_l     + ( u(k,j,i) + dist_yz(k,j,1) )  &
1229                                                         * dzw(k) * dy         &
1230                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1231                                              BTEST( wall_flags_0(k,j,i), 1 ) )
1232             ENDDO
1233          ENDDO
1234#if defined( __parallel )
1235          CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                      &
1236                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1237          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                          &
1238                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1239#else
1240          volume_flow = volume_flow_l
1241          mc_factor   = mc_factor_l
1242#endif
1243
1244          mc_factor = volume_flow / mc_factor
1245
1246!
1247!--       Add disturbances
1248          IF ( myidx == id_stg_left  )  THEN
1249
1250             DO  j = nysg, nyng
1251                DO  k = nzb, nzt+1
1252                   u(k,j,-nbgp+1:0) = ( u(k,j,-nbgp+1:0) + dist_yz(k,j,1) )    &
1253                                        * mc_factor
1254                   v(k,j,-nbgp:-1)  = ( v(k,j,-nbgp:-1)  + dist_yz(k,j,2)  )   &
1255                                        * mc_factor
1256                   w(k,j,-nbgp:-1)  = ( w(k,j,-nbgp:-1)  + dist_yz(k,j,3)  )   &
1257                                        * mc_factor
1258                ENDDO
1259             ENDDO
1260          ENDIF
1261          IF ( myidx == id_stg_right  )  THEN
1262
1263             DO  j = nysg, nyng
1264                DO  k = nzb, nzt+1
1265                   u(k,j,nxr+1:nxr+nbgp) = ( u(k,j,nxr+1:nxr+nbgp) +           &
1266                                             dist_yz(k,j,1) ) * mc_factor
1267                   v(k,j,nxr+1:nxr+nbgp) = ( v(k,j,nxr+1:nxr+nbgp) +           &
1268                                             dist_yz(k,j,2) ) * mc_factor
1269                   w(k,j,nxr+1:nxr+nbgp) = ( w(k,j,nxr+1:nxr+nbgp) +           &
1270                                             dist_yz(k,j,3) ) * mc_factor
1271                ENDDO
1272             ENDDO
1273          ENDIF
1274       ENDIF
1275
1276    ENDIF
1277!
1278!-- Turbulence generation at north and south boundary
1279    IF ( myidy == id_stg_north  .OR.  myidy == id_stg_south )  THEN
1280
1281       DO  i = nxlg, nxrg
1282          DO  k = nzb, nzt + 1
1283!
1284!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1285             IF ( tu(k) == 0.0_wp )  THEN
1286                fu_xz(k,i) = fuo_xz(k,i)
1287             ELSE
1288                fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
1289                         fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
1290             ENDIF
1291
1292             IF ( tv(k) == 0.0_wp )  THEN
1293                fv_xz(k,i) = fvo_xz(k,i)
1294             ELSE
1295                fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
1296                         fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
1297             ENDIF
1298
1299             IF ( tw(k) == 0.0_wp )  THEN
1300                fw_xz(k,i) = fwo_xz(k,i)
1301             ELSE
1302                fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
1303                         fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
1304             ENDIF
1305!
1306!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
1307!--          Additional factors are added to improve the variance of v and w
1308             IF( k == 0 )  THEN
1309                dist_xz(k,i,1) = 0.0_wp
1310                dist_xz(k,i,2) = 0.0_wp
1311                dist_xz(k,i,3) = 0.0_wp
1312
1313             ELSE
1314                dist_xz(k,i,1) = a11(k) * fu_xz(k,i)
1315                !experimental test of 1.2
1316                dist_xz(k,i,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
1317                                         * 1.2_wp )                            &
1318                                       * (   a21(k) * fu_xz(k,i)               &
1319                                           + a22(k) * fv_xz(k,i) )
1320                dist_xz(k,i,3) = ( SQRT(a33(k) / MAXVAL(a33) )                 &
1321                                         * 1.3_wp )                            &
1322                                       * (   a31(k) * fu_xz(k,i)               &
1323                                           + a32(k) * fv_xz(k,i)               &
1324                                           + a33(k) * fw_xz(k,i) )
1325             ENDIF
1326
1327          ENDDO
1328       ENDDO
1329!
1330!--    Mass flux correction following Kim et al. (2013)
1331!--    This correction factor insures that the mass flux is preserved at the
1332!--    inflow boundary.
1333!--    First, calculate volume flow at xz boundary
1334       IF ( myidy == id_stg_south  ) j = nys
1335       IF ( myidy == id_stg_north )  j = nyn+1
1336
1337       volume_flow_l = 0.0_wp
1338       volume_flow   = 0.0_wp
1339       mc_factor_l   = 0.0_wp
1340       mc_factor     = 0.0_wp
1341       DO  i = nxl, nxr
1342          DO  k = nzb+1, nzt
1343             volume_flow_l = volume_flow_l + v(k,j,i) * dzw(k) * dx            &
1344                                  * MERGE( 1.0_wp, 0.0_wp,                     &
1345                                           BTEST( wall_flags_0(k,j,i), 2 ) )
1346
1347             mc_factor_l = mc_factor_l     + ( v(k,j,i) + dist_xz(k,i,2) )     &
1348                                                      * dzw(k) * dx            &
1349                                  * MERGE( 1.0_wp, 0.0_wp,                     &
1350                                           BTEST( wall_flags_0(k,j,i), 2 ) )
1351          ENDDO
1352       ENDDO
1353#if defined( __parallel )
1354       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                         &
1355                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1356       CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                             &
1357                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1358#else
1359       volume_flow = volume_flow_l
1360       mc_factor   = mc_factor_l
1361#endif
1362
1363       mc_factor = volume_flow / mc_factor
1364
1365!
1366!--    Add disturbances
1367       IF ( myidy == id_stg_south  )  THEN
1368
1369          DO  i = nxlg, nxrg
1370             DO  k = nzb, nzt+1
1371                u(k,-nbgp:-1,i) = ( u(k,-nbgp:-1,i) + dist_xz(k,i,1) )         &
1372                                     * mc_factor
1373                v(k,-nbgp:0,i)  = ( v(k,-nbgp:0,i)  + dist_xz(k,i,2)  )        &
1374                                     * mc_factor
1375                w(k,-nbgp:-1,i) = ( w(k,-nbgp:-1,i) + dist_xz(k,i,3)  )        &
1376                                     * mc_factor
1377             ENDDO
1378          ENDDO
1379       ENDIF
1380       IF ( myidy == id_stg_north  )  THEN
1381
1382          DO  i = nxlg, nxrg
1383             DO  k = nzb, nzt+1
1384                u(k,nyn+1:nyn+nbgp,i) = ( u(k,nyn+1:nyn+nbgp,i) +              &
1385                                          dist_xz(k,i,1) ) * mc_factor
1386                v(k,nyn+1:nyn+nbgp,i) = ( v(k,nyn+1:nyn+nbgp,i) +              &
1387                                          dist_xz(k,i,2) ) * mc_factor
1388                w(k,nyn+1:nyn+nbgp,i) = ( w(k,nyn+1:nyn+nbgp,i) +              &
1389                                          dist_xz(k,i,3) ) * mc_factor
1390             ENDDO
1391          ENDDO
1392       ENDIF
1393    ENDIF
1394
1395    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1396
1397 END SUBROUTINE stg_main
1398
1399!------------------------------------------------------------------------------!
1400! Description:
1401! ------------
1402!> Generate a set of random number rand_it wich is equal on each PE
1403!> and calculate the velocity seed f_n.
1404!> f_n is splitted in vertical direction by the number of PEs in x-direction and
1405!> and each PE calculates a vertical subsection of f_n. At the the end, all
1406!> parts are collected to form the full array.
1407!------------------------------------------------------------------------------!
1408 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
1409
1410
1411    USE indices,                                                               &
1412        ONLY: ny
1413
1414    IMPLICIT NONE
1415
1416    INTEGER(iwp) :: id          !< core ids at respective boundaries
1417    INTEGER(iwp) :: j           !< loop index in y-direction
1418    INTEGER(iwp) :: jj          !< loop index in y-direction
1419    INTEGER(iwp) :: k           !< loop index in z-direction
1420    INTEGER(iwp) :: kk          !< loop index in z-direction
1421    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
1422
1423    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y    !< length scale in y-direction
1424    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
1425    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y2   !< n_y*2
1426    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
1427
1428    REAL(wp) :: nyz_inv         !< inverse of number of grid points in yz-slice
1429    REAL(wp) :: rand_av         !< average of random number
1430    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
1431
1432    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter func in y-dir
1433    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
1434    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
1435    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)     :: f_n     !<  velocity seed
1436    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
1437
1438
1439!
1440!-- Generate random numbers using a seed generated in stg_init.
1441!-- The set of random numbers are modified to have an average of 0 and
1442!-- unit variance.
1443    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:ny+mergp) )
1444
1445    rand_av        = 0.0_wp
1446    rand_sigma_inv = 0.0_wp
1447    nyz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( ny+1 ), KIND=wp )
1448
1449    DO  j = 0, ny
1450       DO  k = nzb, nzt+1
1451          CALL RANDOM_NUMBER( rand_it(k,j) )
1452          rand_av = rand_av + rand_it(k,j)
1453       ENDDO
1454    ENDDO
1455
1456    rand_av = rand_av * nyz_inv
1457
1458    DO  j = 0, ny
1459       DO  k = nzb, nzt+1
1460          rand_it(k,j)   = rand_it(k,j) - rand_av
1461          rand_sigma_inv = rand_sigma_inv + rand_it(k,j) ** 2
1462       ENDDO
1463    ENDDO
1464
1465    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nyz_inv)
1466
1467    DO  j = 0, ny
1468       DO  k = nzb, nzt+1
1469          rand_it(k,j) = rand_it(k,j) * rand_sigma_inv
1470       ENDDO
1471    ENDDO
1472
1473!
1474!-- Periodic fill of random number in space
1475    DO  j = 0, ny
1476       DO  k = 1, mergp
1477          rand_it(nzb  -k,j) = rand_it(nzt+2-k,j)    ! bottom margin
1478          rand_it(nzt+1+k,j) = rand_it(nzb+k-1,j)    ! top margin
1479       ENDDO
1480    ENDDO
1481    DO  j = 1, mergp
1482       DO  k = nzb-mergp, nzt+1+mergp
1483          rand_it(k,  -j) = rand_it(k,ny-j+1)        ! south margin
1484          rand_it(k,ny+j) = rand_it(k,   j-1)        ! north margin
1485       ENDDO
1486    ENDDO
1487
1488!
1489!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1490    n_y2 = n_y * 2
1491    n_z2 = n_z * 2
1492    f_n_l  = 0.0_wp
1493
1494    DO  j = nysg, nyng
1495       DO  k = nzb_x_stg, nzt_x_stg+1
1496          DO  jj = -n_y2(k), n_y2(k)
1497             DO  kk = -n_z2(k), n_z2(k)
1498                f_n_l(k,j) = f_n_l(k,j)                                        &
1499                           + b_y(jj,k) * b_z(kk,k) * rand_it(k+kk,j+jj)
1500             ENDDO
1501          ENDDO
1502       ENDDO
1503    ENDDO
1504
1505    DEALLOCATE( rand_it )
1506!
1507!-- Gather velocity seeds of full subdomain
1508    send_count = nzt_x_stg - nzb_x_stg + 1
1509    IF ( nzt_x_stg == nzt )  send_count = send_count + 1
1510
1511#if defined( __parallel )
1512    CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,    &
1513                      f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,  &
1514                      id, comm1dx, ierr )
1515#else
1516    f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x_stg:nzt_x_stg+1,nysg:nyng)
1517#endif
1518
1519
1520 END SUBROUTINE stg_generate_seed_yz
1521
1522
1523
1524
1525!------------------------------------------------------------------------------!
1526! Description:
1527! ------------
1528!> Generate a set of random number rand_it wich is equal on each PE
1529!> and calculate the velocity seed f_n.
1530!> f_n is splitted in vertical direction by the number of PEs in y-direction and
1531!> and each PE calculates a vertical subsection of f_n. At the the end, all
1532!> parts are collected to form the full array.
1533!------------------------------------------------------------------------------!
1534 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
1535
1536
1537    USE indices,                                                               &
1538        ONLY: nx
1539
1540
1541    IMPLICIT NONE
1542
1543    INTEGER(iwp) :: id          !< core ids at respective boundaries
1544    INTEGER(iwp) :: i           !< loop index in x-direction
1545    INTEGER(iwp) :: ii          !< loop index in x-direction
1546    INTEGER(iwp) :: k           !< loop index in z-direction
1547    INTEGER(iwp) :: kk          !< loop index in z-direction
1548    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
1549
1550    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x    !< length scale in x-direction
1551    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
1552    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_y*2
1553    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
1554
1555    REAL(wp) :: nxz_inv         !< inverse of number of grid points in xz-slice
1556    REAL(wp) :: rand_av         !< average of random number
1557    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
1558
1559    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter func in y-dir
1560    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
1561    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
1562    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)     :: f_n     !<  velocity seed
1563    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
1564
1565
1566!
1567!-- Generate random numbers using a seed generated in stg_init.
1568!-- The set of random numbers are modified to have an average of 0 and
1569!-- unit variance.
1570    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:nx+mergp) )
1571
1572    rand_av        = 0.0_wp
1573    rand_sigma_inv = 0.0_wp
1574    nxz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( nx+1 ), KIND=wp )
1575
1576    DO  i = 0, nx
1577       DO  k = nzb, nzt+1
1578          CALL RANDOM_NUMBER( rand_it(k,i) )
1579          rand_av = rand_av + rand_it(k,i)
1580       ENDDO
1581    ENDDO
1582
1583    rand_av = rand_av * nxz_inv
1584
1585    DO  i = 0, nx
1586       DO  k = nzb, nzt+1
1587          rand_it(k,i)   = rand_it(k,i) - rand_av
1588          rand_sigma_inv = rand_sigma_inv + rand_it(k,i) ** 2
1589       ENDDO
1590    ENDDO
1591
1592    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nxz_inv)
1593
1594    DO  i = 0, nx
1595       DO  k = nzb, nzt+1
1596          rand_it(k,i) = rand_it(k,i) * rand_sigma_inv
1597       ENDDO
1598    ENDDO
1599
1600!
1601!-- Periodic fill of random number in space
1602    DO  i = 0, nx
1603       DO  k = 1, mergp
1604          rand_it(nzb-k,i)   = rand_it(nzt+2-k,i)    ! bottom margin
1605          rand_it(nzt+1+k,i) = rand_it(nzb+k-1,i)    ! top margin
1606       ENDDO
1607    ENDDO
1608    DO  i = 1, mergp
1609       DO  k = nzb-mergp, nzt+1+mergp
1610          rand_it(k,-i)   = rand_it(k,nx-i+1)        ! left margin
1611          rand_it(k,nx+i) = rand_it(k,i-1)           ! right margin
1612       ENDDO
1613    ENDDO
1614
1615!
1616!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1617    n_x2 = n_x * 2
1618    n_z2 = n_z * 2
1619    f_n_l  = 0.0_wp
1620
1621    DO  i = nxlg, nxrg
1622       DO  k = nzb_y_stg, nzt_y_stg+1
1623          DO  ii = -n_x2(k), n_x2(k)
1624             DO  kk = -n_z2(k), n_z2(k)
1625                f_n_l(k,i) = f_n_l(k,i)                                        &
1626                           + b_x(ii,k) * b_z(kk,k) * rand_it(k+kk,i+ii)
1627             ENDDO
1628          ENDDO
1629       ENDDO
1630    ENDDO
1631
1632    DEALLOCATE( rand_it )
1633
1634!
1635!-- Gather velocity seeds of full subdomain
1636    send_count = nzt_y_stg - nzb_y_stg + 1
1637    IF ( nzt_y_stg == nzt )  send_count = send_count + 1
1638
1639
1640#if defined( __parallel )
1641    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,    &
1642                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz,  &
1643                      id, comm1dy, ierr )
1644#else
1645    f_n(nzb+1:nzt+1,nxlg:nxrg) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg)
1646#endif
1647
1648
1649 END SUBROUTINE stg_generate_seed_xz
1650
1651 END MODULE synthetic_turbulence_generator_mod
Note: See TracBrowser for help on using the repository browser.