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

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

Bugfix in initialization of synthetic turbulence generator and calculation of turbulence scaling parameters in case of topography; checks for synthetic turbulence generator and offline nesting added

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