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

Last change on this file since 3830 was 3775, checked in by gronemeier, 5 years ago

removed unused variables

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