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

Last change on this file since 3065 was 3065, checked in by Giersch, 3 years ago

New vertical stretching procedure has been introduced

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