Home

Context Navigation

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

Last change on this file since 4321 was 4309, checked in by suehring, 5 years ago
Synthetic turbulence generator: Computation of velocity seeds optimized. This implies that random numbers are computed now using the parallel random number generator. Random number are now only computed and normalized locally, while distributed over all mpi ranks afterwards, instead of computing random numbers on a global array. urther, the number of calls for the time-consuming velocity-seed generation is reduced - now the left and right, as well as the north and south boundary share the same velocity-seed matrices.
Property svn:keywords set to `Id`
File size: 105.2 KB

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