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

Last change on this file since 3371 was 3349, checked in by suehring, 6 years ago

Bugfix in initialization of soil properties from dynamic input file

  • Property svn:keywords set to Id
File size: 81.4 KB
Line 
1!> @synthetic_turbulence_generator_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: synthetic_turbulence_generator_mod.f90 3349 2018-10-15 16:39:41Z knoop $
27! Fix for format descriptor
28!
29! 3348 2018-10-15 14:30:51Z suehring
30! - Revise structure of init routine
31! - introduce new parameters to skip STG for some timesteps
32! - introduce time-dependent parametrization of Reynolds-stress tensor
33! - Bugfix in allocation of mean_inflow_profiles
34!
35! 3341 2018-10-15 10:31:27Z suehring
36! Introduce module parameter for number of inflow profiles
37!
38! 3288 2018-09-28 10:23:08Z suehring
39! Modularization of all bulk cloud physics code components
40!
41! 3248 2018-09-14 09:42:06Z sward
42! Minor formating changes
43!
44! 3246 2018-09-13 15:14:50Z sward
45! Added error handling for input namelist via parin_fail_message
46!
47! 3241 2018-09-12 15:02:00Z raasch
48! unused variables removed
49!
50! 3186 2018-07-30 17:07:14Z suehring
51! Mask topography while imposing inflow perturbations at the boundaries; do not
52! impose perturbations at top boundary as well as ghost points
53!
54! 3183 2018-07-27 14:25:55Z suehring
55! Rename variables and extend error message
56! Enable geneartor also for stretched grids
57!
58! 3182 2018-07-27 13:36:03Z suehring
59! Error message related to vertical stretching has been added, dz was replaced
60! by dz(1)
61!
62! 3051 2018-05-30 17:43:55Z suehring
63! Bugfix in calculation of initial Reynolds-stress tensor.
64!
65! 3045 2018-05-28 07:55:41Z Giersch
66! Error messages revised
67!
68! 3044 2018-05-25 10:59:41Z gronemeier
69! Add missing variable descriptions
70!
71! 3038 2018-05-24 10:54:00Z gronemeier
72! updated variable descriptions
73!
74! 2967 2018-04-13 11:22:08Z raasch
75! bugfix: missing parallel cpp-directives added
76!
77! 2946 2018-04-04 17:01:23Z suehring
78! Remove unused module load
79!
80! 2945 2018-04-04 16:27:14Z suehring
81! - Bugfix in parallelization of synthetic turbulence generator in case the
82!   z-dimension is not an integral divisor of the number of processors along
83!   the x- and y-dimension
84! - Revision in control mimic in case of RAN-LES nesting
85!
86! 2938 2018-03-27 15:52:42Z suehring
87! Apply turbulence generator at all non-cyclic lateral boundaries in case of
88! realistic Inifor large-scale forcing or RANS-LES nesting
89!
90! 2936 2018-03-27 14:49:27Z suehring
91! variable named found has been introduced for checking if restart data was found,
92! reading of restart strings has been moved completely to read_restart_data_mod,
93! redundant skipping function has been removed, stg_read/write_restart_data
94! have been renamed to stg_r/wrd_global, stg_rrd_global is called in
95! read_restart_data_mod now, flag syn_turb_gen_prerun and marker *** end stg
96! *** have been removed (Giersch), strings and their respective lengths are
97! written out and read now in case of restart runs to get rid of prescribed
98! character lengths (Giersch), CASE DEFAULT was added if restart data is read
99!
100! 2841 2018-02-27 15:02:57Z suehring
101! Bugfix: wrong placement of include 'mpif.h' corrected
102!
103! 2836 2018-02-26 13:40:05Z Giersch
104! The variables synthetic_turbulence_generator and
105! use_synthetic_turbulence_generator have been abbreviated + syn_turb_gen_prerun
106! flag is used to define if module related parameters were outputted as restart
107! data
108!
109! 2716 2017-12-29 16:35:59Z kanani
110! Corrected "Former revisions" section
111!
112! 2696 2017-12-14 17:12:51Z kanani
113! Change in file header (GPL part)
114!
115! 2669 2017-12-06 16:03:27Z raasch
116! unit number for file containing turbulence generator data changed to 90
117! bugfix: preprocessor directives added for MPI specific code
118!
119! 2576 2017-10-24 13:49:46Z Giersch
120! Definition of a new function called stg_skip_global to skip module
121! parameters during reading restart data
122!
123! 2563 2017-10-19 15:36:10Z Giersch
124! stg_read_restart_data is called in stg_parin in the case of a restart run
125!
126! 2259 2017-06-08 09:09:11Z gronemeier
127! Initial revision
128!
129!
130!
131! Authors:
132! --------
133! @author Tobias Gronemeier, Matthias Suehring, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
134!
135!
136!
137! Description:
138! ------------
139!> The module generates turbulence at the inflow boundary based on a method by
140!> Xie and Castro (2008) utilizing a Lund rotation (Lund, 1998) and a mass-flux
141!> correction by Kim et al. (2013).
142!> The turbulence is correlated based on length scales in y- and z-direction and
143!> a time scale for each velocity component. The profiles of length and time
144!> scales, mean u, v, w, e and pt, and all components of the Reynolds stress
145!> tensor can be either read from file STG_PROFILES, or will be parametrized
146!> within the boundary layer.
147!>
148!> @todo test restart
149!>       enable cyclic_fill
150!>       implement turbulence generation for e and pt
151!> @todo Input of height-constant length scales via namelist
152!> @note <Enter notes on the module>
153!> @bug  Height information from input file is not used. Profiles from input
154!>       must match with current PALM grid.
155!>       In case of restart, velocity seeds differ from precursor run if a11,
156!>       a22, or a33 are zero.
157!------------------------------------------------------------------------------!
158 MODULE synthetic_turbulence_generator_mod
159
160
161    USE arrays_3d,                                                             &
162        ONLY:  mean_inflow_profiles, q, pt, u, v, w, zu, zw
163
164    USE basic_constants_and_equations_mod,                                     &
165        ONLY:  g, kappa, pi
166
167    USE control_parameters,                                                    &
168        ONLY:  initializing_actions, num_mean_inflow_profiles, message_string, &
169               syn_turb_gen
170
171    USE cpulog,                                                                &
172        ONLY:  cpu_log, log_point
173
174    USE indices,                                                               &
175        ONLY:  nbgp, nzb, nzt, nxl, nxlg, nxr, nxrg, nys, nyn, nyng, nysg
176
177    USE kinds
178
179#if defined( __parallel )  &&  !defined( __mpifh )
180    USE MPI
181#endif
182
183    USE nesting_offl_mod,                                                      &
184        ONLY:  zi_ribulk
185
186    USE pegrid,                                                                &
187        ONLY:  comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
188
189    USE transpose_indices,                                                     &
190        ONLY: nzb_x, nzt_x
191
192
193    IMPLICIT NONE
194
195#if defined( __parallel )  &&  defined( __mpifh )
196    INCLUDE "mpif.h"
197#endif
198
199
200    LOGICAL ::  velocity_seed_initialized = .FALSE.     !< true after first call of stg_main
201    LOGICAL ::  parametrize_inflow_turbulence = .FALSE. !< flag indicating that inflow turbulence is either read from file (.FALSE.) or if it parametrized
202    LOGICAL ::  use_syn_turb_gen = .FALSE.              !< switch to use synthetic turbulence generator
203
204    INTEGER(iwp) ::  id_stg_left        !< left lateral boundary core id in case of turbulence generator
205    INTEGER(iwp) ::  id_stg_north       !< north lateral boundary core id in case of turbulence generator
206    INTEGER(iwp) ::  id_stg_right       !< right lateral boundary core id in case of turbulence generator
207    INTEGER(iwp) ::  id_stg_south       !< south lateral boundary core id in case of turbulence generator
208    INTEGER(iwp) ::  stg_type_xz        !< MPI type for full z range
209    INTEGER(iwp) ::  stg_type_xz_small  !< MPI type for small z range
210    INTEGER(iwp) ::  stg_type_yz        !< MPI type for full z range
211    INTEGER(iwp) ::  stg_type_yz_small  !< MPI type for small z range
212    INTEGER(iwp) ::  merg               !< maximum length scale (in gp)
213    INTEGER(iwp) ::  mergp              !< merg + nbgp
214    INTEGER(iwp) ::  nzb_x_stg          !< lower bound of z coordinate (required for transposing z on PEs along x)
215    INTEGER(iwp) ::  nzt_x_stg          !< upper bound of z coordinate (required for transposing z on PEs along x)
216    INTEGER(iwp) ::  nzb_y_stg          !< lower bound of z coordinate (required for transposing z on PEs along y)
217    INTEGER(iwp) ::  nzt_y_stg          !< upper bound of z coordinate (required for transposing z on PEs along y)
218
219    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displs_xz      !< displacement for MPI_GATHERV
220    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  recv_count_xz  !< receive count for MPI_GATHERV
221    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displs_yz      !< displacement for MPI_GATHERV
222    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  recv_count_yz  !< receive count for MPI_GATHERV
223    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nux            !< length scale of u in x direction (in gp)
224    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nuy            !< length scale of u in y direction (in gp)
225    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nuz            !< length scale of u in z direction (in gp)
226    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nvx            !< length scale of v in x direction (in gp)
227    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nvy            !< length scale of v in y direction (in gp)
228    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nvz            !< length scale of v in z direction (in gp)
229    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwx            !< length scale of w in x direction (in gp)
230    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwy            !< length scale of w in y direction (in gp)
231    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwz            !< length scale of w in z direction (in gp)
232    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  seed           !< seed of random number for rn-generator
233
234    REAL(wp) ::  cosmo_ref = 10.0_wp      !< height of first vertical grid level in mesoscale model, used for calculation of scaling parameters
235    REAL(wp) ::  dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics
236    REAL(wp) ::  dt_stg_call = 5.0_wp     !< time interval for calling synthetic turbulence generator
237    REAL(wp) ::  mc_factor                !< mass flux correction factor
238    REAL(wp) ::  scale_l                  !< scaling parameter used for turbulence parametrization - Obukhov length
239    REAL(wp) ::  scale_us                 !< scaling parameter used for turbulence parametrization - friction velocity
240    REAL(wp) ::  scale_wm                 !< scaling parameter used for turbulence parametrization - momentum scale 
241    REAL(wp) ::  time_stg_adjust = 0.0_wp !< time counter for adjusting turbulence information   
242    REAL(wp) ::  time_stg_call = 0.0_wp   !< time counter for calling generator   
243   
244   
245    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r11              !< Reynolds parameter
246    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r21              !< Reynolds parameter
247    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r22              !< Reynolds parameter
248    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r31              !< Reynolds parameter
249    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r32              !< Reynolds parameter
250    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r33              !< Reynolds parameter
251   
252    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a11             !< coefficient for Lund rotation
253    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a21             !< coefficient for Lund rotation
254    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a22             !< coefficient for Lund rotation
255    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a31             !< coefficient for Lund rotation
256    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a32             !< coefficient for Lund rotation
257    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a33             !< coefficient for Lund rotation
258    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tu              !< Lagrangian time scale of u
259    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tv              !< Lagrangian time scale of v
260    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tw              !< Lagrangian time scale of w
261
262    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bux           !< filter function for u in x direction
263    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  buy           !< filter function for u in y direction
264    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  buz           !< filter function for u in z direction
265    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bvx           !< filter function for v in x direction
266    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bvy           !< filter function for v in y direction
267    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bvz           !< filter function for v in z direction
268    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bwx           !< filter function for w in y direction
269    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bwy           !< filter function for w in y direction
270    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bwz           !< filter function for w in z direction
271    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fu_xz         !< velocity seed for u at xz plane
272    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fuo_xz        !< velocity seed for u at xz plane with new random number
273    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fu_yz         !< velocity seed for u at yz plane
274    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fuo_yz        !< velocity seed for u at yz plane with new random number
275    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fv_xz         !< velocity seed for v at xz plane
276    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fvo_xz        !< velocity seed for v at xz plane with new random number
277    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fv_yz         !< velocity seed for v at yz plane
278    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fvo_yz        !< velocity seed for v at yz plane with new random number
279    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fw_xz         !< velocity seed for w at xz plane
280    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fwo_xz        !< velocity seed for w at xz plane with new random number
281    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fw_yz         !< velocity seed for w at yz plane
282    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fwo_yz        !< velocity seed for w at yz plane with new random number
283   
284    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_xz     !< imposed disturbances at north/south boundary
285    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_yz     !< imposed disturbances at north/south boundary
286
287!
288!-- PALM interfaces:
289!-- Adjust time and lenght scales, Reynolds stress, and filter functions
290    INTERFACE stg_adjust
291       MODULE PROCEDURE stg_adjust
292    END INTERFACE stg_adjust
293!
294!-- Input parameter checks to be done in check_parameters
295    INTERFACE stg_check_parameters
296       MODULE PROCEDURE stg_check_parameters
297    END INTERFACE stg_check_parameters
298
299!
300!-- Calculate filter functions
301    INTERFACE stg_filter_func
302       MODULE PROCEDURE stg_filter_func
303    END INTERFACE stg_filter_func
304
305!
306!-- Generate velocity seeds at south and north domain boundary
307    INTERFACE stg_generate_seed_xz
308       MODULE PROCEDURE stg_generate_seed_xz
309    END INTERFACE stg_generate_seed_xz
310!
311!-- Generate velocity seeds at left and/or right domain boundary
312    INTERFACE stg_generate_seed_yz
313       MODULE PROCEDURE stg_generate_seed_yz
314    END INTERFACE stg_generate_seed_yz
315
316!
317!-- Output of information to the header file
318    INTERFACE stg_header
319       MODULE PROCEDURE stg_header
320    END INTERFACE stg_header
321
322!
323!-- Initialization actions
324    INTERFACE stg_init
325       MODULE PROCEDURE stg_init
326    END INTERFACE stg_init
327
328!
329!-- Main procedure of synth. turb. gen.
330    INTERFACE stg_main
331       MODULE PROCEDURE stg_main
332    END INTERFACE stg_main
333
334!
335!-- Reading of NAMELIST parameters
336    INTERFACE stg_parin
337       MODULE PROCEDURE stg_parin
338    END INTERFACE stg_parin
339
340!
341!-- Reading of parameters for restart runs
342    INTERFACE stg_rrd_global
343       MODULE PROCEDURE stg_rrd_global
344    END INTERFACE stg_rrd_global
345
346!
347!-- Writing of binary output for restart runs
348    INTERFACE stg_wrd_global
349       MODULE PROCEDURE stg_wrd_global
350    END INTERFACE stg_wrd_global
351
352    SAVE
353
354    PRIVATE
355
356!
357!-- Public interfaces
358    PUBLIC  stg_adjust, stg_check_parameters, stg_header, stg_init, stg_main,  &
359            stg_parin, stg_rrd_global, stg_wrd_global
360
361!
362!-- Public variables
363    PUBLIC  dt_stg_call, dt_stg_adjust, id_stg_left, id_stg_north,             &
364            id_stg_right, id_stg_south, parametrize_inflow_turbulence,         &
365            time_stg_adjust, time_stg_call, use_syn_turb_gen
366
367
368 CONTAINS
369
370
371!------------------------------------------------------------------------------!
372! Description:
373! ------------
374!> Check parameters routine for synthetic turbulence generator
375!------------------------------------------------------------------------------!
376 SUBROUTINE stg_check_parameters
377
378
379    USE control_parameters,                                                    &
380        ONLY:  bc_lr, bc_ns, child_domain, nesting_offline, rans_mode,         &
381               turbulent_inflow
382
383    USE pmc_interface,                                                         &
384        ONLY : rans_mode_parent
385
386
387    IMPLICIT NONE
388
389    IF ( .NOT. use_syn_turb_gen  .AND.  .NOT. rans_mode  .AND.                 &
390          nesting_offline )  THEN
391       message_string = 'No synthetic turbulence generator is applied. ' //    &
392                        'In case PALM operates in LES mode and lateral ' //    &
393                        'boundary conditions are provided by COSMO model, ' // &
394                        'turbulence may require large adjustment lenght at ' //&
395                        'the lateral inflow boundaries. Please check your ' // &
396                        'results carefully.' 
397       CALL message( 'stg_check_parameters', 'PA0000', 0, 0, 0, 6, 0 )
398    ENDIF
399
400    IF ( .NOT. use_syn_turb_gen  .AND.  child_domain                           &
401         .AND. rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
402       message_string = 'Synthetic turbulence generator has to be applied ' // &
403                        'when nesting is applied and parent operates in '  //  &
404                        'RANS-mode but current child in LES mode.'
405       CALL message( 'stg_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
406    ENDIF
407
408    IF ( use_syn_turb_gen )  THEN
409
410       IF ( .NOT. nesting_offline  .AND.  .NOT. child_domain )  THEN
411       
412          IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0     &
413        .AND.  INDEX( initializing_actions, 'read_restart_data' ) == 0 )  THEN
414             message_string = 'Using synthetic turbulence generator ' //       &
415                              'requires %initializing_actions = '         //   &
416                              '"set_constant_profiles" or "read_restart_data"' //&
417                              ', if not offline nesting is applied.'
418             CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
419          ENDIF
420
421          IF ( bc_lr /= 'dirichlet/radiation' )  THEN
422             message_string = 'Using synthetic turbulence generator ' //       &
423                              'requires &bc_lr = "dirichlet/radiation", ' //   &
424                              'if not offline nesting is applied.'
425             CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
426          ENDIF
427          IF ( bc_ns /= 'cyclic' )  THEN
428             message_string = 'Using synthetic turbulence generator ' //       &
429                              'requires &bc_ns = "cyclic", ' //               &
430                              'if not offline nesting is applied.'
431             CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
432          ENDIF
433
434       ENDIF
435
436       IF ( turbulent_inflow )  THEN
437          message_string = 'Using synthetic turbulence generator ' //          &
438                           'in combination &with turbulent_inflow = .T. '//    &
439                              'is not allowed'
440          CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
441       ENDIF
442
443    ENDIF
444
445 END SUBROUTINE stg_check_parameters
446
447
448!------------------------------------------------------------------------------!
449! Description:
450! ------------
451!> Header output for synthetic turbulence generator
452!------------------------------------------------------------------------------!
453 SUBROUTINE stg_header ( io )
454
455
456    IMPLICIT NONE
457
458    INTEGER(iwp), INTENT(IN) ::  io   !< Unit of the output file
459
460!
461!-- Write synthetic turbulence generator Header
462    WRITE( io, 1 )
463    IF ( use_syn_turb_gen )  THEN
464       WRITE( io, 2 )
465    ELSE
466       WRITE( io, 3 )
467    ENDIF
468   
469    IF ( parametrize_inflow_turbulence )  THEN
470       WRITE( io, 4 ) dt_stg_adjust
471    ELSE
472       WRITE( io, 5 )
473    ENDIF
474
4751   FORMAT (//' Synthetic turbulence generator information:'/                  &
476              ' ------------------------------------------'/)
4772   FORMAT ('    synthetic turbulence generator is switched on')
4783   FORMAT ('    synthetic turbulence generator is switched off')
4794   FORMAT ('    imposed turbulence statistics are parametrized and ajdusted to boundary-layer development each ', F8.2, ' s' )
4805   FORMAT ('    imposed turbulence is read from file' )
481
482 END SUBROUTINE stg_header
483
484
485!------------------------------------------------------------------------------!
486! Description:
487! ------------
488!> Initialization of the synthetic turbulence generator
489!------------------------------------------------------------------------------!
490 SUBROUTINE stg_init
491
492
493    USE arrays_3d,                                                             &
494        ONLY:  dzw, ddzw, u_init, v_init, zu
495
496    USE control_parameters,                                                    &
497        ONLY:  child_domain, coupling_char, e_init, nesting_offline, rans_mode
498
499    USE grid_variables,                                                        &
500        ONLY:  ddx, ddy, dx, dy
501
502    USE indices,                                                               &
503        ONLY:  nz
504
505    USE pmc_interface,                                                         &
506        ONLY : rans_mode_parent
507
508
509    IMPLICIT NONE
510
511    LOGICAL ::  file_stg_exist = .FALSE. !< flag indicating whether parameter file for Reynolds stress and length scales exist
512
513#if defined( __parallel )
514    INTEGER(KIND=MPI_ADDRESS_KIND) :: extent !< extent of new MPI type
515    INTEGER(KIND=MPI_ADDRESS_KIND) :: tob=0  !< dummy variable
516#endif
517
518    INTEGER(iwp) :: i                        !> grid index in x-direction
519    INTEGER(iwp) :: j                        !> loop index
520    INTEGER(iwp) :: k                        !< index
521    INTEGER(iwp) :: newtype                  !< dummy MPI type
522    INTEGER(iwp) :: realsize                 !< size of REAL variables
523    INTEGER(iwp) :: nseed                    !< dimension of random number seed
524    INTEGER(iwp) :: startseed = 1234567890   !< start seed for random number generator
525
526!
527!-- Dummy variables used for reading profiles
528    REAL(wp) :: d1      !< u profile
529    REAL(wp) :: d2      !< v profile
530    REAL(wp) :: d3      !< w profile
531    REAL(wp) :: d5      !< e profile
532    REAL(wp) :: d11     !< vertical interpolation for a11
533    REAL(wp) :: d21     !< vertical interpolation for a21
534    REAL(wp) :: d22     !< vertical interpolation for a22
535    REAL(wp) :: luy     !< length scale for u in y direction
536    REAL(wp) :: luz     !< length scale for u in z direction
537    REAL(wp) :: lvy     !< length scale for v in y direction
538    REAL(wp) :: lvz     !< length scale for v in z direction
539    REAL(wp) :: lwy     !< length scale for w in y direction
540    REAL(wp) :: lwz     !< length scale for w in z direction
541    REAL(wp) :: nnz     !< increment used to determine processor decomposition of z-axis along x and y direction
542    REAL(wp) :: zz      !< height
543
544
545#if defined( __parallel )
546    CALL MPI_BARRIER( comm2d, ierr )
547#endif
548
549    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
550
551#if defined( __parallel )
552!
553!-- Determine processor decomposition of z-axis along x- and y-direction
554    nnz = nz / pdims(1)
555    nzb_x_stg = 1 + myidx * INT( nnz )
556    nzt_x_stg = ( myidx + 1 ) * INT( nnz )
557
558    IF ( MOD( nz , pdims(1) ) /= 0  .AND.  myidx == id_stg_right )             &
559       nzt_x_stg = nzt_x_stg + myidx * ( nnz - INT( nnz ) )
560
561    IF ( nesting_offline   .OR.  ( child_domain  .AND.  rans_mode_parent       &
562                            .AND.  .NOT.  rans_mode ) )  THEN
563       nnz = nz / pdims(2)
564       nzb_y_stg = 1 + myidy * INT( nnz )
565       nzt_y_stg = ( myidy + 1 ) * INT( nnz )
566
567       IF ( MOD( nz , pdims(2) ) /= 0  .AND.  myidy == id_stg_north )          &
568          nzt_y_stg = nzt_y_stg + myidy * ( nnz - INT( nnz ) )
569    ENDIF
570
571!
572!-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted
573!-- velocity seeds
574    CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
575    extent = 1 * realsize
576!
577!-- Set-up MPI datatyp to involve all cores for turbulence generation at yz
578!-- layer
579!-- stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
580    CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyng-nysg+1],                 &
581            [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
582    CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz, ierr )
583    CALL MPI_TYPE_COMMIT( stg_type_yz, ierr )
584    CALL MPI_TYPE_FREE( newtype, ierr )
585
586    ! stg_type_yz_small: yz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
587    CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyng-nysg+1],     &
588            [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
589    CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr )
590    CALL MPI_TYPE_COMMIT( stg_type_yz_small, ierr )
591    CALL MPI_TYPE_FREE( newtype, ierr )
592
593    ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
594    ALLOCATE( recv_count_yz(pdims(1)), displs_yz(pdims(1)) )
595
596    recv_count_yz           = nzt_x_stg-nzb_x_stg + 1
597    recv_count_yz(pdims(1)) = recv_count_yz(pdims(1)) + 1
598
599    DO  j = 1, pdims(1)
600       displs_yz(j) = 0 + (nzt_x_stg-nzb_x_stg+1) * (j-1)
601    ENDDO
602!
603!-- Set-up MPI datatyp to involve all cores for turbulence generation at xz
604!-- layer
605!-- stg_type_xz: xz-slice with vertical bounds nzb:nzt+1
606    IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
607                           .AND.  .NOT.  rans_mode ) )  THEN
608       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nxrg-nxlg+1],              &
609               [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
610       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz, ierr )
611       CALL MPI_TYPE_COMMIT( stg_type_xz, ierr )
612       CALL MPI_TYPE_FREE( newtype, ierr )
613
614       ! stg_type_yz_small: xz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
615       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_y_stg-nzb_y_stg+2,nxrg-nxlg+1],  &
616               [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
617       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz_small, ierr )
618       CALL MPI_TYPE_COMMIT( stg_type_xz_small, ierr )
619       CALL MPI_TYPE_FREE( newtype, ierr )
620
621       ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
622       ALLOCATE( recv_count_xz(pdims(2)), displs_xz(pdims(2)) )
623
624       recv_count_xz           = nzt_y_stg-nzb_y_stg + 1
625       recv_count_xz(pdims(2)) = recv_count_xz(pdims(2)) + 1
626
627       DO  j = 1, pdims(2)
628          displs_xz(j) = 0 + (nzt_y_stg-nzb_y_stg+1) * (j-1)
629       ENDDO
630
631    ENDIF
632
633#endif
634!
635!-- Define seed of random number
636    CALL RANDOM_SEED()
637    CALL RANDOM_SEED( size=nseed )
638    ALLOCATE( seed(1:nseed) )
639    DO  j = 1, nseed
640       seed(j) = startseed + j
641    ENDDO
642    CALL RANDOM_SEED(put=seed)
643!
644!-- Allocate required arrays
645!-- mean_inflow profiles must not be allocated in offline nesting
646    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
647       IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )                          &
648          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
649    ENDIF
650
651    ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                 &
652               a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1),                 &
653               nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1),                 &
654               nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1),                 &
655               nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
656               r11(nzb:nzt+1), r21(nzb:nzt+1), r22(nzb:nzt+1),                 &
657               r31(nzb:nzt+1), r32(nzb:nzt+1), r33(nzb:nzt+1),                 &
658               tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
659               
660    ALLOCATE ( dist_xz(nzb:nzt+1,nxlg:nxrg,3) )
661    ALLOCATE ( dist_yz(nzb:nzt+1,nysg:nyng,3) )
662    dist_xz = 0.0_wp
663    dist_yz = 0.0_wp
664!
665!-- Read inflow profile
666!-- Height levels of profiles in input profile is as follows:
667!-- zu: luy, luz, tu, lvy, lvz, tv, r11, r21, r22, d1, d2, d5
668!-- zw: lwy, lwz, tw, r31, r32, r33, d3
669!-- WARNING: zz is not used at the moment
670    INQUIRE( FILE = 'STG_PROFILES' // TRIM( coupling_char ),                   &
671             EXIST = file_stg_exist  )
672
673    IF ( file_stg_exist )  THEN
674
675       OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD',     &
676                      FORM='FORMATTED')
677!
678!--    Skip header
679       READ( 90, * )
680
681       DO  k = nzb+1, nzt+1
682          READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k), &
683                        r11(k), r21(k), r22(k), r31(k), r32(k), r33(k),        &
684                        d1, d2, d3, d5
685
686!
687!--       Convert length scales from meter to number of grid points.
688          nuy(k) = INT( luy * ddy )
689          nuz(k) = INT( luz * ddzw(k) )
690          nvy(k) = INT( lvy * ddy )
691          nvz(k) = INT( lvz * ddzw(k) )
692          nwy(k) = INT( lwy * ddy )
693          nwz(k) = INT( lwz * ddzw(k) )
694!
695!--       Workaround, assume isotropic turbulence
696          nwx(k) = nwy(k)
697          nvx(k) = nvy(k)
698          nux(k) = nuy(k)
699!
700!--       Save Mean inflow profiles
701          IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
702             mean_inflow_profiles(k,1) = d1
703             mean_inflow_profiles(k,2) = d2
704            !  mean_inflow_profiles(k,4) = d4
705             mean_inflow_profiles(k,5) = d5
706          ENDIF
707       ENDDO
708!
709!--    Set lenght scales at surface grid point
710       nuy(nzb) = nuy(nzb+1) 
711       nuz(nzb) = nuz(nzb+1)
712       nvy(nzb) = nvy(nzb+1)
713       nvz(nzb) = nvz(nzb+1)
714       nwy(nzb) = nwy(nzb+1)
715       nwz(nzb) = nwz(nzb+1)
716       
717       CLOSE( 90 )
718!
719!--    Calculate coefficient matrix from Reynolds stress tensor 
720!--    (Lund rotation)
721       CALL calc_coeff_matrix
722!
723!-- No information about turbulence and its length scales are available.
724!-- Instead, parametrize turbulence which is imposed at the boundaries.
725!-- Set flag which indicates that turbulence is parametrized, which is done
726!-- later when energy-balance models are already initialized. This is
727!-- because the STG needs information about surface properties such as
728!-- roughness to build 'realistic' turbulence profiles.
729    ELSE
730!
731!--    Set flag indicating that turbulence is parametrized
732       parametrize_inflow_turbulence = .TRUE.
733!
734!--    Determine boundary-layer depth, which is used to initialize lenght
735!--    scales
736       CALL calc_scaling_variables
737!
738!--    Initialize lenght and time scales, which in turn are used
739!--    to initialize the filter functions.
740       CALL calc_length_and_time_scale
741           
742    ENDIF
743
744!
745!-- Assign initial profiles. Note, this is only required if turbulent
746!-- inflow from the left is desired, not in case of any of the
747!-- nesting (offline or self nesting) approaches.
748    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
749       u_init = mean_inflow_profiles(:,1)
750       v_init = mean_inflow_profiles(:,2)
751      !pt_init = mean_inflow_profiles(:,4)
752       e_init = MAXVAL( mean_inflow_profiles(:,5) )
753    ENDIF
754   
755!
756!-- Define the size of the filter functions and allocate them.
757    merg = 0
758
759    ! arrays must be large enough to cover the largest length scale
760    DO  k = nzb, nzt+1
761       j = MAX( ABS(nux(k)), ABS(nuy(k)), ABS(nuz(k)), &
762                ABS(nvx(k)), ABS(nvy(k)), ABS(nvz(k)), &
763                ABS(nwx(k)), ABS(nwy(k)), ABS(nwz(k))  )
764       IF ( j > merg )  merg = j
765    ENDDO
766
767    merg  = 2 * merg
768    mergp = merg + nbgp
769
770    ALLOCATE ( bux(-merg:merg,nzb:nzt+1),                                      &
771               buy(-merg:merg,nzb:nzt+1),                                      &
772               buz(-merg:merg,nzb:nzt+1),                                      &
773               bvx(-merg:merg,nzb:nzt+1),                                      &
774               bvy(-merg:merg,nzb:nzt+1),                                      &
775               bvz(-merg:merg,nzb:nzt+1),                                      &
776               bwx(-merg:merg,nzb:nzt+1),                                      &
777               bwy(-merg:merg,nzb:nzt+1),                                      &
778               bwz(-merg:merg,nzb:nzt+1)  )
779
780!
781!-- Allocate velocity seeds for turbulence at xz-layer
782    ALLOCATE ( fu_xz( nzb:nzt+1,nxlg:nxrg), fuo_xz(nzb:nzt+1,nxlg:nxrg),       &
783               fv_xz( nzb:nzt+1,nxlg:nxrg), fvo_xz(nzb:nzt+1,nxlg:nxrg),       &
784               fw_xz( nzb:nzt+1,nxlg:nxrg), fwo_xz(nzb:nzt+1,nxlg:nxrg)  )
785
786!
787!-- Allocate velocity seeds for turbulence at yz-layer
788    ALLOCATE ( fu_yz( nzb:nzt+1,nysg:nyng), fuo_yz(nzb:nzt+1,nysg:nyng),       &
789               fv_yz( nzb:nzt+1,nysg:nyng), fvo_yz(nzb:nzt+1,nysg:nyng),       &
790               fw_yz( nzb:nzt+1,nysg:nyng), fwo_yz(nzb:nzt+1,nysg:nyng)  )
791
792    fu_xz  = 0.0_wp
793    fuo_xz = 0.0_wp
794    fv_xz  = 0.0_wp
795    fvo_xz = 0.0_wp
796    fw_xz  = 0.0_wp
797    fwo_xz = 0.0_wp
798
799    fu_yz  = 0.0_wp
800    fuo_yz = 0.0_wp
801    fv_yz  = 0.0_wp
802    fvo_yz = 0.0_wp
803    fw_yz  = 0.0_wp
804    fwo_yz = 0.0_wp
805
806!
807!-- Create filter functions
808    CALL stg_filter_func( nux, bux ) !filter ux
809    CALL stg_filter_func( nuy, buy ) !filter uy
810    CALL stg_filter_func( nuz, buz ) !filter uz
811    CALL stg_filter_func( nvx, bvx ) !filter vx
812    CALL stg_filter_func( nvy, bvy ) !filter vy
813    CALL stg_filter_func( nvz, bvz ) !filter vz
814    CALL stg_filter_func( nwx, bwx ) !filter wx
815    CALL stg_filter_func( nwy, bwy ) !filter wy
816    CALL stg_filter_func( nwz, bwz ) !filter wz
817
818#if defined( __parallel )
819    CALL MPI_BARRIER( comm2d, ierr )
820#endif
821
822!
823!-- In case of restart, calculate velocity seeds fu, fv, fw from former
824!   time step.
825!   Bug: fu, fv, fw are different in those heights where a11, a22, a33
826!        are 0 compared to the prerun. This is mostly for k=nzt+1.
827    IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
828       IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
829
830          IF ( myidx == id_stg_left  )  i = -1
831          IF ( myidx == id_stg_right )  i = nxr+1
832
833          DO  j = nysg, nyng
834             DO  k = nzb, nzt+1
835
836                IF  ( a11(k) .NE. 0._wp ) THEN
837                   fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
838                ELSE
839                   fu_yz(k,j) = 0._wp
840                ENDIF
841
842                IF  ( a22(k) .NE. 0._wp ) THEN
843                   fv_yz(k,j) = ( v(k,j,i) / mc_factor - a21(k) * fu_yz(k,j) - &
844                               v_init(k) ) / a22(k)
845                ELSE
846                   fv_yz(k,j) = 0._wp
847                ENDIF
848
849                IF  ( a33(k) .NE. 0._wp ) THEN
850                   fw_yz(k,j) = ( w(k,j,i) / mc_factor - a31(k) * fu_yz(k,j) - &
851                               a32(k) * fv_yz(k,j) ) / a33(k)
852                ELSE
853                   fw_yz = 0._wp
854                ENDIF
855
856             ENDDO
857          ENDDO
858       ENDIF
859    ENDIF
860
861    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
862
863 END SUBROUTINE stg_init
864
865
866!------------------------------------------------------------------------------!
867! Description:
868! ------------
869!> Calculate filter function bxx from length scale nxx following Eg.9 and 10
870!> (Xie and Castro, 2008)
871!------------------------------------------------------------------------------!
872 SUBROUTINE stg_filter_func( nxx, bxx )
873
874
875    IMPLICIT NONE
876
877    INTEGER(iwp) :: k         !< loop index
878    INTEGER(iwp) :: n_k       !< length scale nXX in height k
879    INTEGER(iwp) :: n_k2      !< n_k * 2
880    INTEGER(iwp) :: nf        !< index for length scales
881
882    REAL(wp) :: bdenom        !< denominator for filter functions bXX
883    REAL(wp) :: qsi = 1.0_wp  !< minimization factor
884
885    INTEGER(iwp), DIMENSION(:) :: nxx(nzb:nzt+1)           !< length scale (in gp)
886
887    REAL(wp), DIMENSION(:,:) :: bxx(-merg:merg,nzb:nzt+1)  !< filter function
888
889
890    bxx = 0.0_wp
891
892    DO  k = nzb, nzt+1
893       bdenom = 0.0_wp
894       n_k    = nxx(k)
895       IF ( n_k /= 0 )  THEN
896          n_k2 = n_k * 2
897
898!
899!--       ( Eq.10 )^2
900          DO  nf = -n_k2, n_k2
901             bdenom = bdenom + EXP( -qsi * pi * ABS(nf) / n_k )**2
902          ENDDO
903
904!
905!--       ( Eq.9 )
906          bdenom = SQRT( bdenom )
907          DO  nf = -n_k2, n_k2
908             bxx(nf,k) = EXP( -qsi * pi * ABS(nf) / n_k ) / bdenom
909          ENDDO
910       ENDIF
911    ENDDO
912
913 END SUBROUTINE stg_filter_func
914
915
916!------------------------------------------------------------------------------!
917! Description:
918! ------------
919!> Parin for &stg_par for synthetic turbulence generator
920!------------------------------------------------------------------------------!
921 SUBROUTINE stg_parin
922
923
924    IMPLICIT NONE
925
926    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
927
928
929    NAMELIST /stg_par/  dt_stg_adjust, dt_stg_call, use_syn_turb_gen
930
931    line = ' '
932
933!
934!-- Try to find stg package
935    REWIND ( 11 )
936    line = ' '
937    DO WHILE ( INDEX( line, '&stg_par' ) == 0 )
938       READ ( 11, '(A)', END=20 )  line
939    ENDDO
940    BACKSPACE ( 11 )
941
942!
943!-- Read namelist
944    READ ( 11, stg_par, ERR = 10, END = 20 )
945
946!
947!-- Set flag that indicates that the synthetic turbulence generator is switched
948!-- on
949    syn_turb_gen = .TRUE.
950    GOTO 20
951
952 10 BACKSPACE( 11 )
953    READ( 11 , '(A)') line
954    CALL parin_fail_message( 'stg_par', line )
955
956 20 CONTINUE
957
958 END SUBROUTINE stg_parin
959
960
961!------------------------------------------------------------------------------!
962! Description:
963! ------------
964!> This routine reads the respective restart data.
965!------------------------------------------------------------------------------!
966 SUBROUTINE stg_rrd_global( found )
967
968
969    USE control_parameters,                                                    &
970        ONLY: length, restart_string
971
972
973    IMPLICIT NONE
974
975    LOGICAL, INTENT(OUT)  ::  found !< flag indicating if variable was found
976
977
978    found = .TRUE.
979
980
981    SELECT CASE ( restart_string(1:length) )
982
983       CASE ( 'mc_factor' )
984          READ ( 13 )  mc_factor
985       CASE ( 'use_syn_turb_gen' )
986          READ ( 13 )  use_syn_turb_gen
987
988       CASE DEFAULT
989
990          found = .FALSE.
991
992    END SELECT
993
994
995 END SUBROUTINE stg_rrd_global
996
997
998!------------------------------------------------------------------------------!
999! Description:
1000! ------------
1001!> This routine writes the respective restart data.
1002!------------------------------------------------------------------------------!
1003 SUBROUTINE stg_wrd_global
1004
1005
1006    IMPLICIT NONE
1007
1008    CALL wrd_write_string( 'mc_factor' )
1009    WRITE ( 14 )  mc_factor
1010
1011    CALL wrd_write_string( 'use_syn_turb_gen' )
1012    WRITE ( 14 )  use_syn_turb_gen
1013
1014
1015 END SUBROUTINE stg_wrd_global
1016
1017
1018!------------------------------------------------------------------------------!
1019! Description:
1020! ------------
1021!> Create turbulent inflow fields for u, v, w with prescribed length scales and
1022!> Reynolds stress tensor after a method of Xie and Castro (2008), modified
1023!> following suggestions of Kim et al. (2013), and using a Lund rotation
1024!> (Lund, 1998).
1025!------------------------------------------------------------------------------!
1026 SUBROUTINE stg_main
1027
1028
1029    USE arrays_3d,                                                             &
1030        ONLY:  dzw
1031
1032    USE control_parameters,                                                    &
1033        ONLY:  child_domain, dt_3d, intermediate_timestep_count,               &
1034               nesting_offline, rans_mode, simulated_time, volume_flow_initial
1035
1036    USE grid_variables,                                                        &
1037        ONLY:  dx, dy
1038
1039    USE indices,                                                               &
1040        ONLY:  wall_flags_0
1041
1042    USE statistics,                                                            &
1043        ONLY:  weight_substep
1044
1045    USE pmc_interface,                                                         &
1046        ONLY : rans_mode_parent
1047
1048
1049    IMPLICIT NONE
1050
1051    INTEGER(iwp) :: i           !< grid index in x-direction
1052    INTEGER(iwp) :: j           !< loop index in y-direction
1053    INTEGER(iwp) :: k           !< loop index in z-direction
1054
1055    REAL(wp) :: dt_stg          !< wheighted subtimestep
1056    REAL(wp) :: mc_factor_l     !< local mass flux correction factor
1057    REAL(wp) :: volume_flow     !< mass flux through lateral boundary
1058    REAL(wp) :: volume_flow_l   !< local mass flux through lateral boundary
1059
1060    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1061
1062!
1063!-- Calculate time step which is needed for filter functions
1064    dt_stg = MAX( dt_3d, dt_stg_call ) !* weight_substep(intermediate_timestep_count)
1065!
1066!-- Initial value of fu, fv, fw
1067    IF ( simulated_time == 0.0_wp .AND. .NOT. velocity_seed_initialized )  THEN
1068       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
1069       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
1070       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
1071
1072       IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent      &
1073                                .AND.  .NOT.  rans_mode ) )  THEN
1074!
1075!--       Generate turbulence at right boundary
1076          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_right )
1077          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_right )
1078          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_right )
1079!
1080!--       Generate turbulence at north boundary
1081          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_north )
1082          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_north )
1083          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_north )
1084!
1085!--       Generate turbulence at south boundary
1086          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_south )
1087          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_south )
1088          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_south )
1089       ENDIF
1090       velocity_seed_initialized = .TRUE.
1091    ENDIF
1092   
1093!
1094!-- New set of fu, fv, fw
1095    CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
1096    CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
1097    CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
1098   
1099    IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
1100                          .AND.  .NOT.  rans_mode ) )  THEN
1101!
1102!--    Generate turbulence at right boundary
1103       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
1104       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
1105       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
1106!
1107!--    Generate turbulence at north boundary
1108       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
1109       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
1110       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
1111!
1112!--    Generate turbulence at south boundary
1113       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
1114       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
1115       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
1116    ENDIF
1117   
1118!
1119!-- Turbulence generation at left and or right boundary
1120    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
1121   
1122       DO  j = nysg, nyng
1123          DO  k = nzb, nzt + 1
1124!
1125!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1126             IF ( tu(k) == 0.0_wp )  THEN
1127                fu_yz(k,j) = fuo_yz(k,j)
1128             ELSE
1129                fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
1130                         fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
1131             ENDIF
1132
1133             IF ( tv(k) == 0.0_wp )  THEN
1134                fv_yz(k,j) = fvo_yz(k,j)
1135             ELSE
1136                fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
1137                         fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
1138             ENDIF
1139
1140             IF ( tw(k) == 0.0_wp )  THEN
1141                fw_yz(k,j) = fwo_yz(k,j)
1142             ELSE
1143                fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
1144                         fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
1145             ENDIF
1146!
1147!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
1148!--          Additional factors are added to improve the variance of v and w
1149             IF( k == 0 )  THEN
1150                dist_yz(k,j,1) = 0.0_wp
1151                dist_yz(k,j,2) = 0.0_wp
1152                dist_yz(k,j,3) = 0.0_wp
1153             ELSE
1154                dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp )
1155                !experimental test of 1.2
1156                dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
1157                                      * 1.2_wp )                               &
1158                                     * (   a21(k) * fu_yz(k,j)                 &
1159                                         + a22(k) * fv_yz(k,j) ), 3.0_wp )
1160                dist_yz(k,j,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
1161                                      * 1.3_wp )                               &
1162                                    * (   a31(k) * fu_yz(k,j)                  &
1163                                        + a32(k) * fv_yz(k,j)                  &
1164                                        + a33(k) * fw_yz(k,j) ), 3.0_wp )
1165             ENDIF
1166
1167          ENDDO
1168       ENDDO
1169
1170!
1171!--    Mass flux correction following Kim et al. (2013)
1172!--    This correction factor insures that the mass flux is preserved at the
1173!--    inflow boundary
1174       IF ( .NOT. nesting_offline  .AND.  .NOT. child_domain )  THEN
1175          mc_factor_l = 0.0_wp
1176          mc_factor   = 0.0_wp
1177          DO  j = nys, nyn
1178             DO  k = nzb+1, nzt
1179                mc_factor_l = mc_factor_l + dzw(k)  *                          &
1180                              ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) )
1181             ENDDO
1182          ENDDO
1183
1184#if defined( __parallel )
1185          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 1, MPI_REAL, MPI_SUM,    &
1186                              comm1dy, ierr )
1187#else
1188          mc_factor = mc_factor_l
1189#endif
1190
1191          mc_factor = volume_flow_initial(1) / mc_factor
1192
1193!
1194!--       Add disturbance at the inflow
1195          DO  j = nysg, nyng
1196             DO  k = nzb, nzt+1
1197                 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +              &
1198                                      dist_yz(k,j,1)             ) * mc_factor
1199                 v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +              &
1200                                      dist_yz(k,j,2)             ) * mc_factor
1201                 w(k,j,-nbgp:-1)  =   dist_yz(k,j,3)               * mc_factor
1202             ENDDO
1203          ENDDO
1204
1205       ELSE
1206!
1207!--       First, calculate volume flow at yz boundary
1208          IF ( myidx == id_stg_left  )  i = nxl
1209          IF ( myidx == id_stg_right )  i = nxr+1
1210
1211          volume_flow_l = 0.0_wp
1212          volume_flow   = 0.0_wp
1213          mc_factor_l   = 0.0_wp
1214          mc_factor     = 0.0_wp
1215          DO  j = nys, nyn
1216             DO  k = nzb+1, nzt
1217                volume_flow_l = volume_flow_l + u(k,j,i) * dzw(k) * dy         &
1218                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1219                                              BTEST( wall_flags_0(k,j,i), 1 ) )
1220
1221                mc_factor_l = mc_factor_l     + ( u(k,j,i) + dist_yz(k,j,1) )  &
1222                                                         * dzw(k) * dy         &
1223                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1224                                              BTEST( wall_flags_0(k,j,i), 1 ) )
1225             ENDDO
1226          ENDDO
1227#if defined( __parallel )
1228          CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                      &
1229                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1230          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                          &
1231                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1232#else
1233          volume_flow = volume_flow_l
1234          mc_factor   = mc_factor_l
1235#endif
1236
1237          mc_factor = volume_flow / mc_factor
1238
1239!
1240!--       Add disturbances
1241          IF ( myidx == id_stg_left  )  THEN
1242             DO  j = nys, nyn
1243                DO  k = nzb+1, nzt
1244                   u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) )                    &
1245                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
1246                                            BTEST( wall_flags_0(k,j,0), 1 ) )
1247                   u(k,j,-1) = u(k,j,0)
1248                   v(k,j,-1)  = ( v(k,j,-1)  + dist_yz(k,j,2)  )               &
1249                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
1250                                            BTEST( wall_flags_0(k,j,-1), 2 ) )
1251                   w(k,j,-1)  = ( w(k,j,-1)  + dist_yz(k,j,3)  )               &
1252                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
1253                                            BTEST( wall_flags_0(k,j,-1), 3 ) )
1254                ENDDO
1255             ENDDO
1256          ENDIF
1257          IF ( myidx == id_stg_right  )  THEN
1258             DO  j = nys, nyn
1259                DO  k = nzb+1, nzt
1260                   u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) )            &
1261                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1262                                          BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
1263                   v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) )            &
1264                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1265                                          BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
1266                   w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) )            &
1267                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1268                                          BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
1269                ENDDO
1270             ENDDO
1271          ENDIF
1272       ENDIF
1273
1274    ENDIF
1275!
1276!-- Turbulence generation at north and south boundary
1277    IF ( myidy == id_stg_north  .OR.  myidy == id_stg_south )  THEN
1278   
1279       DO  i = nxlg, nxrg
1280          DO  k = nzb, nzt + 1
1281!
1282!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1283             IF ( tu(k) == 0.0_wp )  THEN
1284                fu_xz(k,i) = fuo_xz(k,i)
1285             ELSE
1286                fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
1287                         fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
1288             ENDIF
1289
1290             IF ( tv(k) == 0.0_wp )  THEN
1291                fv_xz(k,i) = fvo_xz(k,i)
1292             ELSE
1293                fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
1294                      fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
1295             ENDIF
1296           
1297             IF ( tw(k) == 0.0_wp )  THEN
1298                fw_xz(k,i) = fwo_xz(k,i)
1299             ELSE
1300                fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
1301                      fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
1302             ENDIF
1303!
1304!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
1305!--          Additional factors are added to improve the variance of v and w
1306             IF( k == 0 )  THEN
1307                dist_xz(k,i,1) = 0.0_wp
1308                dist_xz(k,i,2) = 0.0_wp
1309                dist_xz(k,i,3) = 0.0_wp
1310
1311             ELSE
1312                dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp )
1313                !experimental test of 1.2
1314                dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
1315                                      * 1.2_wp )                               &
1316                                     * (   a21(k) * fu_xz(k,i)                 &
1317                                         + a22(k) * fv_xz(k,i) ), 3.0_wp )
1318                dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
1319                                      * 1.3_wp )                               &
1320                                    * (   a31(k) * fu_xz(k,i)                  &
1321                                        + a32(k) * fv_xz(k,i)                  &
1322                                        + a33(k) * fw_xz(k,i) ), 3.0_wp )
1323             ENDIF
1324
1325          ENDDO
1326       ENDDO
1327
1328!
1329!--    Mass flux correction following Kim et al. (2013)
1330!--    This correction factor insures that the mass flux is preserved at the
1331!--    inflow boundary.
1332!--    First, calculate volume flow at xz boundary
1333       IF ( myidy == id_stg_south  ) j = nys
1334       IF ( myidy == id_stg_north )  j = nyn+1
1335
1336       volume_flow_l = 0.0_wp
1337       volume_flow   = 0.0_wp
1338       mc_factor_l   = 0.0_wp
1339       mc_factor     = 0.0_wp
1340       DO  i = nxl, nxr
1341          DO  k = nzb+1, nzt
1342             volume_flow_l = volume_flow_l + v(k,j,i) * dzw(k) * dx            &
1343                                  * MERGE( 1.0_wp, 0.0_wp,                     &
1344                                           BTEST( wall_flags_0(k,j,i), 2 ) )
1345
1346             mc_factor_l = mc_factor_l     + ( v(k,j,i) + dist_xz(k,i,2) )     &
1347                                                      * dzw(k) * dx            &
1348                                  * MERGE( 1.0_wp, 0.0_wp,                     &
1349                                           BTEST( wall_flags_0(k,j,i), 2 ) )
1350          ENDDO
1351       ENDDO
1352#if defined( __parallel )
1353       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                         &
1354                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1355       CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                             &
1356                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1357#else
1358       volume_flow = volume_flow_l
1359       mc_factor   = mc_factor_l
1360#endif
1361
1362       mc_factor = volume_flow / mc_factor
1363
1364!
1365!--    Add disturbances
1366       IF ( myidy == id_stg_south  )  THEN
1367
1368          DO  i = nxl, nxr
1369             DO  k = nzb+1, nzt
1370                u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) )                     &
1371                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
1372                                           BTEST( wall_flags_0(k,-1,i), 1 ) )
1373                v(k,0,i)  = ( v(k,0,i)  + dist_xz(k,i,2)  )                    &
1374                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
1375                                           BTEST( wall_flags_0(k,0,i), 2 ) )
1376                v(k,-1,i) = v(k,0,i)
1377                w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3)  )                    &
1378                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
1379                                           BTEST( wall_flags_0(k,-1,i), 3 ) )
1380             ENDDO
1381          ENDDO
1382       ENDIF
1383       IF ( myidy == id_stg_north  )  THEN
1384
1385          DO  i = nxl, nxr
1386             DO  k = nzb+1, nzt
1387                u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )               &
1388                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1389                                          BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
1390                v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) )               &
1391                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1392                                          BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
1393                w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) )               &
1394                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1395                                          BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
1396             ENDDO
1397          ENDDO
1398       ENDIF
1399    ENDIF
1400!
1401!-- Finally, set time counter for calling STG to zero
1402    time_stg_call = 0.0_wp
1403
1404    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1405
1406 END SUBROUTINE stg_main
1407
1408!------------------------------------------------------------------------------!
1409! Description:
1410! ------------
1411!> Generate a set of random number rand_it wich is equal on each PE
1412!> and calculate the velocity seed f_n.
1413!> f_n is splitted in vertical direction by the number of PEs in x-direction and
1414!> and each PE calculates a vertical subsection of f_n. At the the end, all
1415!> parts are collected to form the full array.
1416!------------------------------------------------------------------------------!
1417 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
1418
1419
1420    USE indices,                                                               &
1421        ONLY: ny
1422
1423    IMPLICIT NONE
1424
1425    INTEGER(iwp) :: id          !< core ids at respective boundaries
1426    INTEGER(iwp) :: j           !< loop index in y-direction
1427    INTEGER(iwp) :: jj          !< loop index in y-direction
1428    INTEGER(iwp) :: k           !< loop index in z-direction
1429    INTEGER(iwp) :: kk          !< loop index in z-direction
1430    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
1431
1432    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y    !< length scale in y-direction
1433    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
1434    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y2   !< n_y*2
1435    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
1436
1437    REAL(wp) :: nyz_inv         !< inverse of number of grid points in yz-slice
1438    REAL(wp) :: rand_av         !< average of random number
1439    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
1440
1441    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter function in y-direction
1442    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
1443    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
1444    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)     :: f_n     !<  velocity seed
1445    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
1446
1447
1448!
1449!-- Generate random numbers using a seed generated in stg_init.
1450!-- The set of random numbers are modified to have an average of 0 and
1451!-- unit variance.
1452    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:ny+mergp) )
1453
1454    rand_av        = 0.0_wp
1455    rand_sigma_inv = 0.0_wp
1456    nyz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( ny+1 ), KIND=wp )
1457
1458    DO  j = 0, ny
1459       DO  k = nzb, nzt+1
1460          CALL RANDOM_NUMBER( rand_it(k,j) )
1461          rand_av = rand_av + rand_it(k,j)
1462       ENDDO
1463    ENDDO
1464
1465    rand_av = rand_av * nyz_inv
1466
1467    DO  j = 0, ny
1468       DO  k = nzb, nzt+1
1469          rand_it(k,j)   = rand_it(k,j) - rand_av
1470          rand_sigma_inv = rand_sigma_inv + rand_it(k,j) ** 2
1471       ENDDO
1472    ENDDO
1473
1474    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nyz_inv)
1475
1476    DO  j = 0, ny
1477       DO  k = nzb, nzt+1
1478          rand_it(k,j) = rand_it(k,j) * rand_sigma_inv
1479       ENDDO
1480    ENDDO
1481
1482!
1483!-- Periodic fill of random number in space
1484    DO  j = 0, ny
1485       DO  k = 1, mergp
1486          rand_it(nzb  -k,j) = rand_it(nzt+2-k,j)    ! bottom margin
1487          rand_it(nzt+1+k,j) = rand_it(nzb+k-1,j)    ! top margin
1488       ENDDO
1489    ENDDO
1490    DO  j = 1, mergp
1491       DO  k = nzb-mergp, nzt+1+mergp
1492          rand_it(k,  -j) = rand_it(k,ny-j+1)        ! south margin
1493          rand_it(k,ny+j) = rand_it(k,   j-1)        ! north margin
1494       ENDDO
1495    ENDDO
1496
1497!
1498!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1499    n_y2 = n_y * 2
1500    n_z2 = n_z * 2
1501    f_n_l  = 0.0_wp
1502
1503    DO  j = nysg, nyng
1504       DO  k = nzb_x_stg, nzt_x_stg+1
1505          DO  jj = -n_y2(k), n_y2(k)
1506             DO  kk = -n_z2(k), n_z2(k)
1507                f_n_l(k,j) = f_n_l(k,j)                                        &
1508                           + b_y(jj,k) * b_z(kk,k) * rand_it(k+kk,j+jj)
1509             ENDDO
1510          ENDDO
1511       ENDDO
1512    ENDDO
1513
1514    DEALLOCATE( rand_it )
1515!
1516!-- Gather velocity seeds of full subdomain
1517    send_count = nzt_x_stg - nzb_x_stg + 1
1518    IF ( nzt_x_stg == nzt )  send_count = send_count + 1
1519
1520#if defined( __parallel )
1521    CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,    &
1522                      f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,  &
1523                      id, comm1dx, ierr )
1524#else
1525    f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x_stg:nzt_x_stg+1,nysg:nyng)
1526#endif
1527
1528
1529 END SUBROUTINE stg_generate_seed_yz
1530
1531
1532
1533
1534!------------------------------------------------------------------------------!
1535! Description:
1536! ------------
1537!> Generate a set of random number rand_it wich is equal on each PE
1538!> and calculate the velocity seed f_n.
1539!> f_n is splitted in vertical direction by the number of PEs in y-direction and
1540!> and each PE calculates a vertical subsection of f_n. At the the end, all
1541!> parts are collected to form the full array.
1542!------------------------------------------------------------------------------!
1543 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
1544
1545
1546    USE indices,                                                               &
1547        ONLY: nx
1548
1549
1550    IMPLICIT NONE
1551
1552    INTEGER(iwp) :: id          !< core ids at respective boundaries
1553    INTEGER(iwp) :: i           !< loop index in x-direction
1554    INTEGER(iwp) :: ii          !< loop index in x-direction
1555    INTEGER(iwp) :: k           !< loop index in z-direction
1556    INTEGER(iwp) :: kk          !< loop index in z-direction
1557    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
1558
1559    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x    !< length scale in x-direction
1560    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
1561    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_x*2
1562    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
1563
1564    REAL(wp) :: nxz_inv         !< inverse of number of grid points in xz-slice
1565    REAL(wp) :: rand_av         !< average of random number
1566    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
1567
1568    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter function in x-direction
1569    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
1570    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
1571    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)     :: f_n     !<  velocity seed
1572    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
1573
1574
1575!
1576!-- Generate random numbers using a seed generated in stg_init.
1577!-- The set of random numbers are modified to have an average of 0 and
1578!-- unit variance.
1579    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:nx+mergp) )
1580
1581    rand_av        = 0.0_wp
1582    rand_sigma_inv = 0.0_wp
1583    nxz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( nx+1 ), KIND=wp )
1584
1585    DO  i = 0, nx
1586       DO  k = nzb, nzt+1
1587          CALL RANDOM_NUMBER( rand_it(k,i) )
1588          rand_av = rand_av + rand_it(k,i)
1589       ENDDO
1590    ENDDO
1591
1592    rand_av = rand_av * nxz_inv
1593
1594    DO  i = 0, nx
1595       DO  k = nzb, nzt+1
1596          rand_it(k,i)   = rand_it(k,i) - rand_av
1597          rand_sigma_inv = rand_sigma_inv + rand_it(k,i) ** 2
1598       ENDDO
1599    ENDDO
1600
1601    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nxz_inv)
1602
1603    DO  i = 0, nx
1604       DO  k = nzb, nzt+1
1605          rand_it(k,i) = rand_it(k,i) * rand_sigma_inv
1606       ENDDO
1607    ENDDO
1608
1609!
1610!-- Periodic fill of random number in space
1611    DO  i = 0, nx
1612       DO  k = 1, mergp
1613          rand_it(nzb-k,i)   = rand_it(nzt+2-k,i)    ! bottom margin
1614          rand_it(nzt+1+k,i) = rand_it(nzb+k-1,i)    ! top margin
1615       ENDDO
1616    ENDDO
1617    DO  i = 1, mergp
1618       DO  k = nzb-mergp, nzt+1+mergp
1619          rand_it(k,-i)   = rand_it(k,nx-i+1)        ! left margin
1620          rand_it(k,nx+i) = rand_it(k,i-1)           ! right margin
1621       ENDDO
1622    ENDDO
1623
1624!
1625!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1626    n_x2 = n_x * 2
1627    n_z2 = n_z * 2
1628    f_n_l  = 0.0_wp
1629
1630    DO  i = nxlg, nxrg
1631       DO  k = nzb_y_stg, nzt_y_stg+1
1632          DO  ii = -n_x2(k), n_x2(k)
1633             DO  kk = -n_z2(k), n_z2(k)
1634                f_n_l(k,i) = f_n_l(k,i)                                        &
1635                           + b_x(ii,k) * b_z(kk,k) * rand_it(k+kk,i+ii)
1636             ENDDO
1637          ENDDO
1638       ENDDO
1639    ENDDO
1640
1641    DEALLOCATE( rand_it )
1642
1643!
1644!-- Gather velocity seeds of full subdomain
1645    send_count = nzt_y_stg - nzb_y_stg + 1
1646    IF ( nzt_y_stg == nzt )  send_count = send_count + 1
1647
1648
1649#if defined( __parallel )
1650    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,    &
1651                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz,  &
1652                      id, comm1dy, ierr )
1653#else
1654    f_n(nzb+1:nzt+1,nxlg:nxrg) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg)
1655#endif
1656
1657
1658 END SUBROUTINE stg_generate_seed_xz
1659
1660!------------------------------------------------------------------------------!
1661! Description:
1662! ------------
1663!> Parametrization of the Reynolds stress tensor, following the parametrization
1664!> described in Rotach et al. (1996), which is applied in state-of-the-art
1665!> dispserion modelling. Please note, the parametrization does not distinguish
1666!> between along-wind and cross-wind turbulence.
1667!------------------------------------------------------------------------------!
1668 SUBROUTINE parametrize_reynolds_stress
1669
1670    USE arrays_3d,                                                             &
1671        ONLY:  zu
1672
1673    IMPLICIT NONE
1674
1675    INTEGER(iwp) :: k            !< loop index in z-direction
1676   
1677    REAL(wp)     ::  zzi         !< ratio of z/zi
1678   
1679!
1680!--
1681    DO  k = nzb+1, nzt+1
1682         
1683       IF ( zu(k) <= zi_ribulk )  THEN
1684!
1685!--       Determine normalized height coordinate
1686          zzi = zu(k) / zi_ribulk
1687!
1688!--       u'u' and v'v'. Assume isotropy. Note, add a small negative number
1689!--       to the denominator, else the merge-function can crash if scale_l is
1690!--       zero.
1691          r11(k) = scale_us**2 * (                                             &
1692                      MERGE( 0.35_wp * (                                       &
1693                           - zi_ribulk / ( kappa * scale_l - 10E-4_wp )        &
1694                                       )**( 2.0_wp / 3.0_wp ),                 &
1695                             0.0_wp,                                           &
1696                             scale_l < 0.0_wp )                                &
1697                    + 5.0_wp - 4.0_wp * zzi                                    &
1698                                 )
1699                                 
1700          r22(k) = r11(k)
1701!
1702!--       w'w'
1703          r33(k) = scale_wm**2 * (                                             &
1704                      1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi ) &
1705                    + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2            &                     
1706                                 )
1707!
1708!--       u'w' and v'w'. Assume isotropy.
1709          r31(k) = - scale_us**2 * (                                           &
1710                         1.0_wp - EXP( 3.0_wp * ( zzi - 1.0_wp ) )             &
1711                                   )
1712                             
1713          r32(k) = r31(k)
1714!
1715!--       For u'v' no parametrization exist so far - ?. For simplicity assume
1716!--       a similar profile as for u'w'.
1717          r21(k) = r31(k)
1718!
1719!--    Above the boundary layer, assmume laminar flow conditions.
1720       ELSE
1721          r11(k) = 10E-8_wp
1722          r22(k) = 10E-8_wp
1723          r33(k) = 10E-8_wp
1724          r21(k) = 10E-8_wp
1725          r31(k) = 10E-8_wp
1726          r32(k) = 10E-8_wp
1727       ENDIF
1728!        write(9,*) zu(k), r11(k), r33(k), r31(k), zi_ribulk, scale_us, scale_wm, scale_l
1729    ENDDO
1730
1731!
1732!-- Set bottom boundary condition   
1733    r11(nzb) = r11(nzb+1)
1734    r22(nzb) = r22(nzb+1)
1735    r33(nzb) = r33(nzb+1)
1736
1737    r21(nzb) = r11(nzb+1)
1738    r31(nzb) = r31(nzb+1)
1739    r32(nzb) = r32(nzb+1)
1740   
1741
1742 END SUBROUTINE parametrize_reynolds_stress 
1743 
1744!------------------------------------------------------------------------------!
1745! Description:
1746! ------------
1747!> Calculate the coefficient matrix from the Lund rotation.
1748!------------------------------------------------------------------------------!
1749 SUBROUTINE calc_coeff_matrix
1750
1751    IMPLICIT NONE
1752
1753    INTEGER(iwp) :: k   !< loop index in z-direction
1754   
1755!
1756!-- Calculate coefficient matrix. Split loops to allow for loop vectorization.
1757    DO  k = nzb+1, nzt+1
1758       IF ( r11(k) > 0.0_wp )  THEN
1759          a11(k) = SQRT( r11(k) ) 
1760          a21(k) = r21(k) / a11(k)
1761          a31(k) = r31(k) / a11(k)
1762       ELSE
1763          a11(k) = 10E-8_wp
1764          a21(k) = 10E-8_wp
1765          a31(k) = 10E-8_wp
1766       ENDIF
1767    ENDDO
1768    DO  k = nzb+1, nzt+1
1769       a22(k) = r22(k) - a21(k)**2 
1770       IF ( a22(k) > 0.0_wp )  THEN
1771          a22(k) = SQRT( a22(k) )
1772          a32(k) = r32(k) - a21(k) * a31(k) / a22(k)
1773       ELSE
1774          a22(k) = 10E-8_wp
1775          a32(k) = 10E-8_wp
1776       ENDIF
1777    ENDDO
1778    DO  k = nzb+1, nzt+1
1779       a33(k) = r33(k) - a31(k)**2 - a32(k)**2
1780       IF ( a33(k) > 0.0_wp )  THEN
1781          a33(k) =  SQRT( a33(k) )
1782       ELSE
1783          a33(k) = 10E-8_wp
1784       ENDIF
1785    ENDDO   
1786!
1787!-- Set bottom boundary condition
1788    a11(nzb) = a11(nzb+1)
1789    a22(nzb) = a22(nzb+1)
1790    a21(nzb) = a21(nzb+1)
1791    a33(nzb) = a33(nzb+1)   
1792    a31(nzb) = a31(nzb+1)
1793    a32(nzb) = a32(nzb+1)   
1794
1795 END SUBROUTINE calc_coeff_matrix
1796 
1797!------------------------------------------------------------------------------!
1798! Description:
1799! ------------
1800!> This routine controls the re-adjustment of the turbulence statistics used
1801!> for generating turbulence at the lateral boundaries. 
1802!------------------------------------------------------------------------------!
1803 SUBROUTINE stg_adjust
1804
1805    IMPLICIT NONE
1806   
1807    INTEGER(iwp) ::  k
1808
1809    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1810!
1811!-- Compute mean boundary layer height according to Richardson-Bulk
1812!-- criterion using the inflow profiles. Further velocity scale as well as
1813!-- mean friction velocity are calculated.
1814    CALL calc_scaling_variables
1815!
1816!-- Set length and time scales depending on boundary-layer height
1817    CALL calc_length_and_time_scale
1818!
1819!-- Parametrize Reynolds-stress tensor, diagonal elements as well
1820!-- as r21 (v'u'), r31 (w'u'), r32 (w'v'). Parametrization follows
1821!-- Rotach et al. (1996) and is based on boundary-layer depth,
1822!-- friction velocity and velocity scale.
1823    CALL parametrize_reynolds_stress
1824!
1825!-- Calculate coefficient matrix from Reynolds stress tensor 
1826!-- (Lund rotation)
1827    CALL calc_coeff_matrix
1828!
1829!-- Determine filter functions on basis of updated length scales
1830    CALL stg_filter_func( nux, bux ) !filter ux
1831    CALL stg_filter_func( nuy, buy ) !filter uy
1832    CALL stg_filter_func( nuz, buz ) !filter uz
1833    CALL stg_filter_func( nvx, bvx ) !filter vx
1834    CALL stg_filter_func( nvy, bvy ) !filter vy
1835    CALL stg_filter_func( nvz, bvz ) !filter vz
1836    CALL stg_filter_func( nwx, bwx ) !filter wx
1837    CALL stg_filter_func( nwy, bwy ) !filter wy
1838    CALL stg_filter_func( nwz, bwz ) !filter wz
1839!
1840!-- Reset time counter for controlling next adjustment to zero
1841    time_stg_adjust = 0.0_wp
1842   
1843    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1844   
1845 END SUBROUTINE stg_adjust 
1846 
1847 
1848!------------------------------------------------------------------------------!
1849! Description:
1850! ------------
1851!> Calculates turbuluent length and time scales if these are not available
1852!> from measurements.
1853!------------------------------------------------------------------------------!
1854 SUBROUTINE calc_length_and_time_scale
1855
1856    USE arrays_3d,                                                             &
1857        ONLY:  dzw, ddzw, u_init, v_init, zu, zw
1858
1859    USE grid_variables,                                                        &
1860        ONLY:  ddx, ddy, dx, dy
1861       
1862    USE statistics,                                                            &
1863        ONLY:  hom
1864
1865    IMPLICIT NONE
1866
1867
1868    INTEGER(iwp) :: k            !< loop index in z-direction
1869
1870    REAL(wp) ::  length_scale    !< typical length scale
1871   
1872!
1873!-- In initial call the boundary-layer depth can be zero. This case, set
1874!-- minimum value for boundary-layer depth, to setup length scales correctly.
1875    zi_ribulk = MAX( zi_ribulk, zw(nzb+2) )
1876!
1877!-- Set-up default turbulent length scales. From the numerical point of
1878!-- view the imposed perturbations should not be immediately dissipated
1879!-- by the numerics. The numerical dissipation, however, acts on scales
1880!-- up to 8 x the grid spacing. For this reason, set the turbulence
1881!-- length scale to 8 time the grid spacing. Further, above the boundary
1882!-- layer height, set turbulence lenght scales to zero (equivalent to not
1883!-- imposing any perturbations) in order to save computational costs.
1884!-- Typical time scales are derived by assuming Taylors's hypothesis,
1885!-- using the length scales and the mean profiles of u- and v-component.
1886    DO  k = nzb+1, nzt+1
1887   
1888       length_scale = 8.0_wp * MIN( dx, dy, dzw(k) )
1889         
1890       IF ( zu(k) <= zi_ribulk )  THEN 
1891!
1892!--       Assume isotropic turbulence length scales
1893          nux(k) = MAX( INT( length_scale * ddx     ), 1 )
1894          nuy(k) = MAX( INT( length_scale * ddy     ), 1 )
1895          nuz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
1896          nvx(k) = MAX( INT( length_scale * ddx     ), 1 )
1897          nvy(k) = MAX( INT( length_scale * ddy     ), 1 )
1898          nvz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
1899          nwx(k) = MAX( INT( length_scale * ddx     ), 1 )
1900          nwy(k) = MAX( INT( length_scale * ddy     ), 1 )
1901          nwz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
1902!
1903!--       Limit time scales, else they become very larger for low wind speed,
1904!--       imposing long-living inflow perturbations which in turn propagate
1905!--       further into the model domain. Use u_init and v_init to calculate
1906!--       the time scales, which will be equal to the inflow profiles, both,
1907!--       in offline nesting mode or in dirichlet/radiation mode.
1908          tu(k)  = MIN( dt_stg_adjust, length_scale /                          &
1909                                  ( ABS( u_init(k) ) + 0.1_wp ) )
1910          tv(k)  = MIN( dt_stg_adjust, length_scale /                          &
1911                                  ( ABS( v_init(k) ) + 0.1_wp ) )
1912!
1913!--       Time scale of w-component is a mixture from u- and v-component.
1914          tw(k)  = SQRT( tu(k)**2 + tv(k)**2 )
1915!
1916!--    Above the boundary layer length scales are zero, i.e. imposed turbulence
1917!--    is not correlated in space and time, just white noise. This saves
1918!--    computations power.
1919       ELSE
1920          nux(k) = 0.0_wp
1921          nuy(k) = 0.0_wp
1922          nuz(k) = 0.0_wp
1923          nvx(k) = 0.0_wp
1924          nvy(k) = 0.0_wp
1925          nvz(k) = 0.0_wp
1926          nwx(k) = 0.0_wp
1927          nwy(k) = 0.0_wp
1928          nwz(k) = 0.0_wp
1929         
1930          tu(k)  = 0.0_wp
1931          tv(k)  = 0.0_wp
1932          tw(k)  = 0.0_wp
1933       ENDIF
1934    ENDDO
1935!
1936!-- Set bottom boundary condition for the length and time scales
1937    nux(nzb) = nux(nzb+1)
1938    nuy(nzb) = nuy(nzb+1)
1939    nuz(nzb) = nuz(nzb+1)
1940    nvx(nzb) = nvx(nzb+1)
1941    nvy(nzb) = nvy(nzb+1)
1942    nvz(nzb) = nvz(nzb+1)
1943    nwx(nzb) = nwx(nzb+1)
1944    nwy(nzb) = nwy(nzb+1)
1945    nwz(nzb) = nwz(nzb+1)
1946
1947    tu(nzb)  = tu(nzb+1)
1948    tv(nzb)  = tv(nzb+1)
1949    tw(nzb)  = tw(nzb+1)
1950
1951
1952 END SUBROUTINE calc_length_and_time_scale 
1953
1954!------------------------------------------------------------------------------!
1955! Description:
1956! ------------
1957!> Calculate scaling variables which are used for turbulence parametrization
1958!> according to Rotach et al. (1996). Scaling variables are: friction velocity,
1959!> boundary-layer depth, momentum velocity scale, and Obukhov length.
1960!------------------------------------------------------------------------------!
1961 SUBROUTINE calc_scaling_variables
1962
1963
1964    USE control_parameters,                                                    &
1965        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
1966               humidity, pt_surface
1967             
1968    USE indices,                                                               &
1969        ONLY:  nx, ny 
1970       
1971    USE surface_mod,                                                           &
1972        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
1973
1974    IMPLICIT NONE
1975
1976    INTEGER(iwp) :: i            !< loop index in x-direction
1977    INTEGER(iwp) :: j            !< loop index in y-direction
1978    INTEGER(iwp) :: k            !< loop index in z-direction
1979    INTEGER(iwp) :: k_ref        !< index in z-direction for reference height
1980    INTEGER(iwp) :: m            !< surface element index
1981
1982    REAL(wp) ::  friction_vel_l         !< mean friction veloctiy on subdomain
1983    REAL(wp) ::  shf_mean               !< mean surface sensible heat flux
1984    REAL(wp) ::  shf_mean_l             !< mean surface sensible heat flux on subdomain
1985    REAL(wp) ::  u_int                  !< u-component
1986    REAL(wp) ::  v_int                  !< v-component
1987    REAL(wp) ::  vpt_surface            !< near-surface virtual potential temperature
1988    REAL(wp) ::  w_convective           !< convective velocity scale
1989    REAL(wp) ::  z0_mean                !< mean roughness length
1990    REAL(wp) ::  z0_mean_l              !< mean roughness length on subdomain
1991   
1992!
1993!-- Mean friction velocity and velocity scale. Therefore,
1994!-- pre-calculate mean roughness length and surface sensible heat flux
1995!-- in the model domain, which are further used to estimate friction
1996!-- velocity and velocity scale. Note, for z0 linear averaging is applied,
1997!-- even though this is known to unestimate the effective roughness.
1998!-- This need to be revised later.
1999    z0_mean_l  = 0.0_wp
2000    shf_mean_l = 0.0_wp
2001    DO  m = 1, surf_def_h(0)%ns
2002       z0_mean_l  = z0_mean_l  + surf_def_h(0)%z0(m)
2003       shf_mean_l = shf_mean_l + surf_def_h(0)%shf(m)
2004    ENDDO   
2005    DO  m = 1, surf_lsm_h%ns
2006       z0_mean_l  = z0_mean_l  + surf_lsm_h%z0(m)
2007       shf_mean_l = shf_mean_l + surf_lsm_h%shf(m)
2008    ENDDO
2009    DO  m = 1, surf_usm_h%ns
2010       z0_mean_l  = z0_mean_l  + surf_usm_h%z0(m)
2011       shf_mean_l = shf_mean_l + surf_usm_h%shf(m)
2012    ENDDO
2013   
2014#if defined( __parallel )
2015    CALL MPI_ALLREDUCE( z0_mean_l, z0_mean, 1, MPI_REAL, MPI_SUM,              &
2016                        comm2d, ierr )
2017    CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM,            &
2018                        comm2d, ierr )
2019#else
2020    z0_mean  = z0_mean_l
2021    shf_mean = shf_mean_l
2022#endif
2023    z0_mean  = z0_mean  / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
2024    shf_mean = shf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
2025       
2026!
2027!-- Note, Inifor does not use logarithmic interpolation of the
2028!-- velocity components near the ground, so that near-surface
2029!-- wind speed and thus the friction velocity is overestimated.
2030!-- However, friction velocity is used for turbulence
2031!-- parametrization, so that more physically meaningful values are important.
2032!-- Hence, derive friction velocity from wind speed at a reference height.
2033!-- For a first guess use 20 m, which is in the range of the first
2034!-- COSMO vertical level.
2035    k_ref = MINLOC( ABS( zu - cosmo_ref ), DIM = 1 ) - 1
2036   
2037    friction_vel_l = 0.0_wp
2038    IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
2039       
2040       i = MERGE( -1, nxr + 1, bc_dirichlet_l )
2041
2042       DO  j = nys, nyn
2043!
2044!--       Note, in u- and v- component the imposed perturbations
2045!--       from the STG are already included. Check whether this
2046!--       makes any difference compared to using the pure-mean
2047!--       inflow profiles.
2048          u_int = MERGE( u(k_ref,j,i+1), u(k_ref,j,i), bc_dirichlet_l )
2049          v_int = v(k_ref,j,i)
2050!
2051!--       Calculate friction velocity and sum-up. Therefore, assume
2052!--       neutral condtions.
2053          friction_vel_l = friction_vel_l + kappa *                            &
2054                            SQRT( u_int * u_int + v_int * v_int )  /           &
2055                            LOG( zu(k_ref) / z0_mean )
2056         
2057       ENDDO
2058       
2059    ENDIF
2060   
2061    IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
2062       
2063       j = MERGE( -1, nyn + 1, bc_dirichlet_s )
2064   
2065       DO  i = nxl, nxr
2066         
2067          u_int = u(k_ref,j,i)
2068          v_int = MERGE( v(k_ref,j+1,i), v(k_ref,j,i), bc_dirichlet_s )
2069 
2070          friction_vel_l = friction_vel_l + kappa *                            &
2071                            SQRT( u_int * u_int + v_int * v_int )  /           &
2072                            LOG( zu(k_ref) / z0_mean )                       
2073         
2074       ENDDO
2075       
2076    ENDIF
2077   
2078#if defined( __parallel )
2079    CALL MPI_ALLREDUCE( friction_vel_l, scale_us, 1, MPI_REAL, MPI_SUM,        &
2080                        comm2d, ierr )
2081#else
2082    scale_us = friction_vel_l
2083#endif
2084    scale_us = scale_us / REAL( 2 * nx + 2 * ny, KIND = wp )
2085   
2086!
2087!-- Compute mean Obukhov length
2088    IF ( shf_mean > 0.0_wp )  THEN
2089       scale_l = - pt_surface / ( kappa * g ) * scale_us**3 / shf_mean
2090    ELSE
2091       scale_l = 0.0_wp
2092    ENDIF
2093   
2094!
2095!-- Compute mean convective velocity scale. Note, in case the mean heat flux
2096!-- is negative, set convective velocity scale to zero.
2097    IF ( shf_mean > 0.0_wp )  THEN
2098       w_convective = ( g * shf_mean * zi_ribulk / pt_surface )**( 1.0_wp / 3.0_wp )
2099    ELSE
2100       w_convective = 0.0_wp
2101    ENDIF
2102!
2103!-- Finally, in order to consider also neutral or stable stratification,
2104!-- compute momentum velocity scale from u* and convective velocity scale,
2105!-- according to Rotach et al. (1996).
2106    scale_wm = ( scale_us**3 + 0.6_wp * w_convective**3 )**( 1.0_wp / 3.0_wp )
2107   
2108 END SUBROUTINE calc_scaling_variables
2109
2110 END MODULE synthetic_turbulence_generator_mod
Note: See TracBrowser for help on using the repository browser.