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

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

Bugfix in parallelization of synthetic turbulence generator; revision in control mimic of synthetic turbulence generator in case of RAN-LES nesting; checks for dynamic input file added; control mimic for topography input slightly revised.

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