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

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

Remaining error messages revised, comments extended

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