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

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

updated variable description

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