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

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

Bugfix in allocation of mean_inflow_profiles in case of restarts

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