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

Last change on this file since 3994 was 3987, checked in by kanani, 5 years ago

clean up location, debug and error messages

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