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

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

Bugfix in initialization of turbulence generator in case of restart runs and offline nesting

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