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

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