Home

Context Navigation

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

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