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

Last change on this file since 3246 was 3246, checked in by sward, 3 years ago

Added error handling for wrong input parameters

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