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

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

Minor bugfix in STG

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