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

Last change on this file since 3885 was 3885, checked in by kanani, 2 years ago

restructure/add location/debug messages

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