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

Last change on this file since 3248 was 3248, checked in by sward, 6 years ago

Minor format changes

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