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

Last change on this file since 4144 was 4144, checked in by raasch, 5 years ago

relational operators .EQ., .NE., etc. replaced by ==, /=, etc.

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