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

Last change on this file since 3318 was 3289, checked in by suehring, 6 years ago

last commit documented

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