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

Last change on this file since 3937 was 3937, checked in by suehring, 2 years ago

Move initialization of STG behind initialization of offline nesting; Bugfix in STG in case of very early restart; calculation of scaling parameters used for parametrization of synthetic turbulence profiles improved; in offline nesting, set boundary value at upper-left and upper-south grid points for u- and v-component, respectively

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