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

Last change on this file since 2967 was 2967, checked in by raasch, 3 years ago

bugfix: missing parallel cpp-directives added

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