Home

Context Navigation

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

Last change on this file since 4562 was 4562, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 105.8 KB

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