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

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

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

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