Home

Context Navigation

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

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