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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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