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

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

Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode

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