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

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

various changes to avoid compiler warnings (mainly removal of unused variables)

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