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

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

Correct and clean-up cpu_logs, some overlapping counts (chemistry_model_mod, disturb_heatflux, large_scale_forcing_nudging_mod, ocean_mod, palm, prognostic_equations, synthetic_turbulence_generator_mod, time_integration, time_integration_spinup, turbulence_closure_mod)

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