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

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

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

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