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

Last change on this file since 3690 was 3646, checked in by kanani, 5 years ago

Bugfix: replace simulated_time by time_since_reference_point where required

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