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

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

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