Home

Context Navigation

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

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