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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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