Home

Context Navigation

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

Last change on this file since 4338 was 4335, checked in by suehring, 5 years ago
Plant canopy: implement fix for LAD on building edges also for the vectorl-optimized branch
Property svn:keywords set to `Id`
File size: 106.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |