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

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

Bugfix for revision 2945, perform checks only if dynamic input file is available; remove unused module load

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