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

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

Mask topography while imposing inflow perturbations at the boundaries; do not impose perturbations at top boundary as well as ghost points; Remove print statement; Read zsoil dimension lenght only if soil variables are provided in dynamic driver

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