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

Last change on this file since 3347 was 3347, checked in by suehring, 5 years ago

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

  • 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!
23!
24! Former revisions:
25! -----------------
26! $Id: synthetic_turbulence_generator_mod.f90 3347 2018-10-15 14:21:08Z 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 ' \\       &
477            'ajdusted to boundary-layer development each ', F8.2, ' s' )
4785   FORMAT ('    imposed turbulence is read from file' )
479
480 END SUBROUTINE stg_header
481
482
483!------------------------------------------------------------------------------!
484! Description:
485! ------------
486!> Initialization of the synthetic turbulence generator
487!------------------------------------------------------------------------------!
488 SUBROUTINE stg_init
489
490
491    USE arrays_3d,                                                             &
492        ONLY:  dzw, ddzw, u_init, v_init, zu
493
494    USE control_parameters,                                                    &
495        ONLY:  child_domain, coupling_char, e_init, nesting_offline, rans_mode
496
497    USE grid_variables,                                                        &
498        ONLY:  ddx, ddy, dx, dy
499
500    USE indices,                                                               &
501        ONLY:  nz
502
503    USE pmc_interface,                                                         &
504        ONLY : rans_mode_parent
505
506
507    IMPLICIT NONE
508
509    LOGICAL ::  file_stg_exist = .FALSE. !< flag indicating whether parameter file for Reynolds stress and length scales exist
510
511#if defined( __parallel )
512    INTEGER(KIND=MPI_ADDRESS_KIND) :: extent !< extent of new MPI type
513    INTEGER(KIND=MPI_ADDRESS_KIND) :: tob=0  !< dummy variable
514#endif
515
516    INTEGER(iwp) :: i                        !> grid index in x-direction
517    INTEGER(iwp) :: j                        !> loop index
518    INTEGER(iwp) :: k                        !< index
519    INTEGER(iwp) :: newtype                  !< dummy MPI type
520    INTEGER(iwp) :: realsize                 !< size of REAL variables
521    INTEGER(iwp) :: nseed                    !< dimension of random number seed
522    INTEGER(iwp) :: startseed = 1234567890   !< start seed for random number generator
523
524!
525!-- Dummy variables used for reading profiles
526    REAL(wp) :: d1      !< u profile
527    REAL(wp) :: d2      !< v profile
528    REAL(wp) :: d3      !< w profile
529    REAL(wp) :: d5      !< e profile
530    REAL(wp) :: d11     !< vertical interpolation for a11
531    REAL(wp) :: d21     !< vertical interpolation for a21
532    REAL(wp) :: d22     !< vertical interpolation for a22
533    REAL(wp) :: luy     !< length scale for u in y direction
534    REAL(wp) :: luz     !< length scale for u in z direction
535    REAL(wp) :: lvy     !< length scale for v in y direction
536    REAL(wp) :: lvz     !< length scale for v in z direction
537    REAL(wp) :: lwy     !< length scale for w in y direction
538    REAL(wp) :: lwz     !< length scale for w in z direction
539    REAL(wp) :: nnz     !< increment used to determine processor decomposition of z-axis along x and y direction
540    REAL(wp) :: zz      !< height
541
542
543#if defined( __parallel )
544    CALL MPI_BARRIER( comm2d, ierr )
545#endif
546
547    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
548
549#if defined( __parallel )
550!
551!-- Determine processor decomposition of z-axis along x- and y-direction
552    nnz = nz / pdims(1)
553    nzb_x_stg = 1 + myidx * INT( nnz )
554    nzt_x_stg = ( myidx + 1 ) * INT( nnz )
555
556    IF ( MOD( nz , pdims(1) ) /= 0  .AND.  myidx == id_stg_right )             &
557       nzt_x_stg = nzt_x_stg + myidx * ( nnz - INT( nnz ) )
558
559    IF ( nesting_offline   .OR.  ( child_domain  .AND.  rans_mode_parent       &
560                            .AND.  .NOT.  rans_mode ) )  THEN
561       nnz = nz / pdims(2)
562       nzb_y_stg = 1 + myidy * INT( nnz )
563       nzt_y_stg = ( myidy + 1 ) * INT( nnz )
564
565       IF ( MOD( nz , pdims(2) ) /= 0  .AND.  myidy == id_stg_north )          &
566          nzt_y_stg = nzt_y_stg + myidy * ( nnz - INT( nnz ) )
567    ENDIF
568
569!
570!-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted
571!-- velocity seeds
572    CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
573    extent = 1 * realsize
574!
575!-- Set-up MPI datatyp to involve all cores for turbulence generation at yz
576!-- layer
577!-- stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
578    CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyng-nysg+1],                 &
579            [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
580    CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz, ierr )
581    CALL MPI_TYPE_COMMIT( stg_type_yz, ierr )
582    CALL MPI_TYPE_FREE( newtype, ierr )
583
584    ! stg_type_yz_small: yz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
585    CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyng-nysg+1],     &
586            [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
587    CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr )
588    CALL MPI_TYPE_COMMIT( stg_type_yz_small, ierr )
589    CALL MPI_TYPE_FREE( newtype, ierr )
590
591    ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
592    ALLOCATE( recv_count_yz(pdims(1)), displs_yz(pdims(1)) )
593
594    recv_count_yz           = nzt_x_stg-nzb_x_stg + 1
595    recv_count_yz(pdims(1)) = recv_count_yz(pdims(1)) + 1
596
597    DO  j = 1, pdims(1)
598       displs_yz(j) = 0 + (nzt_x_stg-nzb_x_stg+1) * (j-1)
599    ENDDO
600!
601!-- Set-up MPI datatyp to involve all cores for turbulence generation at xz
602!-- layer
603!-- stg_type_xz: xz-slice with vertical bounds nzb:nzt+1
604    IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
605                           .AND.  .NOT.  rans_mode ) )  THEN
606       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nxrg-nxlg+1],              &
607               [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
608       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz, ierr )
609       CALL MPI_TYPE_COMMIT( stg_type_xz, ierr )
610       CALL MPI_TYPE_FREE( newtype, ierr )
611
612       ! stg_type_yz_small: xz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
613       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_y_stg-nzb_y_stg+2,nxrg-nxlg+1],  &
614               [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
615       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz_small, ierr )
616       CALL MPI_TYPE_COMMIT( stg_type_xz_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_xz(pdims(2)), displs_xz(pdims(2)) )
621
622       recv_count_xz           = nzt_y_stg-nzb_y_stg + 1
623       recv_count_xz(pdims(2)) = recv_count_xz(pdims(2)) + 1
624
625       DO  j = 1, pdims(2)
626          displs_xz(j) = 0 + (nzt_y_stg-nzb_y_stg+1) * (j-1)
627       ENDDO
628
629    ENDIF
630
631#endif
632!
633!-- Define seed of random number
634    CALL RANDOM_SEED()
635    CALL RANDOM_SEED( size=nseed )
636    ALLOCATE( seed(1:nseed) )
637    DO  j = 1, nseed
638       seed(j) = startseed + j
639    ENDDO
640    CALL RANDOM_SEED(put=seed)
641!
642!-- Allocate required arrays
643!-- mean_inflow profiles must not be allocated in offline nesting
644    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
645       IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )                          &
646          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
647    ENDIF
648
649    ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                 &
650               a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1),                 &
651               nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1),                 &
652               nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1),                 &
653               nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
654               r11(nzb:nzt+1), r21(nzb:nzt+1), r22(nzb:nzt+1),                 &
655               r31(nzb:nzt+1), r32(nzb:nzt+1), r33(nzb:nzt+1),                 &
656               tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
657               
658    ALLOCATE ( dist_xz(nzb:nzt+1,nxlg:nxrg,3) )
659    ALLOCATE ( dist_yz(nzb:nzt+1,nysg:nyng,3) )
660    dist_xz = 0.0_wp
661    dist_yz = 0.0_wp
662!
663!-- Read inflow profile
664!-- Height levels of profiles in input profile is as follows:
665!-- zu: luy, luz, tu, lvy, lvz, tv, r11, r21, r22, d1, d2, d5
666!-- zw: lwy, lwz, tw, r31, r32, r33, d3
667!-- WARNING: zz is not used at the moment
668    INQUIRE( FILE = 'STG_PROFILES' // TRIM( coupling_char ),                   &
669             EXIST = file_stg_exist  )
670
671    IF ( file_stg_exist )  THEN
672
673       OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD',     &
674                      FORM='FORMATTED')
675!
676!--    Skip header
677       READ( 90, * )
678
679       DO  k = nzb+1, nzt+1
680          READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k), &
681                        r11(k), r21(k), r22(k), r31(k), r32(k), r33(k),        &
682                        d1, d2, d3, d5
683
684!
685!--       Convert length scales from meter to number of grid points.
686          nuy(k) = INT( luy * ddy )
687          nuz(k) = INT( luz * ddzw(k) )
688          nvy(k) = INT( lvy * ddy )
689          nvz(k) = INT( lvz * ddzw(k) )
690          nwy(k) = INT( lwy * ddy )
691          nwz(k) = INT( lwz * ddzw(k) )
692!
693!--       Workaround, assume isotropic turbulence
694          nwx(k) = nwy(k)
695          nvx(k) = nvy(k)
696          nux(k) = nuy(k)
697!
698!--       Save Mean inflow profiles
699          IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
700             mean_inflow_profiles(k,1) = d1
701             mean_inflow_profiles(k,2) = d2
702            !  mean_inflow_profiles(k,4) = d4
703             mean_inflow_profiles(k,5) = d5
704          ENDIF
705       ENDDO
706!
707!--    Set lenght scales at surface grid point
708       nuy(nzb) = nuy(nzb+1) 
709       nuz(nzb) = nuz(nzb+1)
710       nvy(nzb) = nvy(nzb+1)
711       nvz(nzb) = nvz(nzb+1)
712       nwy(nzb) = nwy(nzb+1)
713       nwz(nzb) = nwz(nzb+1)
714       
715       CLOSE( 90 )
716!
717!--    Calculate coefficient matrix from Reynolds stress tensor 
718!--    (Lund rotation)
719       CALL calc_coeff_matrix
720!
721!-- No information about turbulence and its length scales are available.
722!-- Instead, parametrize turbulence which is imposed at the boundaries.
723!-- Set flag which indicates that turbulence is parametrized, which is done
724!-- later when energy-balance models are already initialized. This is
725!-- because the STG needs information about surface properties such as
726!-- roughness to build 'realistic' turbulence profiles.
727    ELSE
728!
729!--    Set flag indicating that turbulence is parametrized
730       parametrize_inflow_turbulence = .TRUE.
731!
732!--    Determine boundary-layer depth, which is used to initialize lenght
733!--    scales
734       CALL calc_scaling_variables
735!
736!--    Initialize lenght and time scales, which in turn are used
737!--    to initialize the filter functions.
738       CALL calc_length_and_time_scale
739           
740    ENDIF
741
742!
743!-- Assign initial profiles. Note, this is only required if turbulent
744!-- inflow from the left is desired, not in case of any of the
745!-- nesting (offline or self nesting) approaches.
746    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
747       u_init = mean_inflow_profiles(:,1)
748       v_init = mean_inflow_profiles(:,2)
749      !pt_init = mean_inflow_profiles(:,4)
750       e_init = MAXVAL( mean_inflow_profiles(:,5) )
751    ENDIF
752   
753!
754!-- Define the size of the filter functions and allocate them.
755    merg = 0
756
757    ! arrays must be large enough to cover the largest length scale
758    DO  k = nzb, nzt+1
759       j = MAX( ABS(nux(k)), ABS(nuy(k)), ABS(nuz(k)), &
760                ABS(nvx(k)), ABS(nvy(k)), ABS(nvz(k)), &
761                ABS(nwx(k)), ABS(nwy(k)), ABS(nwz(k))  )
762       IF ( j > merg )  merg = j
763    ENDDO
764
765    merg  = 2 * merg
766    mergp = merg + nbgp
767
768    ALLOCATE ( bux(-merg:merg,nzb:nzt+1),                                      &
769               buy(-merg:merg,nzb:nzt+1),                                      &
770               buz(-merg:merg,nzb:nzt+1),                                      &
771               bvx(-merg:merg,nzb:nzt+1),                                      &
772               bvy(-merg:merg,nzb:nzt+1),                                      &
773               bvz(-merg:merg,nzb:nzt+1),                                      &
774               bwx(-merg:merg,nzb:nzt+1),                                      &
775               bwy(-merg:merg,nzb:nzt+1),                                      &
776               bwz(-merg:merg,nzb:nzt+1)  )
777
778!
779!-- Allocate velocity seeds for turbulence at xz-layer
780    ALLOCATE ( fu_xz( nzb:nzt+1,nxlg:nxrg), fuo_xz(nzb:nzt+1,nxlg:nxrg),       &
781               fv_xz( nzb:nzt+1,nxlg:nxrg), fvo_xz(nzb:nzt+1,nxlg:nxrg),       &
782               fw_xz( nzb:nzt+1,nxlg:nxrg), fwo_xz(nzb:nzt+1,nxlg:nxrg)  )
783
784!
785!-- Allocate velocity seeds for turbulence at yz-layer
786    ALLOCATE ( fu_yz( nzb:nzt+1,nysg:nyng), fuo_yz(nzb:nzt+1,nysg:nyng),       &
787               fv_yz( nzb:nzt+1,nysg:nyng), fvo_yz(nzb:nzt+1,nysg:nyng),       &
788               fw_yz( nzb:nzt+1,nysg:nyng), fwo_yz(nzb:nzt+1,nysg:nyng)  )
789
790    fu_xz  = 0.0_wp
791    fuo_xz = 0.0_wp
792    fv_xz  = 0.0_wp
793    fvo_xz = 0.0_wp
794    fw_xz  = 0.0_wp
795    fwo_xz = 0.0_wp
796
797    fu_yz  = 0.0_wp
798    fuo_yz = 0.0_wp
799    fv_yz  = 0.0_wp
800    fvo_yz = 0.0_wp
801    fw_yz  = 0.0_wp
802    fwo_yz = 0.0_wp
803
804!
805!-- Create filter functions
806    CALL stg_filter_func( nux, bux ) !filter ux
807    CALL stg_filter_func( nuy, buy ) !filter uy
808    CALL stg_filter_func( nuz, buz ) !filter uz
809    CALL stg_filter_func( nvx, bvx ) !filter vx
810    CALL stg_filter_func( nvy, bvy ) !filter vy
811    CALL stg_filter_func( nvz, bvz ) !filter vz
812    CALL stg_filter_func( nwx, bwx ) !filter wx
813    CALL stg_filter_func( nwy, bwy ) !filter wy
814    CALL stg_filter_func( nwz, bwz ) !filter wz
815
816#if defined( __parallel )
817    CALL MPI_BARRIER( comm2d, ierr )
818#endif
819
820!
821!-- In case of restart, calculate velocity seeds fu, fv, fw from former
822!   time step.
823!   Bug: fu, fv, fw are different in those heights where a11, a22, a33
824!        are 0 compared to the prerun. This is mostly for k=nzt+1.
825    IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
826       IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
827
828          IF ( myidx == id_stg_left  )  i = -1
829          IF ( myidx == id_stg_right )  i = nxr+1
830
831          DO  j = nysg, nyng
832             DO  k = nzb, nzt+1
833
834                IF  ( a11(k) .NE. 0._wp ) THEN
835                   fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
836                ELSE
837                   fu_yz(k,j) = 0._wp
838                ENDIF
839
840                IF  ( a22(k) .NE. 0._wp ) THEN
841                   fv_yz(k,j) = ( v(k,j,i) / mc_factor - a21(k) * fu_yz(k,j) - &
842                               v_init(k) ) / a22(k)
843                ELSE
844                   fv_yz(k,j) = 0._wp
845                ENDIF
846
847                IF  ( a33(k) .NE. 0._wp ) THEN
848                   fw_yz(k,j) = ( w(k,j,i) / mc_factor - a31(k) * fu_yz(k,j) - &
849                               a32(k) * fv_yz(k,j) ) / a33(k)
850                ELSE
851                   fw_yz = 0._wp
852                ENDIF
853
854             ENDDO
855          ENDDO
856       ENDIF
857    ENDIF
858
859    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
860
861 END SUBROUTINE stg_init
862
863
864!------------------------------------------------------------------------------!
865! Description:
866! ------------
867!> Calculate filter function bxx from length scale nxx following Eg.9 and 10
868!> (Xie and Castro, 2008)
869!------------------------------------------------------------------------------!
870 SUBROUTINE stg_filter_func( nxx, bxx )
871
872
873    IMPLICIT NONE
874
875    INTEGER(iwp) :: k         !< loop index
876    INTEGER(iwp) :: n_k       !< length scale nXX in height k
877    INTEGER(iwp) :: n_k2      !< n_k * 2
878    INTEGER(iwp) :: nf        !< index for length scales
879
880    REAL(wp) :: bdenom        !< denominator for filter functions bXX
881    REAL(wp) :: qsi = 1.0_wp  !< minimization factor
882
883    INTEGER(iwp), DIMENSION(:) :: nxx(nzb:nzt+1)           !< length scale (in gp)
884
885    REAL(wp), DIMENSION(:,:) :: bxx(-merg:merg,nzb:nzt+1)  !< filter function
886
887
888    bxx = 0.0_wp
889
890    DO  k = nzb, nzt+1
891       bdenom = 0.0_wp
892       n_k    = nxx(k)
893       IF ( n_k /= 0 )  THEN
894          n_k2 = n_k * 2
895
896!
897!--       ( Eq.10 )^2
898          DO  nf = -n_k2, n_k2
899             bdenom = bdenom + EXP( -qsi * pi * ABS(nf) / n_k )**2
900          ENDDO
901
902!
903!--       ( Eq.9 )
904          bdenom = SQRT( bdenom )
905          DO  nf = -n_k2, n_k2
906             bxx(nf,k) = EXP( -qsi * pi * ABS(nf) / n_k ) / bdenom
907          ENDDO
908       ENDIF
909    ENDDO
910
911 END SUBROUTINE stg_filter_func
912
913
914!------------------------------------------------------------------------------!
915! Description:
916! ------------
917!> Parin for &stg_par for synthetic turbulence generator
918!------------------------------------------------------------------------------!
919 SUBROUTINE stg_parin
920
921
922    IMPLICIT NONE
923
924    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
925
926
927    NAMELIST /stg_par/  dt_stg_adjust, dt_stg_call, use_syn_turb_gen
928
929    line = ' '
930
931!
932!-- Try to find stg package
933    REWIND ( 11 )
934    line = ' '
935    DO WHILE ( INDEX( line, '&stg_par' ) == 0 )
936       READ ( 11, '(A)', END=20 )  line
937    ENDDO
938    BACKSPACE ( 11 )
939
940!
941!-- Read namelist
942    READ ( 11, stg_par, ERR = 10, END = 20 )
943
944!
945!-- Set flag that indicates that the synthetic turbulence generator is switched
946!-- on
947    syn_turb_gen = .TRUE.
948    GOTO 20
949
950 10 BACKSPACE( 11 )
951    READ( 11 , '(A)') line
952    CALL parin_fail_message( 'stg_par', line )
953
954 20 CONTINUE
955
956 END SUBROUTINE stg_parin
957
958
959!------------------------------------------------------------------------------!
960! Description:
961! ------------
962!> This routine reads the respective restart data.
963!------------------------------------------------------------------------------!
964 SUBROUTINE stg_rrd_global( found )
965
966
967    USE control_parameters,                                                    &
968        ONLY: length, restart_string
969
970
971    IMPLICIT NONE
972
973    LOGICAL, INTENT(OUT)  ::  found !< flag indicating if variable was found
974
975
976    found = .TRUE.
977
978
979    SELECT CASE ( restart_string(1:length) )
980
981       CASE ( 'mc_factor' )
982          READ ( 13 )  mc_factor
983       CASE ( 'use_syn_turb_gen' )
984          READ ( 13 )  use_syn_turb_gen
985
986       CASE DEFAULT
987
988          found = .FALSE.
989
990    END SELECT
991
992
993 END SUBROUTINE stg_rrd_global
994
995
996!------------------------------------------------------------------------------!
997! Description:
998! ------------
999!> This routine writes the respective restart data.
1000!------------------------------------------------------------------------------!
1001 SUBROUTINE stg_wrd_global
1002
1003
1004    IMPLICIT NONE
1005
1006    CALL wrd_write_string( 'mc_factor' )
1007    WRITE ( 14 )  mc_factor
1008
1009    CALL wrd_write_string( 'use_syn_turb_gen' )
1010    WRITE ( 14 )  use_syn_turb_gen
1011
1012
1013 END SUBROUTINE stg_wrd_global
1014
1015
1016!------------------------------------------------------------------------------!
1017! Description:
1018! ------------
1019!> Create turbulent inflow fields for u, v, w with prescribed length scales and
1020!> Reynolds stress tensor after a method of Xie and Castro (2008), modified
1021!> following suggestions of Kim et al. (2013), and using a Lund rotation
1022!> (Lund, 1998).
1023!------------------------------------------------------------------------------!
1024 SUBROUTINE stg_main
1025
1026
1027    USE arrays_3d,                                                             &
1028        ONLY:  dzw
1029
1030    USE control_parameters,                                                    &
1031        ONLY:  child_domain, dt_3d, intermediate_timestep_count,               &
1032               nesting_offline, rans_mode, simulated_time, volume_flow_initial
1033
1034    USE grid_variables,                                                        &
1035        ONLY:  dx, dy
1036
1037    USE indices,                                                               &
1038        ONLY:  wall_flags_0
1039
1040    USE statistics,                                                            &
1041        ONLY:  weight_substep
1042
1043    USE pmc_interface,                                                         &
1044        ONLY : rans_mode_parent
1045
1046
1047    IMPLICIT NONE
1048
1049    INTEGER(iwp) :: i           !< grid index in x-direction
1050    INTEGER(iwp) :: j           !< loop index in y-direction
1051    INTEGER(iwp) :: k           !< loop index in z-direction
1052
1053    REAL(wp) :: dt_stg          !< wheighted subtimestep
1054    REAL(wp) :: mc_factor_l     !< local mass flux correction factor
1055    REAL(wp) :: volume_flow     !< mass flux through lateral boundary
1056    REAL(wp) :: volume_flow_l   !< local mass flux through lateral boundary
1057
1058    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1059
1060!
1061!-- Calculate time step which is needed for filter functions
1062    dt_stg = MAX( dt_3d, dt_stg_call ) !* weight_substep(intermediate_timestep_count)
1063!
1064!-- Initial value of fu, fv, fw
1065    IF ( simulated_time == 0.0_wp .AND. .NOT. velocity_seed_initialized )  THEN
1066       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
1067       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
1068       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
1069
1070       IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent      &
1071                                .AND.  .NOT.  rans_mode ) )  THEN
1072!
1073!--       Generate turbulence at right boundary
1074          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_right )
1075          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_right )
1076          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_right )
1077!
1078!--       Generate turbulence at north boundary
1079          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_north )
1080          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_north )
1081          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_north )
1082!
1083!--       Generate turbulence at south boundary
1084          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_south )
1085          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_south )
1086          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_south )
1087       ENDIF
1088       velocity_seed_initialized = .TRUE.
1089    ENDIF
1090   
1091!
1092!-- New set of fu, fv, fw
1093    CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
1094    CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
1095    CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
1096   
1097    IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
1098                          .AND.  .NOT.  rans_mode ) )  THEN
1099!
1100!--    Generate turbulence at right boundary
1101       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
1102       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
1103       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
1104!
1105!--    Generate turbulence at north boundary
1106       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
1107       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
1108       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
1109!
1110!--    Generate turbulence at south boundary
1111       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
1112       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
1113       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
1114    ENDIF
1115   
1116!
1117!-- Turbulence generation at left and or right boundary
1118    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
1119   
1120       DO  j = nysg, nyng
1121          DO  k = nzb, nzt + 1
1122!
1123!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1124             IF ( tu(k) == 0.0_wp )  THEN
1125                fu_yz(k,j) = fuo_yz(k,j)
1126             ELSE
1127                fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
1128                         fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
1129             ENDIF
1130
1131             IF ( tv(k) == 0.0_wp )  THEN
1132                fv_yz(k,j) = fvo_yz(k,j)
1133             ELSE
1134                fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
1135                         fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
1136             ENDIF
1137
1138             IF ( tw(k) == 0.0_wp )  THEN
1139                fw_yz(k,j) = fwo_yz(k,j)
1140             ELSE
1141                fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
1142                         fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
1143             ENDIF
1144!
1145!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
1146!--          Additional factors are added to improve the variance of v and w
1147             IF( k == 0 )  THEN
1148                dist_yz(k,j,1) = 0.0_wp
1149                dist_yz(k,j,2) = 0.0_wp
1150                dist_yz(k,j,3) = 0.0_wp
1151             ELSE
1152                dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp )
1153                !experimental test of 1.2
1154                dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
1155                                      * 1.2_wp )                               &
1156                                     * (   a21(k) * fu_yz(k,j)                 &
1157                                         + a22(k) * fv_yz(k,j) ), 3.0_wp )
1158                dist_yz(k,j,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
1159                                      * 1.3_wp )                               &
1160                                    * (   a31(k) * fu_yz(k,j)                  &
1161                                        + a32(k) * fv_yz(k,j)                  &
1162                                        + a33(k) * fw_yz(k,j) ), 3.0_wp )
1163             ENDIF
1164
1165          ENDDO
1166       ENDDO
1167
1168!
1169!--    Mass flux correction following Kim et al. (2013)
1170!--    This correction factor insures that the mass flux is preserved at the
1171!--    inflow boundary
1172       IF ( .NOT. nesting_offline  .AND.  .NOT. child_domain )  THEN
1173          mc_factor_l = 0.0_wp
1174          mc_factor   = 0.0_wp
1175          DO  j = nys, nyn
1176             DO  k = nzb+1, nzt
1177                mc_factor_l = mc_factor_l + dzw(k)  *                          &
1178                              ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) )
1179             ENDDO
1180          ENDDO
1181
1182#if defined( __parallel )
1183          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 1, MPI_REAL, MPI_SUM,    &
1184                              comm1dy, ierr )
1185#else
1186          mc_factor = mc_factor_l
1187#endif
1188
1189          mc_factor = volume_flow_initial(1) / mc_factor
1190
1191!
1192!--       Add disturbance at the inflow
1193          DO  j = nysg, nyng
1194             DO  k = nzb, nzt+1
1195                 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +              &
1196                                      dist_yz(k,j,1)             ) * mc_factor
1197                 v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +              &
1198                                      dist_yz(k,j,2)             ) * mc_factor
1199                 w(k,j,-nbgp:-1)  =   dist_yz(k,j,3)               * mc_factor
1200             ENDDO
1201          ENDDO
1202
1203       ELSE
1204!
1205!--       First, calculate volume flow at yz boundary
1206          IF ( myidx == id_stg_left  )  i = nxl
1207          IF ( myidx == id_stg_right )  i = nxr+1
1208
1209          volume_flow_l = 0.0_wp
1210          volume_flow   = 0.0_wp
1211          mc_factor_l   = 0.0_wp
1212          mc_factor     = 0.0_wp
1213          DO  j = nys, nyn
1214             DO  k = nzb+1, nzt
1215                volume_flow_l = volume_flow_l + u(k,j,i) * dzw(k) * dy         &
1216                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1217                                              BTEST( wall_flags_0(k,j,i), 1 ) )
1218
1219                mc_factor_l = mc_factor_l     + ( u(k,j,i) + dist_yz(k,j,1) )  &
1220                                                         * dzw(k) * dy         &
1221                                     * MERGE( 1.0_wp, 0.0_wp,                  &
1222                                              BTEST( wall_flags_0(k,j,i), 1 ) )
1223             ENDDO
1224          ENDDO
1225#if defined( __parallel )
1226          CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                      &
1227                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1228          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                          &
1229                              1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1230#else
1231          volume_flow = volume_flow_l
1232          mc_factor   = mc_factor_l
1233#endif
1234
1235          mc_factor = volume_flow / mc_factor
1236
1237!
1238!--       Add disturbances
1239          IF ( myidx == id_stg_left  )  THEN
1240             DO  j = nys, nyn
1241                DO  k = nzb+1, nzt
1242                   u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) )                    &
1243                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
1244                                            BTEST( wall_flags_0(k,j,0), 1 ) )
1245                   u(k,j,-1) = u(k,j,0)
1246                   v(k,j,-1)  = ( v(k,j,-1)  + dist_yz(k,j,2)  )               &
1247                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
1248                                            BTEST( wall_flags_0(k,j,-1), 2 ) )
1249                   w(k,j,-1)  = ( w(k,j,-1)  + dist_yz(k,j,3)  )               &
1250                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
1251                                            BTEST( wall_flags_0(k,j,-1), 3 ) )
1252                ENDDO
1253             ENDDO
1254          ENDIF
1255          IF ( myidx == id_stg_right  )  THEN
1256             DO  j = nys, nyn
1257                DO  k = nzb+1, nzt
1258                   u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) )            &
1259                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1260                                          BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
1261                   v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) )            &
1262                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1263                                          BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
1264                   w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) )            &
1265                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1266                                          BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
1267                ENDDO
1268             ENDDO
1269          ENDIF
1270       ENDIF
1271
1272    ENDIF
1273!
1274!-- Turbulence generation at north and south boundary
1275    IF ( myidy == id_stg_north  .OR.  myidy == id_stg_south )  THEN
1276   
1277       DO  i = nxlg, nxrg
1278          DO  k = nzb, nzt + 1
1279!
1280!--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1281             IF ( tu(k) == 0.0_wp )  THEN
1282                fu_xz(k,i) = fuo_xz(k,i)
1283             ELSE
1284                fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
1285                         fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
1286             ENDIF
1287
1288             IF ( tv(k) == 0.0_wp )  THEN
1289                fv_xz(k,i) = fvo_xz(k,i)
1290             ELSE
1291                fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
1292                      fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
1293             ENDIF
1294           
1295             IF ( tw(k) == 0.0_wp )  THEN
1296                fw_xz(k,i) = fwo_xz(k,i)
1297             ELSE
1298                fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
1299                      fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
1300             ENDIF
1301!
1302!--          Lund rotation following Eq. 17 in Xie and Castro (2008).
1303!--          Additional factors are added to improve the variance of v and w
1304             IF( k == 0 )  THEN
1305                dist_xz(k,i,1) = 0.0_wp
1306                dist_xz(k,i,2) = 0.0_wp
1307                dist_xz(k,i,3) = 0.0_wp
1308
1309             ELSE
1310                dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp )
1311                !experimental test of 1.2
1312                dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
1313                                      * 1.2_wp )                               &
1314                                     * (   a21(k) * fu_xz(k,i)                 &
1315                                         + a22(k) * fv_xz(k,i) ), 3.0_wp )
1316                dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
1317                                      * 1.3_wp )                               &
1318                                    * (   a31(k) * fu_xz(k,i)                  &
1319                                        + a32(k) * fv_xz(k,i)                  &
1320                                        + a33(k) * fw_xz(k,i) ), 3.0_wp )
1321             ENDIF
1322
1323          ENDDO
1324       ENDDO
1325
1326!
1327!--    Mass flux correction following Kim et al. (2013)
1328!--    This correction factor insures that the mass flux is preserved at the
1329!--    inflow boundary.
1330!--    First, calculate volume flow at xz boundary
1331       IF ( myidy == id_stg_south  ) j = nys
1332       IF ( myidy == id_stg_north )  j = nyn+1
1333
1334       volume_flow_l = 0.0_wp
1335       volume_flow   = 0.0_wp
1336       mc_factor_l   = 0.0_wp
1337       mc_factor     = 0.0_wp
1338       DO  i = nxl, nxr
1339          DO  k = nzb+1, nzt
1340             volume_flow_l = volume_flow_l + v(k,j,i) * dzw(k) * dx            &
1341                                  * MERGE( 1.0_wp, 0.0_wp,                     &
1342                                           BTEST( wall_flags_0(k,j,i), 2 ) )
1343
1344             mc_factor_l = mc_factor_l     + ( v(k,j,i) + dist_xz(k,i,2) )     &
1345                                                      * dzw(k) * dx            &
1346                                  * MERGE( 1.0_wp, 0.0_wp,                     &
1347                                           BTEST( wall_flags_0(k,j,i), 2 ) )
1348          ENDDO
1349       ENDDO
1350#if defined( __parallel )
1351       CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                         &
1352                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1353       CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                             &
1354                           1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1355#else
1356       volume_flow = volume_flow_l
1357       mc_factor   = mc_factor_l
1358#endif
1359
1360       mc_factor = volume_flow / mc_factor
1361
1362!
1363!--    Add disturbances
1364       IF ( myidy == id_stg_south  )  THEN
1365
1366          DO  i = nxl, nxr
1367             DO  k = nzb+1, nzt
1368                u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) )                     &
1369                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
1370                                           BTEST( wall_flags_0(k,-1,i), 1 ) )
1371                v(k,0,i)  = ( v(k,0,i)  + dist_xz(k,i,2)  )                    &
1372                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
1373                                           BTEST( wall_flags_0(k,0,i), 2 ) )
1374                v(k,-1,i) = v(k,0,i)
1375                w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3)  )                    &
1376                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
1377                                           BTEST( wall_flags_0(k,-1,i), 3 ) )
1378             ENDDO
1379          ENDDO
1380       ENDIF
1381       IF ( myidy == id_stg_north  )  THEN
1382
1383          DO  i = nxl, nxr
1384             DO  k = nzb+1, nzt
1385                u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )               &
1386                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1387                                          BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
1388                v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) )               &
1389                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1390                                          BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
1391                w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) )               &
1392                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
1393                                          BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
1394             ENDDO
1395          ENDDO
1396       ENDIF
1397    ENDIF
1398!
1399!-- Finally, set time counter for calling STG to zero
1400    time_stg_call = 0.0_wp
1401
1402    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1403
1404 END SUBROUTINE stg_main
1405
1406!------------------------------------------------------------------------------!
1407! Description:
1408! ------------
1409!> Generate a set of random number rand_it wich is equal on each PE
1410!> and calculate the velocity seed f_n.
1411!> f_n is splitted in vertical direction by the number of PEs in x-direction and
1412!> and each PE calculates a vertical subsection of f_n. At the the end, all
1413!> parts are collected to form the full array.
1414!------------------------------------------------------------------------------!
1415 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
1416
1417
1418    USE indices,                                                               &
1419        ONLY: ny
1420
1421    IMPLICIT NONE
1422
1423    INTEGER(iwp) :: id          !< core ids at respective boundaries
1424    INTEGER(iwp) :: j           !< loop index in y-direction
1425    INTEGER(iwp) :: jj          !< loop index in y-direction
1426    INTEGER(iwp) :: k           !< loop index in z-direction
1427    INTEGER(iwp) :: kk          !< loop index in z-direction
1428    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
1429
1430    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y    !< length scale in y-direction
1431    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
1432    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y2   !< n_y*2
1433    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
1434
1435    REAL(wp) :: nyz_inv         !< inverse of number of grid points in yz-slice
1436    REAL(wp) :: rand_av         !< average of random number
1437    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
1438
1439    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter function in y-direction
1440    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
1441    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
1442    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)     :: f_n     !<  velocity seed
1443    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
1444
1445
1446!
1447!-- Generate random numbers using a seed generated in stg_init.
1448!-- The set of random numbers are modified to have an average of 0 and
1449!-- unit variance.
1450    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:ny+mergp) )
1451
1452    rand_av        = 0.0_wp
1453    rand_sigma_inv = 0.0_wp
1454    nyz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( ny+1 ), KIND=wp )
1455
1456    DO  j = 0, ny
1457       DO  k = nzb, nzt+1
1458          CALL RANDOM_NUMBER( rand_it(k,j) )
1459          rand_av = rand_av + rand_it(k,j)
1460       ENDDO
1461    ENDDO
1462
1463    rand_av = rand_av * nyz_inv
1464
1465    DO  j = 0, ny
1466       DO  k = nzb, nzt+1
1467          rand_it(k,j)   = rand_it(k,j) - rand_av
1468          rand_sigma_inv = rand_sigma_inv + rand_it(k,j) ** 2
1469       ENDDO
1470    ENDDO
1471
1472    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nyz_inv)
1473
1474    DO  j = 0, ny
1475       DO  k = nzb, nzt+1
1476          rand_it(k,j) = rand_it(k,j) * rand_sigma_inv
1477       ENDDO
1478    ENDDO
1479
1480!
1481!-- Periodic fill of random number in space
1482    DO  j = 0, ny
1483       DO  k = 1, mergp
1484          rand_it(nzb  -k,j) = rand_it(nzt+2-k,j)    ! bottom margin
1485          rand_it(nzt+1+k,j) = rand_it(nzb+k-1,j)    ! top margin
1486       ENDDO
1487    ENDDO
1488    DO  j = 1, mergp
1489       DO  k = nzb-mergp, nzt+1+mergp
1490          rand_it(k,  -j) = rand_it(k,ny-j+1)        ! south margin
1491          rand_it(k,ny+j) = rand_it(k,   j-1)        ! north margin
1492       ENDDO
1493    ENDDO
1494
1495!
1496!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1497    n_y2 = n_y * 2
1498    n_z2 = n_z * 2
1499    f_n_l  = 0.0_wp
1500
1501    DO  j = nysg, nyng
1502       DO  k = nzb_x_stg, nzt_x_stg+1
1503          DO  jj = -n_y2(k), n_y2(k)
1504             DO  kk = -n_z2(k), n_z2(k)
1505                f_n_l(k,j) = f_n_l(k,j)                                        &
1506                           + b_y(jj,k) * b_z(kk,k) * rand_it(k+kk,j+jj)
1507             ENDDO
1508          ENDDO
1509       ENDDO
1510    ENDDO
1511
1512    DEALLOCATE( rand_it )
1513!
1514!-- Gather velocity seeds of full subdomain
1515    send_count = nzt_x_stg - nzb_x_stg + 1
1516    IF ( nzt_x_stg == nzt )  send_count = send_count + 1
1517
1518#if defined( __parallel )
1519    CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,    &
1520                      f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,  &
1521                      id, comm1dx, ierr )
1522#else
1523    f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x_stg:nzt_x_stg+1,nysg:nyng)
1524#endif
1525
1526
1527 END SUBROUTINE stg_generate_seed_yz
1528
1529
1530
1531
1532!------------------------------------------------------------------------------!
1533! Description:
1534! ------------
1535!> Generate a set of random number rand_it wich is equal on each PE
1536!> and calculate the velocity seed f_n.
1537!> f_n is splitted in vertical direction by the number of PEs in y-direction and
1538!> and each PE calculates a vertical subsection of f_n. At the the end, all
1539!> parts are collected to form the full array.
1540!------------------------------------------------------------------------------!
1541 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
1542
1543
1544    USE indices,                                                               &
1545        ONLY: nx
1546
1547
1548    IMPLICIT NONE
1549
1550    INTEGER(iwp) :: id          !< core ids at respective boundaries
1551    INTEGER(iwp) :: i           !< loop index in x-direction
1552    INTEGER(iwp) :: ii          !< loop index in x-direction
1553    INTEGER(iwp) :: k           !< loop index in z-direction
1554    INTEGER(iwp) :: kk          !< loop index in z-direction
1555    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
1556
1557    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x    !< length scale in x-direction
1558    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
1559    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_x*2
1560    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
1561
1562    REAL(wp) :: nxz_inv         !< inverse of number of grid points in xz-slice
1563    REAL(wp) :: rand_av         !< average of random number
1564    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
1565
1566    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter function in x-direction
1567    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
1568    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
1569    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)     :: f_n     !<  velocity seed
1570    REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
1571
1572
1573!
1574!-- Generate random numbers using a seed generated in stg_init.
1575!-- The set of random numbers are modified to have an average of 0 and
1576!-- unit variance.
1577    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:nx+mergp) )
1578
1579    rand_av        = 0.0_wp
1580    rand_sigma_inv = 0.0_wp
1581    nxz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( nx+1 ), KIND=wp )
1582
1583    DO  i = 0, nx
1584       DO  k = nzb, nzt+1
1585          CALL RANDOM_NUMBER( rand_it(k,i) )
1586          rand_av = rand_av + rand_it(k,i)
1587       ENDDO
1588    ENDDO
1589
1590    rand_av = rand_av * nxz_inv
1591
1592    DO  i = 0, nx
1593       DO  k = nzb, nzt+1
1594          rand_it(k,i)   = rand_it(k,i) - rand_av
1595          rand_sigma_inv = rand_sigma_inv + rand_it(k,i) ** 2
1596       ENDDO
1597    ENDDO
1598
1599    rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nxz_inv)
1600
1601    DO  i = 0, nx
1602       DO  k = nzb, nzt+1
1603          rand_it(k,i) = rand_it(k,i) * rand_sigma_inv
1604       ENDDO
1605    ENDDO
1606
1607!
1608!-- Periodic fill of random number in space
1609    DO  i = 0, nx
1610       DO  k = 1, mergp
1611          rand_it(nzb-k,i)   = rand_it(nzt+2-k,i)    ! bottom margin
1612          rand_it(nzt+1+k,i) = rand_it(nzb+k-1,i)    ! top margin
1613       ENDDO
1614    ENDDO
1615    DO  i = 1, mergp
1616       DO  k = nzb-mergp, nzt+1+mergp
1617          rand_it(k,-i)   = rand_it(k,nx-i+1)        ! left margin
1618          rand_it(k,nx+i) = rand_it(k,i-1)           ! right margin
1619       ENDDO
1620    ENDDO
1621
1622!
1623!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1624    n_x2 = n_x * 2
1625    n_z2 = n_z * 2
1626    f_n_l  = 0.0_wp
1627
1628    DO  i = nxlg, nxrg
1629       DO  k = nzb_y_stg, nzt_y_stg+1
1630          DO  ii = -n_x2(k), n_x2(k)
1631             DO  kk = -n_z2(k), n_z2(k)
1632                f_n_l(k,i) = f_n_l(k,i)                                        &
1633                           + b_x(ii,k) * b_z(kk,k) * rand_it(k+kk,i+ii)
1634             ENDDO
1635          ENDDO
1636       ENDDO
1637    ENDDO
1638
1639    DEALLOCATE( rand_it )
1640
1641!
1642!-- Gather velocity seeds of full subdomain
1643    send_count = nzt_y_stg - nzb_y_stg + 1
1644    IF ( nzt_y_stg == nzt )  send_count = send_count + 1
1645
1646
1647#if defined( __parallel )
1648    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,    &
1649                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz,  &
1650                      id, comm1dy, ierr )
1651#else
1652    f_n(nzb+1:nzt+1,nxlg:nxrg) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg)
1653#endif
1654
1655
1656 END SUBROUTINE stg_generate_seed_xz
1657
1658!------------------------------------------------------------------------------!
1659! Description:
1660! ------------
1661!> Parametrization of the Reynolds stress tensor, following the parametrization
1662!> described in Rotach et al. (1996), which is applied in state-of-the-art
1663!> dispserion modelling. Please note, the parametrization does not distinguish
1664!> between along-wind and cross-wind turbulence.
1665!------------------------------------------------------------------------------!
1666 SUBROUTINE parametrize_reynolds_stress
1667
1668    USE arrays_3d,                                                             &
1669        ONLY:  zu
1670
1671    IMPLICIT NONE
1672
1673    INTEGER(iwp) :: k            !< loop index in z-direction
1674   
1675    REAL(wp)     ::  zzi         !< ratio of z/zi
1676   
1677!
1678!--
1679    DO  k = nzb+1, nzt+1
1680         
1681       IF ( zu(k) <= zi_ribulk )  THEN
1682!
1683!--       Determine normalized height coordinate
1684          zzi = zu(k) / zi_ribulk
1685!
1686!--       u'u' and v'v'. Assume isotropy. Note, add a small negative number
1687!--       to the denominator, else the merge-function can crash if scale_l is
1688!--       zero.
1689          r11(k) = scale_us**2 * (                                             &
1690                      MERGE( 0.35_wp * (                                       &
1691                           - zi_ribulk / ( kappa * scale_l - 10E-4_wp )        &
1692                                       )**( 2.0_wp / 3.0_wp ),                 &
1693                             0.0_wp,                                           &
1694                             scale_l < 0.0_wp )                                &
1695                    + 5.0_wp - 4.0_wp * zzi                                    &
1696                                 )
1697                                 
1698          r22(k) = r11(k)
1699!
1700!--       w'w'
1701          r33(k) = scale_wm**2 * (                                             &
1702                      1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi ) &
1703                    + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2            &                     
1704                                 )
1705!
1706!--       u'w' and v'w'. Assume isotropy.
1707          r31(k) = - scale_us**2 * (                                           &
1708                         1.0_wp - EXP( 3.0_wp * ( zzi - 1.0_wp ) )             &
1709                                   )
1710                             
1711          r32(k) = r31(k)
1712!
1713!--       For u'v' no parametrization exist so far - ?. For simplicity assume
1714!--       a similar profile as for u'w'.
1715          r21(k) = r31(k)
1716!
1717!--    Above the boundary layer, assmume laminar flow conditions.
1718       ELSE
1719          r11(k) = 10E-8_wp
1720          r22(k) = 10E-8_wp
1721          r33(k) = 10E-8_wp
1722          r21(k) = 10E-8_wp
1723          r31(k) = 10E-8_wp
1724          r32(k) = 10E-8_wp
1725       ENDIF
1726!        write(9,*) zu(k), r11(k), r33(k), r31(k), zi_ribulk, scale_us, scale_wm, scale_l
1727    ENDDO
1728
1729!
1730!-- Set bottom boundary condition   
1731    r11(nzb) = r11(nzb+1)
1732    r22(nzb) = r22(nzb+1)
1733    r33(nzb) = r33(nzb+1)
1734
1735    r21(nzb) = r11(nzb+1)
1736    r31(nzb) = r31(nzb+1)
1737    r32(nzb) = r32(nzb+1)
1738   
1739
1740 END SUBROUTINE parametrize_reynolds_stress 
1741 
1742!------------------------------------------------------------------------------!
1743! Description:
1744! ------------
1745!> Calculate the coefficient matrix from the Lund rotation.
1746!------------------------------------------------------------------------------!
1747 SUBROUTINE calc_coeff_matrix
1748
1749    IMPLICIT NONE
1750
1751    INTEGER(iwp) :: k   !< loop index in z-direction
1752   
1753!
1754!-- Calculate coefficient matrix. Split loops to allow for loop vectorization.
1755    DO  k = nzb+1, nzt+1
1756       IF ( r11(k) > 0.0_wp )  THEN
1757          a11(k) = SQRT( r11(k) ) 
1758          a21(k) = r21(k) / a11(k)
1759          a31(k) = r31(k) / a11(k)
1760       ELSE
1761          a11(k) = 10E-8_wp
1762          a21(k) = 10E-8_wp
1763          a31(k) = 10E-8_wp
1764       ENDIF
1765    ENDDO
1766    DO  k = nzb+1, nzt+1
1767       a22(k) = r22(k) - a21(k)**2 
1768       IF ( a22(k) > 0.0_wp )  THEN
1769          a22(k) = SQRT( a22(k) )
1770          a32(k) = r32(k) - a21(k) * a31(k) / a22(k)
1771       ELSE
1772          a22(k) = 10E-8_wp
1773          a32(k) = 10E-8_wp
1774       ENDIF
1775    ENDDO
1776    DO  k = nzb+1, nzt+1
1777       a33(k) = r33(k) - a31(k)**2 - a32(k)**2
1778       IF ( a33(k) > 0.0_wp )  THEN
1779          a33(k) =  SQRT( a33(k) )
1780       ELSE
1781          a33(k) = 10E-8_wp
1782       ENDIF
1783    ENDDO   
1784!
1785!-- Set bottom boundary condition
1786    a11(nzb) = a11(nzb+1)
1787    a22(nzb) = a22(nzb+1)
1788    a21(nzb) = a21(nzb+1)
1789    a33(nzb) = a33(nzb+1)   
1790    a31(nzb) = a31(nzb+1)
1791    a32(nzb) = a32(nzb+1)   
1792
1793 END SUBROUTINE calc_coeff_matrix
1794 
1795!------------------------------------------------------------------------------!
1796! Description:
1797! ------------
1798!> This routine controls the re-adjustment of the turbulence statistics used
1799!> for generating turbulence at the lateral boundaries. 
1800!------------------------------------------------------------------------------!
1801 SUBROUTINE stg_adjust
1802
1803    IMPLICIT NONE
1804   
1805    INTEGER(iwp) ::  k
1806
1807    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1808!
1809!-- Compute mean boundary layer height according to Richardson-Bulk
1810!-- criterion using the inflow profiles. Further velocity scale as well as
1811!-- mean friction velocity are calculated.
1812    CALL calc_scaling_variables
1813!
1814!-- Set length and time scales depending on boundary-layer height
1815    CALL calc_length_and_time_scale
1816!
1817!-- Parametrize Reynolds-stress tensor, diagonal elements as well
1818!-- as r21 (v'u'), r31 (w'u'), r32 (w'v'). Parametrization follows
1819!-- Rotach et al. (1996) and is based on boundary-layer depth,
1820!-- friction velocity and velocity scale.
1821    CALL parametrize_reynolds_stress
1822!
1823!-- Calculate coefficient matrix from Reynolds stress tensor 
1824!-- (Lund rotation)
1825    CALL calc_coeff_matrix
1826!
1827!-- Determine filter functions on basis of updated length scales
1828    CALL stg_filter_func( nux, bux ) !filter ux
1829    CALL stg_filter_func( nuy, buy ) !filter uy
1830    CALL stg_filter_func( nuz, buz ) !filter uz
1831    CALL stg_filter_func( nvx, bvx ) !filter vx
1832    CALL stg_filter_func( nvy, bvy ) !filter vy
1833    CALL stg_filter_func( nvz, bvz ) !filter vz
1834    CALL stg_filter_func( nwx, bwx ) !filter wx
1835    CALL stg_filter_func( nwy, bwy ) !filter wy
1836    CALL stg_filter_func( nwz, bwz ) !filter wz
1837!
1838!-- Reset time counter for controlling next adjustment to zero
1839    time_stg_adjust = 0.0_wp
1840   
1841    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1842   
1843 END SUBROUTINE stg_adjust 
1844 
1845 
1846!------------------------------------------------------------------------------!
1847! Description:
1848! ------------
1849!> Calculates turbuluent length and time scales if these are not available
1850!> from measurements.
1851!------------------------------------------------------------------------------!
1852 SUBROUTINE calc_length_and_time_scale
1853
1854    USE arrays_3d,                                                             &
1855        ONLY:  dzw, ddzw, u_init, v_init, zu, zw
1856
1857    USE grid_variables,                                                        &
1858        ONLY:  ddx, ddy, dx, dy
1859       
1860    USE statistics,                                                            &
1861        ONLY:  hom
1862
1863    IMPLICIT NONE
1864
1865
1866    INTEGER(iwp) :: k            !< loop index in z-direction
1867
1868    REAL(wp) ::  length_scale    !< typical length scale
1869   
1870!
1871!-- In initial call the boundary-layer depth can be zero. This case, set
1872!-- minimum value for boundary-layer depth, to setup length scales correctly.
1873    zi_ribulk = MAX( zi_ribulk, zw(nzb+2) )
1874!
1875!-- Set-up default turbulent length scales. From the numerical point of
1876!-- view the imposed perturbations should not be immediately dissipated
1877!-- by the numerics. The numerical dissipation, however, acts on scales
1878!-- up to 8 x the grid spacing. For this reason, set the turbulence
1879!-- length scale to 8 time the grid spacing. Further, above the boundary
1880!-- layer height, set turbulence lenght scales to zero (equivalent to not
1881!-- imposing any perturbations) in order to save computational costs.
1882!-- Typical time scales are derived by assuming Taylors's hypothesis,
1883!-- using the length scales and the mean profiles of u- and v-component.
1884    DO  k = nzb+1, nzt+1
1885   
1886       length_scale = 8.0_wp * MIN( dx, dy, dzw(k) )
1887         
1888       IF ( zu(k) <= zi_ribulk )  THEN 
1889!
1890!--       Assume isotropic turbulence length scales
1891          nux(k) = MAX( INT( length_scale * ddx     ), 1 )
1892          nuy(k) = MAX( INT( length_scale * ddy     ), 1 )
1893          nuz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
1894          nvx(k) = MAX( INT( length_scale * ddx     ), 1 )
1895          nvy(k) = MAX( INT( length_scale * ddy     ), 1 )
1896          nvz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
1897          nwx(k) = MAX( INT( length_scale * ddx     ), 1 )
1898          nwy(k) = MAX( INT( length_scale * ddy     ), 1 )
1899          nwz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
1900!
1901!--       Limit time scales, else they become very larger for low wind speed,
1902!--       imposing long-living inflow perturbations which in turn propagate
1903!--       further into the model domain. Use u_init and v_init to calculate
1904!--       the time scales, which will be equal to the inflow profiles, both,
1905!--       in offline nesting mode or in dirichlet/radiation mode.
1906          tu(k)  = MIN( dt_stg_adjust, length_scale /                          &
1907                                  ( ABS( u_init(k) ) + 0.1_wp ) )
1908          tv(k)  = MIN( dt_stg_adjust, length_scale /                          &
1909                                  ( ABS( v_init(k) ) + 0.1_wp ) )
1910!
1911!--       Time scale of w-component is a mixture from u- and v-component.
1912          tw(k)  = SQRT( tu(k)**2 + tv(k)**2 )
1913!
1914!--    Above the boundary layer length scales are zero, i.e. imposed turbulence
1915!--    is not correlated in space and time, just white noise. This saves
1916!--    computations power.
1917       ELSE
1918          nux(k) = 0.0_wp
1919          nuy(k) = 0.0_wp
1920          nuz(k) = 0.0_wp
1921          nvx(k) = 0.0_wp
1922          nvy(k) = 0.0_wp
1923          nvz(k) = 0.0_wp
1924          nwx(k) = 0.0_wp
1925          nwy(k) = 0.0_wp
1926          nwz(k) = 0.0_wp
1927         
1928          tu(k)  = 0.0_wp
1929          tv(k)  = 0.0_wp
1930          tw(k)  = 0.0_wp
1931       ENDIF
1932    ENDDO
1933!
1934!-- Set bottom boundary condition for the length and time scales
1935    nux(nzb) = nux(nzb+1)
1936    nuy(nzb) = nuy(nzb+1)
1937    nuz(nzb) = nuz(nzb+1)
1938    nvx(nzb) = nvx(nzb+1)
1939    nvy(nzb) = nvy(nzb+1)
1940    nvz(nzb) = nvz(nzb+1)
1941    nwx(nzb) = nwx(nzb+1)
1942    nwy(nzb) = nwy(nzb+1)
1943    nwz(nzb) = nwz(nzb+1)
1944
1945    tu(nzb)  = tu(nzb+1)
1946    tv(nzb)  = tv(nzb+1)
1947    tw(nzb)  = tw(nzb+1)
1948
1949
1950 END SUBROUTINE calc_length_and_time_scale 
1951
1952!------------------------------------------------------------------------------!
1953! Description:
1954! ------------
1955!> Calculate scaling variables which are used for turbulence parametrization
1956!> according to Rotach et al. (1996). Scaling variables are: friction velocity,
1957!> boundary-layer depth, momentum velocity scale, and Obukhov length.
1958!------------------------------------------------------------------------------!
1959 SUBROUTINE calc_scaling_variables
1960
1961
1962    USE control_parameters,                                                    &
1963        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
1964               humidity, pt_surface
1965             
1966    USE indices,                                                               &
1967        ONLY:  nx, ny 
1968       
1969    USE surface_mod,                                                           &
1970        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
1971
1972    IMPLICIT NONE
1973
1974    INTEGER(iwp) :: i            !< loop index in x-direction
1975    INTEGER(iwp) :: j            !< loop index in y-direction
1976    INTEGER(iwp) :: k            !< loop index in z-direction
1977    INTEGER(iwp) :: k_ref        !< index in z-direction for reference height
1978    INTEGER(iwp) :: m            !< surface element index
1979
1980    REAL(wp) ::  friction_vel_l         !< mean friction veloctiy on subdomain
1981    REAL(wp) ::  shf_mean               !< mean surface sensible heat flux
1982    REAL(wp) ::  shf_mean_l             !< mean surface sensible heat flux on subdomain
1983    REAL(wp) ::  u_int                  !< u-component
1984    REAL(wp) ::  v_int                  !< v-component
1985    REAL(wp) ::  vpt_surface            !< near-surface virtual potential temperature
1986    REAL(wp) ::  w_convective           !< convective velocity scale
1987    REAL(wp) ::  z0_mean                !< mean roughness length
1988    REAL(wp) ::  z0_mean_l              !< mean roughness length on subdomain
1989   
1990!
1991!-- Mean friction velocity and velocity scale. Therefore,
1992!-- pre-calculate mean roughness length and surface sensible heat flux
1993!-- in the model domain, which are further used to estimate friction
1994!-- velocity and velocity scale. Note, for z0 linear averaging is applied,
1995!-- even though this is known to unestimate the effective roughness.
1996!-- This need to be revised later.
1997    z0_mean_l  = 0.0_wp
1998    shf_mean_l = 0.0_wp
1999    DO  m = 1, surf_def_h(0)%ns
2000       z0_mean_l  = z0_mean_l  + surf_def_h(0)%z0(m)
2001       shf_mean_l = shf_mean_l + surf_def_h(0)%shf(m)
2002    ENDDO   
2003    DO  m = 1, surf_lsm_h%ns
2004       z0_mean_l  = z0_mean_l  + surf_lsm_h%z0(m)
2005       shf_mean_l = shf_mean_l + surf_lsm_h%shf(m)
2006    ENDDO
2007    DO  m = 1, surf_usm_h%ns
2008       z0_mean_l  = z0_mean_l  + surf_usm_h%z0(m)
2009       shf_mean_l = shf_mean_l + surf_usm_h%shf(m)
2010    ENDDO
2011   
2012#if defined( __parallel )
2013    CALL MPI_ALLREDUCE( z0_mean_l, z0_mean, 1, MPI_REAL, MPI_SUM,              &
2014                        comm2d, ierr )
2015    CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM,            &
2016                        comm2d, ierr )
2017#else
2018    z0_mean  = z0_mean_l
2019    shf_mean = shf_mean_l
2020#endif
2021    z0_mean  = z0_mean  / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
2022    shf_mean = shf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
2023       
2024!
2025!-- Note, Inifor does not use logarithmic interpolation of the
2026!-- velocity components near the ground, so that near-surface
2027!-- wind speed and thus the friction velocity is overestimated.
2028!-- However, friction velocity is used for turbulence
2029!-- parametrization, so that more physically meaningful values are important.
2030!-- Hence, derive friction velocity from wind speed at a reference height.
2031!-- For a first guess use 20 m, which is in the range of the first
2032!-- COSMO vertical level.
2033    k_ref = MINLOC( ABS( zu - cosmo_ref ), DIM = 1 ) - 1
2034   
2035    friction_vel_l = 0.0_wp
2036    IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
2037       
2038       i = MERGE( -1, nxr + 1, bc_dirichlet_l )
2039
2040       DO  j = nys, nyn
2041!
2042!--       Note, in u- and v- component the imposed perturbations
2043!--       from the STG are already included. Check whether this
2044!--       makes any difference compared to using the pure-mean
2045!--       inflow profiles.
2046          u_int = MERGE( u(k_ref,j,i+1), u(k_ref,j,i), bc_dirichlet_l )
2047          v_int = v(k_ref,j,i)
2048!
2049!--       Calculate friction velocity and sum-up. Therefore, assume
2050!--       neutral condtions.
2051          friction_vel_l = friction_vel_l + kappa *                            &
2052                            SQRT( u_int * u_int + v_int * v_int )  /           &
2053                            LOG( zu(k_ref) / z0_mean )
2054         
2055       ENDDO
2056       
2057    ENDIF
2058   
2059    IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
2060       
2061       j = MERGE( -1, nyn + 1, bc_dirichlet_s )
2062   
2063       DO  i = nxl, nxr
2064         
2065          u_int = u(k_ref,j,i)
2066          v_int = MERGE( v(k_ref,j+1,i), v(k_ref,j,i), bc_dirichlet_s )
2067 
2068          friction_vel_l = friction_vel_l + kappa *                            &
2069                            SQRT( u_int * u_int + v_int * v_int )  /           &
2070                            LOG( zu(k_ref) / z0_mean )                       
2071         
2072       ENDDO
2073       
2074    ENDIF
2075   
2076#if defined( __parallel )
2077    CALL MPI_ALLREDUCE( friction_vel_l, scale_us, 1, MPI_REAL, MPI_SUM,        &
2078                        comm2d, ierr )
2079#else
2080    scale_us = friction_vel_l
2081#endif
2082    scale_us = scale_us / REAL( 2 * nx + 2 * ny, KIND = wp )
2083   
2084!
2085!-- Compute mean Obukhov length
2086    IF ( shf_mean > 0.0_wp )  THEN
2087       scale_l = - pt_surface / ( kappa * g ) * scale_us**3 / shf_mean
2088    ELSE
2089       scale_l = 0.0_wp
2090    ENDIF
2091   
2092!
2093!-- Compute mean convective velocity scale. Note, in case the mean heat flux
2094!-- is negative, set convective velocity scale to zero.
2095    IF ( shf_mean > 0.0_wp )  THEN
2096       w_convective = ( g * shf_mean * zi_ribulk / pt_surface )**( 1.0_wp / 3.0_wp )
2097    ELSE
2098       w_convective = 0.0_wp
2099    ENDIF
2100!
2101!-- Finally, in order to consider also neutral or stable stratification,
2102!-- compute momentum velocity scale from u* and convective velocity scale,
2103!-- according to Rotach et al. (1996).
2104    scale_wm = ( scale_us**3 + 0.6_wp * w_convective**3 )**( 1.0_wp / 3.0_wp )
2105   
2106 END SUBROUTINE calc_scaling_variables
2107
2108 END MODULE synthetic_turbulence_generator_mod
Note: See TracBrowser for help on using the repository browser.