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

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

Fix for format descriptor

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