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

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

last commit documented

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