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

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

Remove unused variables

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