Home

Context Navigation

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

Last change on this file since 3274 was 3274, checked in by knoop, 6 years ago
Modularization of all bulk cloud physics code components
Property svn:keywords set to `Id`
File size: 64.6 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 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: synthetic_turbulence_generator_mod.f90 3274 2018-09-24 15:42:55Z knoop $
27	! Modularization of all bulk cloud physics code components
28	!
29	! 3248 2018-09-14 09:42:06Z sward
30	! Minor formating changes
31	!
32	! 3246 2018-09-13 15:14:50Z sward
33	! Added error handling for input namelist via parin_fail_message
34	!
35	! 3241 2018-09-12 15:02:00Z raasch
36	! unused variables removed
37	!
38	! 3186 2018-07-30 17:07:14Z suehring
39	! Mask topography while imposing inflow perturbations at the boundaries; do not
40	! impose perturbations at top boundary as well as ghost points
41	!
42	! 3183 2018-07-27 14:25:55Z suehring
43	! Rename variables and extend error message
44	! Enable geneartor also for stretched grids
45	!
46	! 3182 2018-07-27 13:36:03Z suehring
47	! Error message related to vertical stretching has been added, dz was replaced
48	! by dz(1)
49	!
50	! 3051 2018-05-30 17:43:55Z suehring
51	! Bugfix in calculation of initial Reynolds-stress tensor.
52	!
53	! 3045 2018-05-28 07:55:41Z Giersch
54	! Error messages revised
55	!
56	! 3044 2018-05-25 10:59:41Z gronemeier
57	! Add missing variable descriptions
58	!
59	! 3038 2018-05-24 10:54:00Z gronemeier
60	! updated variable descriptions
61	!
62	! 2967 2018-04-13 11:22:08Z raasch
63	! bugfix: missing parallel cpp-directives added
64	!
65	! 2946 2018-04-04 17:01:23Z suehring
66	! Remove unused module load
67	!
68	! 2945 2018-04-04 16:27:14Z suehring
69	! - Bugfix in parallelization of synthetic turbulence generator in case the
70	! z-dimension is not an integral divisor of the number of processors along
71	! the x- and y-dimension
72	! - Revision in control mimic in case of RAN-LES nesting
73	!
74	! 2938 2018-03-27 15:52:42Z suehring
75	! Apply turbulence generator at all non-cyclic lateral boundaries in case of
76	! realistic Inifor large-scale forcing or RANS-LES nesting
77	!
78	! 2936 2018-03-27 14:49:27Z suehring
79	! variable named found has been introduced for checking if restart data was found,
80	! reading of restart strings has been moved completely to read_restart_data_mod,
81	! redundant skipping function has been removed, stg_read/write_restart_data
82	! have been renamed to stg_r/wrd_global, stg_rrd_global is called in
83	! read_restart_data_mod now, flag syn_turb_gen_prerun and marker *** end stg
84	! *** have been removed (Giersch), strings and their respective lengths are
85	! written out and read now in case of restart runs to get rid of prescribed
86	! character lengths (Giersch), CASE DEFAULT was added if restart data is read
87	!
88	! 2841 2018-02-27 15:02:57Z suehring
89	! Bugfix: wrong placement of include 'mpif.h' corrected
90	!
91	! 2836 2018-02-26 13:40:05Z Giersch
92	! The variables synthetic_turbulence_generator and
93	! use_synthetic_turbulence_generator have been abbreviated + syn_turb_gen_prerun
94	! flag is used to define if module related parameters were outputted as restart
95	! data
96	!
97	! 2716 2017-12-29 16:35:59Z kanani
98	! Corrected "Former revisions" section
99	!
100	! 2696 2017-12-14 17:12:51Z kanani
101	! Change in file header (GPL part)
102	!
103	! 2669 2017-12-06 16:03:27Z raasch
104	! unit number for file containing turbulence generator data changed to 90
105	! bugfix: preprocessor directives added for MPI specific code
106	!
107	! 2576 2017-10-24 13:49:46Z Giersch
108	! Definition of a new function called stg_skip_global to skip module
109	! parameters during reading restart data
110	!
111	! 2563 2017-10-19 15:36:10Z Giersch
112	! stg_read_restart_data is called in stg_parin in the case of a restart run
113	!
114	! 2259 2017-06-08 09:09:11Z gronemeier
115	! Initial revision
116	!
117	!
118	!
119	! Authors:
120	! --------
121	! @author Tobias Gronemeier, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
122	!
123	!
124	! Description:
125	! ------------
126	!> The module generates turbulence at the inflow boundary based on a method by
127	!> Xie and Castro (2008) utilizing a Lund rotation (Lund, 1998) and a mass-flux
128	!> correction by Kim et al. (2013).
129	!> The turbulence is correlated based on length scales in y- and z-direction and
130	!> a time scale for each velocity component. The profiles of length and time
131	!> scales, mean u, v, w, e and pt, and all components of the Reynolds stress
132	!> tensor are read from file STG_PROFILES.
133	!>
134	!> @todo test restart
135	!> enable cyclic_fill
136	!> implement turbulence generation for e and pt
137	!> @todo Input of height-constant length scales via namelist
138	!> @note <Enter notes on the module>
139	!> @bug Height information from input file is not used. Profiles from input
140	!> must match with current PALM grid.
141	!> In case of restart, velocity seeds differ from precursor run if a11,
142	!> a22, or a33 are zero.
143	!------------------------------------------------------------------------------!
144	MODULE synthetic_turbulence_generator_mod
145
146
147	USE arrays_3d, &
148	ONLY: mean_inflow_profiles, u, v, w
149
150	USE basic_constants_and_equations_mod, &
151	ONLY: pi
152
153	USE control_parameters, &
154	ONLY: initializing_actions, message_string, syn_turb_gen
155
156	USE cpulog, &
157	ONLY: cpu_log, log_point, log_point_s
158
159	USE indices, &
160	ONLY: nbgp, nzb, nzt, nxl, nxlg, nxr, nxrg, nys, nyn, nyng, nysg
161
162	USE kinds
163
164	#if defined( __parallel ) && !defined( __mpifh )
165	USE MPI
166	#endif
167
168	USE pegrid, &
169	ONLY: comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
170
171	USE transpose_indices, &
172	ONLY: nzb_x, nzt_x
173
174
175	IMPLICIT NONE
176
177	#if defined( __parallel ) && defined( __mpifh )
178	INCLUDE "mpif.h"
179	#endif
180
181
182	LOGICAL :: velocity_seed_initialized = .FALSE. !< true after first call of stg_main
183	LOGICAL :: use_syn_turb_gen = .FALSE. !< switch to use synthetic turbulence generator
184
185	INTEGER(iwp) :: id_stg_left !< left lateral boundary core id in case of turbulence generator
186	INTEGER(iwp) :: id_stg_north !< north lateral boundary core id in case of turbulence generator
187	INTEGER(iwp) :: id_stg_right !< right lateral boundary core id in case of turbulence generator
188	INTEGER(iwp) :: id_stg_south !< south lateral boundary core id in case of turbulence generator
189	INTEGER(iwp) :: stg_type_xz !< MPI type for full z range
190	INTEGER(iwp) :: stg_type_xz_small !< MPI type for small z range
191	INTEGER(iwp) :: stg_type_yz !< MPI type for full z range
192	INTEGER(iwp) :: stg_type_yz_small !< MPI type for small z range
193	INTEGER(iwp) :: merg !< maximum length scale (in gp)
194	INTEGER(iwp) :: mergp !< merg + nbgp
195	INTEGER(iwp) :: nzb_x_stg !< lower bound of z coordinate (required for transposing z on PEs along x)
196	INTEGER(iwp) :: nzt_x_stg !< upper bound of z coordinate (required for transposing z on PEs along x)
197	INTEGER(iwp) :: nzb_y_stg !< lower bound of z coordinate (required for transposing z on PEs along y)
198	INTEGER(iwp) :: nzt_y_stg !< upper bound of z coordinate (required for transposing z on PEs along y)
199
200	REAL(wp) :: mc_factor !< mass flux correction factor
201
202	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_xz !< displacement for MPI_GATHERV
203	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_xz !< receive count for MPI_GATHERV
204	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_yz !< displacement for MPI_GATHERV
205	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_yz !< receive count for MPI_GATHERV
206	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nux !< length scale of u in x direction (in gp)
207	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy !< length scale of u in y direction (in gp)
208	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz !< length scale of u in z direction (in gp)
209	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvx !< length scale of v in x direction (in gp)
210	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy !< length scale of v in y direction (in gp)
211	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz !< length scale of v in z direction (in gp)
212	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwx !< length scale of w in x direction (in gp)
213	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy !< length scale of w in y direction (in gp)
214	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz !< length scale of w in z direction (in gp)
215
216	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed !< seed of random number for rn-generator
217
218	REAL(wp), DIMENSION(:), ALLOCATABLE :: a11 !< coefficient for Lund rotation
219	REAL(wp), DIMENSION(:), ALLOCATABLE :: a21 !< coefficient for Lund rotation
220	REAL(wp), DIMENSION(:), ALLOCATABLE :: a22 !< coefficient for Lund rotation
221	REAL(wp), DIMENSION(:), ALLOCATABLE :: a31 !< coefficient for Lund rotation
222	REAL(wp), DIMENSION(:), ALLOCATABLE :: a32 !< coefficient for Lund rotation
223	REAL(wp), DIMENSION(:), ALLOCATABLE :: a33 !< coefficient for Lund rotation
224	REAL(wp), DIMENSION(:), ALLOCATABLE :: tu !< Lagrangian time scale of u
225	REAL(wp), DIMENSION(:), ALLOCATABLE :: tv !< Lagrangian time scale of v
226	REAL(wp), DIMENSION(:), ALLOCATABLE :: tw !< Lagrangian time scale of w
227
228	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bux !< filter function for u in x direction
229	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy !< filter function for u in y direction
230	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz !< filter function for u in z direction
231	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvx !< filter function for v in x direction
232	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy !< filter function for v in y direction
233	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz !< filter function for v in z direction
234	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwx !< filter function for w in y direction
235	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy !< filter function for w in y direction
236	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz !< filter function for w in z direction
237	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_xz !< velocity seed for u at xz plane
238	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_xz !< velocity seed for u at xz plane with new random number
239	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_yz !< velocity seed for u at yz plane
240	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_yz !< velocity seed for u at yz plane with new random number
241	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_xz !< velocity seed for v at xz plane
242	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_xz !< velocity seed for v at xz plane with new random number
243	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_yz !< velocity seed for v at yz plane
244	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_yz !< velocity seed for v at yz plane with new random number
245	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_xz !< velocity seed for w at xz plane
246	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_xz !< velocity seed for w at xz plane with new random number
247	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_yz !< velocity seed for w at yz plane
248	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz !< velocity seed for w at yz plane with new random number
249
250
251	!
252	!-- PALM interfaces:
253	!-- Input parameter checks to be done in check_parameters
254	INTERFACE stg_check_parameters
255	MODULE PROCEDURE stg_check_parameters
256	END INTERFACE stg_check_parameters
257
258	!
259	!-- Calculate filter functions
260	INTERFACE stg_filter_func
261	MODULE PROCEDURE stg_filter_func
262	END INTERFACE stg_filter_func
263
264	!
265	!-- Generate velocity seeds at south and north domain boundary
266	INTERFACE stg_generate_seed_xz
267	MODULE PROCEDURE stg_generate_seed_xz
268	END INTERFACE stg_generate_seed_xz
269	!
270	!-- Generate velocity seeds at left and/or right domain boundary
271	INTERFACE stg_generate_seed_yz
272	MODULE PROCEDURE stg_generate_seed_yz
273	END INTERFACE stg_generate_seed_yz
274
275	!
276	!-- Output of information to the header file
277	INTERFACE stg_header
278	MODULE PROCEDURE stg_header
279	END INTERFACE stg_header
280
281	!
282	!-- Initialization actions
283	INTERFACE stg_init
284	MODULE PROCEDURE stg_init
285	END INTERFACE stg_init
286
287	!
288	!-- Main procedure of synth. turb. gen.
289	INTERFACE stg_main
290	MODULE PROCEDURE stg_main
291	END INTERFACE stg_main
292
293	!
294	!-- Reading of NAMELIST parameters
295	INTERFACE stg_parin
296	MODULE PROCEDURE stg_parin
297	END INTERFACE stg_parin
298
299	!
300	!-- Reading of parameters for restart runs
301	INTERFACE stg_rrd_global
302	MODULE PROCEDURE stg_rrd_global
303	END INTERFACE stg_rrd_global
304
305	!
306	!-- Writing of binary output for restart runs
307	INTERFACE stg_wrd_global
308	MODULE PROCEDURE stg_wrd_global
309	END INTERFACE stg_wrd_global
310
311	SAVE
312
313	PRIVATE
314
315	!
316	!-- Public interfaces
317	PUBLIC stg_check_parameters, stg_header, stg_init, stg_main, stg_parin, &
318	stg_wrd_global, stg_rrd_global
319
320	!
321	!-- Public variables
322	PUBLIC id_stg_left, id_stg_north, id_stg_right, id_stg_south, &
323	use_syn_turb_gen
324
325
326	CONTAINS
327
328
329	!------------------------------------------------------------------------------!
330	! Description:
331	! ------------
332	!> Check parameters routine for synthetic turbulence generator
333	!------------------------------------------------------------------------------!
334	SUBROUTINE stg_check_parameters
335
336
337	USE control_parameters, &
338	ONLY: bc_lr, bc_ns, child_domain, nesting_offline, rans_mode, &
339	turbulent_inflow
340
341	USE pmc_interface, &
342	ONLY : rans_mode_parent
343
344
345	IMPLICIT NONE
346
347	IF ( .NOT. use_syn_turb_gen .AND. .NOT. rans_mode .AND. &
348	nesting_offline ) THEN
349	message_string = 'No synthetic turbulence generator is applied. ' // &
350	'In case PALM operates in LES mode and lateral ' // &
351	'boundary conditions are provided by COSMO model, ' // &
352	'turbulence may require large adjustment lenght at ' //&
353	'the lateral inflow boundaries. Please check your ' // &
354	'results carefully.'
355	CALL message( 'stg_check_parameters', 'PA0000', 0, 0, 0, 6, 0 )
356	ENDIF
357
358	IF ( .NOT. use_syn_turb_gen .AND. child_domain &
359	.AND. rans_mode_parent .AND. .NOT. rans_mode ) THEN
360	message_string = 'Synthetic turbulence generator has to be applied ' // &
361	'when nesting is applied and parent operates in ' // &
362	'RANS-mode but current child in LES mode.'
363	CALL message( 'stg_check_parameters', 'PA0000', 1, 2, 0, 6, 0 )
364	ENDIF
365
366	IF ( use_syn_turb_gen ) THEN
367
368	IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN
369
370	IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0 &
371	.AND. INDEX( initializing_actions, 'read_restart_data' ) == 0 ) THEN
372	message_string = 'Using synthetic turbulence generator ' // &
373	'requires %initializing_actions = ' // &
374	'"set_constant_profiles" or "read_restart_data"'
375	CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
376	ENDIF
377
378	IF ( bc_lr /= 'dirichlet/radiation' ) THEN
379	message_string = 'Using synthetic turbulence generator ' // &
380	'requires &bc_lr = "dirichlet/radiation"'
381	CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
382	ENDIF
383	IF ( bc_ns /= 'cyclic' ) THEN
384	message_string = 'Using synthetic turbulence generator ' // &
385	'requires &bc_ns = "cyclic"'
386	CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
387	ENDIF
388
389	ENDIF
390
391	IF ( turbulent_inflow ) THEN
392	message_string = 'Using synthetic turbulence generator ' // &
393	'in combination &with turbulent_inflow = .T. '// &
394	'is not allowed'
395	CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
396	ENDIF
397
398	ENDIF
399
400	END SUBROUTINE stg_check_parameters
401
402
403	!------------------------------------------------------------------------------!
404	! Description:
405	! ------------
406	!> Header output for synthetic turbulence generator
407	!------------------------------------------------------------------------------!
408	SUBROUTINE stg_header ( io )
409
410
411	IMPLICIT NONE
412
413	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
414
415	!
416	!-- Write synthetic turbulence generator Header
417	WRITE( io, 1 )
418	IF ( use_syn_turb_gen ) THEN
419	WRITE( io, 2 )
420	ELSE
421	WRITE( io, 3 )
422	ENDIF
423
424	1 FORMAT (//' Synthetic turbulence generator information:'/ &
425	' ------------------------------------------'/)
426	2 FORMAT (' synthetic turbulence generator switched on')
427	3 FORMAT (' synthetic turbulence generator switched off')
428
429	END SUBROUTINE stg_header
430
431
432	!------------------------------------------------------------------------------!
433	! Description:
434	! ------------
435	!> Initialization of the synthetic turbulence generator
436	!------------------------------------------------------------------------------!
437	SUBROUTINE stg_init
438
439
440	USE arrays_3d, &
441	ONLY: ddzw, u_init, v_init, zu
442
443	USE control_parameters, &
444	ONLY: child_domain, coupling_char, e_init, nesting_offline, rans_mode
445
446	USE grid_variables, &
447	ONLY: ddx, ddy
448
449	USE indices, &
450	ONLY: nz
451
452	USE pmc_interface, &
453	ONLY : rans_mode_parent
454
455
456	IMPLICIT NONE
457
458	LOGICAL :: file_stg_exist = .FALSE. !< flag indication whether parameter file for Reynolds stress and length scales exist
459
460	#if defined( __parallel )
461	INTEGER(KIND=MPI_ADDRESS_KIND) :: extent !< extent of new MPI type
462	INTEGER(KIND=MPI_ADDRESS_KIND) :: tob=0 !< dummy variable
463	#endif
464
465	INTEGER(iwp) :: i !> grid index in x-direction
466	INTEGER(iwp) :: j !> loop index
467	INTEGER(iwp) :: k !< index
468	INTEGER(iwp) :: newtype !< dummy MPI type
469	INTEGER(iwp) :: realsize !< size of REAL variables
470	INTEGER(iwp) :: nseed !< dimension of random number seed
471	INTEGER(iwp) :: startseed = 1234567890 !< start seed for random number generator
472
473	!
474	!-- Dummy variables used for reading profiles
475	REAL(wp) :: d1 !< u profile
476	REAL(wp) :: d2 !< v profile
477	REAL(wp) :: d3 !< w profile
478	REAL(wp) :: d5 !< e profile
479	REAL(wp) :: d11 !< vertical interpolation for a11
480	REAL(wp) :: d21 !< vertical interpolation for a21
481	REAL(wp) :: d22 !< vertical interpolation for a22
482	REAL(wp) :: dum_exp !< dummy variable used for exponential vertical decrease of turbulent length scales
483	REAL(wp) :: luy !< length scale for u in y direction
484	REAL(wp) :: luz !< length scale for u in z direction
485	REAL(wp) :: lvy !< length scale for v in y direction
486	REAL(wp) :: lvz !< length scale for v in z direction
487	REAL(wp) :: lwy !< length scale for w in y direction
488	REAL(wp) :: lwz !< length scale for w in z direction
489	REAL(wp) :: nnz !< increment used to determine processor decomposition of z-axis along x and y direction
490	REAL(wp) :: zz !< height
491
492	REAL(wp) :: length_scale_surface !< typical length scale
493	REAL(wp) :: r_ii_0 !< correction factor
494	REAL(wp) :: time_scale !< typical time scale
495	REAL(wp) :: length_scale_z !< typical length scale
496
497	REAL(wp),DIMENSION(nzb:nzt+1) :: r11 !< Reynolds parameter
498	REAL(wp),DIMENSION(nzb:nzt+1) :: r21 !< Reynolds parameter
499	REAL(wp),DIMENSION(nzb:nzt+1) :: r22 !< Reynolds parameter
500	REAL(wp),DIMENSION(nzb:nzt+1) :: r31 !< Reynolds parameter
501	REAL(wp),DIMENSION(nzb:nzt+1) :: r32 !< Reynolds parameter
502	REAL(wp),DIMENSION(nzb:nzt+1) :: r33 !< Reynolds parameter
503
504
505	#if defined( __parallel )
506	CALL MPI_BARRIER( comm2d, ierr )
507	#endif
508
509	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
510
511	IF ( .NOT. ALLOCATED( mean_inflow_profiles ) ) &
512	ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )
513
514	ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1), &
515	a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1), &
516	nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1), &
517	nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1), &
518	nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1), &
519	tu(nzb:nzt+1), tv(nzb:nzt+1), tw(nzb:nzt+1) )
520
521	#if defined( __parallel )
522	!
523	!-- Determine processor decomposition of z-axis along x- and y-direction
524	nnz = nz / pdims(1)
525	nzb_x_stg = 1 + myidx * INT( nnz )
526	nzt_x_stg = ( myidx + 1 ) * INT( nnz )
527
528	IF ( MOD( nz , pdims(1) ) /= 0 .AND. myidx == id_stg_right ) &
529	nzt_x_stg = nzt_x_stg + myidx * ( nnz - INT( nnz ) )
530	! nzt_x_stg = myidx * nnz + MOD( nz , pdims(1) )
531
532	IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent &
533	.AND. .NOT. rans_mode ) ) THEN
534	nnz = nz / pdims(2)
535	nzb_y_stg = 1 + myidy * INT( nnz )
536	nzt_y_stg = ( myidy + 1 ) * INT( nnz )
537
538	IF ( MOD( nz , pdims(2) ) /= 0 .AND. myidy == id_stg_north ) &
539	nzt_y_stg = nzt_y_stg + myidy * ( nnz - INT( nnz ) )
540	! nzt_y_stg = myidy * nnz + MOD( nz , pdims(2) )
541	ENDIF
542
543	!
544	!-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted
545	!-- velocity seeds
546	CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
547	extent = 1 * realsize
548	!
549	!-- Set-up MPI datatyp to involve all cores for turbulence generation at yz
550	!-- layer
551	!-- stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
552	CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyng-nysg+1], &
553	[1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
554	CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz, ierr )
555	CALL MPI_TYPE_COMMIT( stg_type_yz, ierr )
556	CALL MPI_TYPE_FREE( newtype, ierr )
557
558	! stg_type_yz_small: yz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
559	CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyng-nysg+1], &
560	[1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
561	CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr )
562	CALL MPI_TYPE_COMMIT( stg_type_yz_small, ierr )
563	CALL MPI_TYPE_FREE( newtype, ierr )
564
565	! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
566	ALLOCATE( recv_count_yz(pdims(1)), displs_yz(pdims(1)) )
567
568	recv_count_yz = nzt_x_stg-nzb_x_stg + 1
569	recv_count_yz(pdims(1)) = recv_count_yz(pdims(1)) + 1
570
571	DO j = 1, pdims(1)
572	displs_yz(j) = 0 + (nzt_x_stg-nzb_x_stg+1) * (j-1)
573	ENDDO
574	!
575	!-- Set-up MPI datatyp to involve all cores for turbulence generation at xz
576	!-- layer
577	!-- stg_type_xz: xz-slice with vertical bounds nzb:nzt+1
578	IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent &
579	.AND. .NOT. rans_mode ) ) THEN
580	CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nxrg-nxlg+1], &
581	[1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
582	CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz, ierr )
583	CALL MPI_TYPE_COMMIT( stg_type_xz, ierr )
584	CALL MPI_TYPE_FREE( newtype, ierr )
585
586	! stg_type_yz_small: xz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
587	CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_y_stg-nzb_y_stg+2,nxrg-nxlg+1], &
588	[1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
589	CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz_small, ierr )
590	CALL MPI_TYPE_COMMIT( stg_type_xz_small, ierr )
591	CALL MPI_TYPE_FREE( newtype, ierr )
592
593	! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
594	ALLOCATE( recv_count_xz(pdims(2)), displs_xz(pdims(2)) )
595
596	recv_count_xz = nzt_y_stg-nzb_y_stg + 1
597	recv_count_xz(pdims(2)) = recv_count_xz(pdims(2)) + 1
598
599	DO j = 1, pdims(2)
600	displs_xz(j) = 0 + (nzt_y_stg-nzb_y_stg+1) * (j-1)
601	ENDDO
602
603	ENDIF
604
605	#endif
606	!
607	!-- Define seed of random number
608	CALL RANDOM_SEED()
609	CALL RANDOM_SEED( size=nseed )
610	ALLOCATE( seed(1:nseed) )
611	DO j = 1, nseed
612	seed(j) = startseed + j
613	ENDDO
614	CALL RANDOM_SEED(put=seed)
615
616	!-- Read inflow profile
617	!-- Height levels of profiles in input profile is as follows:
618	!-- zu: luy, luz, tu, lvy, lvz, tv, r11, r21, r22, d1, d2, d5
619	!-- zw: lwy, lwz, tw, r31, r32, r33, d3
620	!-- WARNING: zz is not used at the moment
621	INQUIRE( FILE = 'STG_PROFILES' // TRIM( coupling_char ), &
622	EXIST = file_stg_exist )
623
624	IF ( file_stg_exist ) THEN
625
626	OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD', &
627	FORM='FORMATTED')
628	!
629	!-- Skip header
630	READ( 90, * )
631
632	DO k = nzb+1, nzt+1
633	READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k), &
634	r11(k), r21(k), r22(k), r31(k), r32(k), r33(k), &
635	d1, d2, d3, d5
636
637	!
638	!-- Convert length scales from meter to number of grid points.
639	nuy(k) = INT( luy * ddy )
640	nuz(k) = INT( luz * ddzw(k) )
641	nvy(k) = INT( lvy * ddy )
642	nvz(k) = INT( lvz * ddzw(k) )
643	nwy(k) = INT( lwy * ddy )
644	nwz(k) = INT( lwz * ddzw(k) )
645	!
646	!-- Workaround, assume isotropic turbulence
647	nwx(k) = nwy(k)
648	nvx(k) = nvy(k)
649	nux(k) = nuy(k)
650	!
651	!-- Save Mean inflow profiles
652	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
653	mean_inflow_profiles(k,1) = d1
654	mean_inflow_profiles(k,2) = d2
655	! mean_inflow_profiles(k,4) = d4
656	mean_inflow_profiles(k,5) = d5
657	ENDIF
658	ENDDO
659	!
660	!-- Set lenght scales at surface grid point
661	nuy(nzb) = nuy(nzb+1)
662	nuz(nzb) = nuz(nzb+1)
663	nvy(nzb) = nvy(nzb+1)
664	nvz(nzb) = nvz(nzb+1)
665	nwy(nzb) = nwy(nzb+1)
666	nwz(nzb) = nwz(nzb+1)
667
668	CLOSE( 90 )
669
670	ELSE
671	!
672	!-- Set-up defaul length scales. Assume exponentially decreasing length
673	!-- scales and isotropic turbulence.
674	!-- Typical length (time) scales of 100 m (s) should be a good compromise
675	!-- between all stratrifications. Near-surface variances are fixed to
676	!-- 0.1 m2/s2, vertical fluxes are one order of magnitude smaller.
677	!-- Vertical fluxes
678	length_scale_surface = 100.0_wp
679	r_ii_0 = 0.1_wp
680	time_scale = 100.0_wp
681	DO k = nzb+1, nzt+1
682	dum_exp = MERGE( -zu(k) / ( 0.3* zu(nzt) ), &
683	0.0_wp, &
684	zu(k) > 0.3 * zu(nzt) &
685	)
686	length_scale_z = length_scale_surface * EXP( dum_exp )
687
688	nux(k) = MAX( INT( length_scale_z * ddx ), 1 )
689	nuy(k) = MAX( INT( length_scale_z * ddy ), 1 )
690	nuz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
691	nvx(k) = MAX( INT( length_scale_z * ddx ), 1 )
692	nvy(k) = MAX( INT( length_scale_z * ddy ), 1 )
693	nvz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
694	nwx(k) = MAX( INT( length_scale_z * ddx ), 1 )
695	nwy(k) = MAX( INT( length_scale_z * ddy ), 1 )
696	nwz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
697
698	r11(k) = r_ii_0 * EXP( dum_exp )
699	r22(k) = r_ii_0 * EXP( dum_exp )
700	r33(k) = r_ii_0 * EXP( dum_exp )
701
702	r21(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
703	r31(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
704	r32(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
705
706	tu(k) = time_scale
707	tv(k) = time_scale
708	tw(k) = time_scale
709
710	ENDDO
711	nux(nzb) = nux(nzb+1)
712	nuy(nzb) = nuy(nzb+1)
713	nuz(nzb) = nuz(nzb+1)
714	nvx(nzb) = nvx(nzb+1)
715	nvy(nzb) = nvy(nzb+1)
716	nvz(nzb) = nvz(nzb+1)
717	nwx(nzb) = nwx(nzb+1)
718	nwy(nzb) = nwy(nzb+1)
719	nwz(nzb) = nwz(nzb+1)
720
721	r11(nzb) = r11(nzb+1)
722	r22(nzb) = r22(nzb+1)
723	r33(nzb) = r33(nzb+1)
724
725	r21(nzb) = r11(nzb+1)
726	r31(nzb) = r31(nzb+1)
727	r32(nzb) = r32(nzb+1)
728
729	tu(nzb) = time_scale
730	tv(nzb) = time_scale
731	tw(nzb) = time_scale
732
733	ENDIF
734
735	!
736	!-- Assign initial profiles
737	IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN
738	u_init = mean_inflow_profiles(:,1)
739	v_init = mean_inflow_profiles(:,2)
740	!pt_init = mean_inflow_profiles(:,4)
741	e_init = MAXVAL( mean_inflow_profiles(:,5) )
742	ENDIF
743	!
744	!-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
745	DO k = nzb, nzt+1
746	IF ( r11(k) > 0.0_wp ) THEN
747	a11(k) = SQRT( r11(k) )
748	a21(k) = r21(k) / a11(k)
749	ELSE
750	a11(k) = 0.0_wp
751	a21(k) = 0.0_wp
752	ENDIF
753
754	a22(k) = r22(k) - a21(k)**2
755	IF ( a22(k) > 0.0_wp ) THEN
756	a22(k) = SQRT( a22(k) )
757	ELSE
758	a22(k) = 0.0_wp
759	ENDIF
760
761	!
762	!-- a31, a32, a33 must be calculated with interpolated a11, a21, a22 (d11,
763	!-- d21, d22) because of different vertical grid
764	IF ( k .le. nzt ) THEN
765	d11 = 0.5_wp * ( r11(k) + r11(k+1) )
766	IF ( d11 > 0.0_wp ) THEN
767	d11 = SQRT( d11 )
768	d21 = ( 0.5_wp * ( r21(k) + r21(k+1) ) ) / d11
769	a31(k) = r31(k) / d11
770	ELSE
771	d21 = 0.0_wp
772	a31(k) = 0.0_wp
773	ENDIF
774
775	d22 = 0.5_wp * ( r22(k) + r22(k+1) ) - d21 ** 2
776	IF ( d22 > 0.0_wp ) THEN
777	a32(k) = ( r32(k) - d21 * a31(k) ) / SQRT( d22 )
778	ELSE
779	a32(k) = 0.0_wp
780	ENDIF
781
782	a33(k) = r33(k) - a31(k) 2 - a32(k) 2
783	IF ( a33(k) > 0.0_wp ) THEN
784	a33(k) = SQRT( a33(k) )
785	ELSE
786	a33(k) = 0.0_wp
787	ENDIF
788	ELSE
789	a31(k) = a31(k-1)
790	a32(k) = a32(k-1)
791	a33(k) = a33(k-1)
792	ENDIF
793
794	ENDDO
795	!
796	!-- Define the size of the filter functions and allocate them.
797	merg = 0
798
799	! arrays must be large enough to cover the largest length scale
800	DO k = nzb, nzt+1
801	j = MAX( ABS(nux(k)), ABS(nuy(k)), ABS(nuz(k)), &
802	ABS(nvx(k)), ABS(nvy(k)), ABS(nvz(k)), &
803	ABS(nwx(k)), ABS(nwy(k)), ABS(nwz(k)) )
804	IF ( j > merg ) merg = j
805	ENDDO
806
807	merg = 2 * merg
808	mergp = merg + nbgp
809
810	ALLOCATE ( bux(-merg:merg,nzb:nzt+1), &
811	buy(-merg:merg,nzb:nzt+1), &
812	buz(-merg:merg,nzb:nzt+1), &
813	bvx(-merg:merg,nzb:nzt+1), &
814	bvy(-merg:merg,nzb:nzt+1), &
815	bvz(-merg:merg,nzb:nzt+1), &
816	bwx(-merg:merg,nzb:nzt+1), &
817	bwy(-merg:merg,nzb:nzt+1), &
818	bwz(-merg:merg,nzb:nzt+1) )
819
820	!
821	!-- Allocate velocity seeds for turbulence at xz-layer
822	ALLOCATE ( fu_xz( nzb:nzt+1,nxlg:nxrg), fuo_xz(nzb:nzt+1,nxlg:nxrg), &
823	fv_xz( nzb:nzt+1,nxlg:nxrg), fvo_xz(nzb:nzt+1,nxlg:nxrg), &
824	fw_xz( nzb:nzt+1,nxlg:nxrg), fwo_xz(nzb:nzt+1,nxlg:nxrg) )
825
826	!
827	!-- Allocate velocity seeds for turbulence at yz-layer
828	ALLOCATE ( fu_yz( nzb:nzt+1,nysg:nyng), fuo_yz(nzb:nzt+1,nysg:nyng), &
829	fv_yz( nzb:nzt+1,nysg:nyng), fvo_yz(nzb:nzt+1,nysg:nyng), &
830	fw_yz( nzb:nzt+1,nysg:nyng), fwo_yz(nzb:nzt+1,nysg:nyng) )
831
832	fu_xz = 0.0_wp
833	fuo_xz = 0.0_wp
834	fv_xz = 0.0_wp
835	fvo_xz = 0.0_wp
836	fw_xz = 0.0_wp
837	fwo_xz = 0.0_wp
838
839	fu_yz = 0.0_wp
840	fuo_yz = 0.0_wp
841	fv_yz = 0.0_wp
842	fvo_yz = 0.0_wp
843	fw_yz = 0.0_wp
844	fwo_yz = 0.0_wp
845
846	!
847	!-- Create filter functions
848	CALL stg_filter_func( nux, bux ) !filter ux
849	CALL stg_filter_func( nuy, buy ) !filter uy
850	CALL stg_filter_func( nuz, buz ) !filter uz
851	CALL stg_filter_func( nvx, bvx ) !filter vx
852	CALL stg_filter_func( nvy, bvy ) !filter vy
853	CALL stg_filter_func( nvz, bvz ) !filter vz
854	CALL stg_filter_func( nwx, bwx ) !filter wx
855	CALL stg_filter_func( nwy, bwy ) !filter wy
856	CALL stg_filter_func( nwz, bwz ) !filter wz
857
858	#if defined( __parallel )
859	CALL MPI_BARRIER( comm2d, ierr )
860	#endif
861
862	!
863	!-- In case of restart, calculate velocity seeds fu, fv, fw from former
864	! time step.
865	! Bug: fu, fv, fw are different in those heights where a11, a22, a33
866	! are 0 compared to the prerun. This is mostly for k=nzt+1.
867	IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN
868	IF ( myidx == id_stg_left .OR. myidx == id_stg_right ) THEN
869
870	IF ( myidx == id_stg_left ) i = -1
871	IF ( myidx == id_stg_right ) i = nxr+1
872
873	DO j = nysg, nyng
874	DO k = nzb, nzt+1
875
876	IF ( a11(k) .NE. 0._wp ) THEN
877	fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
878	ELSE
879	fu_yz(k,j) = 0._wp
880	ENDIF
881
882	IF ( a22(k) .NE. 0._wp ) THEN
883	fv_yz(k,j) = ( v(k,j,i) / mc_factor - a21(k) * fu_yz(k,j) - &
884	v_init(k) ) / a22(k)
885	ELSE
886	fv_yz(k,j) = 0._wp
887	ENDIF
888
889	IF ( a33(k) .NE. 0._wp ) THEN
890	fw_yz(k,j) = ( w(k,j,i) / mc_factor - a31(k) * fu_yz(k,j) - &
891	a32(k) * fv_yz(k,j) ) / a33(k)
892	ELSE
893	fw_yz = 0._wp
894	ENDIF
895
896	ENDDO
897	ENDDO
898	ENDIF
899	ENDIF
900
901	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
902
903	END SUBROUTINE stg_init
904
905
906	!------------------------------------------------------------------------------!
907	! Description:
908	! ------------
909	!> Calculate filter function bxx from length scale nxx following Eg.9 and 10
910	!> (Xie and Castro, 2008)
911	!------------------------------------------------------------------------------!
912	SUBROUTINE stg_filter_func( nxx, bxx )
913
914
915	IMPLICIT NONE
916
917	INTEGER(iwp) :: k !< loop index
918	INTEGER(iwp) :: n_k !< length scale nXX in height k
919	INTEGER(iwp) :: n_k2 !< n_k * 2
920	INTEGER(iwp) :: nf !< index for length scales
921
922	REAL(wp) :: bdenom !< denominator for filter functions bXX
923	REAL(wp) :: qsi = 1.0_wp !< minimization factor
924
925	INTEGER(iwp), DIMENSION(:) :: nxx(nzb:nzt+1) !< length scale (in gp)
926
927	REAL(wp), DIMENSION(:,:) :: bxx(-merg:merg,nzb:nzt+1) !< filter function
928
929
930	bxx = 0.0_wp
931
932	DO k = nzb, nzt+1
933	bdenom = 0.0_wp
934	n_k = nxx(k)
935	IF ( n_k /= 0 ) THEN
936	n_k2 = n_k * 2
937
938	!
939	!-- ( Eq.10 )^2
940	DO nf = -n_k2, n_k2
941	bdenom = bdenom + EXP( -qsi * pi * ABS(nf) / n_k )**2
942	ENDDO
943
944	!
945	!-- ( Eq.9 )
946	bdenom = SQRT( bdenom )
947	DO nf = -n_k2, n_k2
948	bxx(nf,k) = EXP( -qsi * pi * ABS(nf) / n_k ) / bdenom
949	ENDDO
950	ENDIF
951	ENDDO
952
953	END SUBROUTINE stg_filter_func
954
955
956	!------------------------------------------------------------------------------!
957	! Description:
958	! ------------
959	!> Parin for &stg_par for synthetic turbulence generator
960	!------------------------------------------------------------------------------!
961	SUBROUTINE stg_parin
962
963
964	IMPLICIT NONE
965
966	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
967
968
969	NAMELIST /stg_par/ use_syn_turb_gen
970
971	line = ' '
972
973	!
974	!-- Try to find stg package
975	REWIND ( 11 )
976	line = ' '
977	DO WHILE ( INDEX( line, '&stg_par' ) == 0 )
978	READ ( 11, '(A)', END=20 ) line
979	ENDDO
980	BACKSPACE ( 11 )
981
982	!
983	!-- Read namelist
984	READ ( 11, stg_par, ERR = 10, END = 20 )
985
986	!
987	!-- Set flag that indicates that the synthetic turbulence generator is switched
988	!-- on
989	syn_turb_gen = .TRUE.
990	GOTO 20
991
992	10 BACKSPACE( 11 )
993	READ( 11 , '(A)') line
994	CALL parin_fail_message( 'stg_par', line )
995
996	20 CONTINUE
997
998	END SUBROUTINE stg_parin
999
1000
1001	!------------------------------------------------------------------------------!
1002	! Description:
1003	! ------------
1004	!> This routine reads the respective restart data.
1005	!------------------------------------------------------------------------------!
1006	SUBROUTINE stg_rrd_global( found )
1007
1008
1009	USE control_parameters, &
1010	ONLY: length, restart_string
1011
1012
1013	IMPLICIT NONE
1014
1015	LOGICAL, INTENT(OUT) :: found !< flag indicating if variable was found
1016
1017
1018	found = .TRUE.
1019
1020
1021	SELECT CASE ( restart_string(1:length) )
1022
1023	CASE ( 'mc_factor' )
1024	READ ( 13 ) mc_factor
1025	CASE ( 'use_syn_turb_gen' )
1026	READ ( 13 ) use_syn_turb_gen
1027
1028	CASE DEFAULT
1029
1030	found = .FALSE.
1031
1032	END SELECT
1033
1034
1035	END SUBROUTINE stg_rrd_global
1036
1037
1038	!------------------------------------------------------------------------------!
1039	! Description:
1040	! ------------
1041	!> This routine writes the respective restart data.
1042	!------------------------------------------------------------------------------!
1043	SUBROUTINE stg_wrd_global
1044
1045
1046	IMPLICIT NONE
1047
1048	CALL wrd_write_string( 'mc_factor' )
1049	WRITE ( 14 ) mc_factor
1050
1051	CALL wrd_write_string( 'use_syn_turb_gen' )
1052	WRITE ( 14 ) use_syn_turb_gen
1053
1054
1055	END SUBROUTINE stg_wrd_global
1056
1057
1058	!------------------------------------------------------------------------------!
1059	! Description:
1060	! ------------
1061	!> Create turbulent inflow fields for u, v, w with prescribed length scales and
1062	!> Reynolds stress tensor after a method of Xie and Castro (2008), modified
1063	!> following suggestions of Kim et al. (2013), and using a Lund rotation
1064	!> (Lund, 1998).
1065	!------------------------------------------------------------------------------!
1066	SUBROUTINE stg_main
1067
1068
1069	USE arrays_3d, &
1070	ONLY: dzw
1071
1072	USE control_parameters, &
1073	ONLY: child_domain, dt_3d, intermediate_timestep_count, &
1074	nesting_offline, rans_mode, simulated_time, volume_flow_initial
1075
1076	USE grid_variables, &
1077	ONLY: dx, dy
1078
1079	USE indices, &
1080	ONLY: wall_flags_0
1081
1082	USE statistics, &
1083	ONLY: weight_substep
1084
1085	USE pmc_interface, &
1086	ONLY : rans_mode_parent
1087
1088
1089	IMPLICIT NONE
1090
1091	INTEGER(iwp) :: i !< grid index in x-direction
1092	INTEGER(iwp) :: j !< loop index in y-direction
1093	INTEGER(iwp) :: k !< loop index in z-direction
1094
1095	REAL(wp) :: dt_stg !< wheighted subtimestep
1096	REAL(wp) :: mc_factor_l !< local mass flux correction factor
1097	REAL(wp) :: volume_flow !< mass flux through lateral boundary
1098	REAL(wp) :: volume_flow_l !< local mass flux through lateral boundary
1099
1100	REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg,5) :: dist_xz !< imposed disturbances at north/south boundary
1101	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5) :: dist_yz !< imposed disturbances at left/right boundary
1102
1103
1104	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
1105
1106	!
1107	!-- Calculate time step which is needed for filter functions
1108	dt_stg = dt_3d * weight_substep(intermediate_timestep_count)
1109
1110	!
1111	!-- Initial value of fu, fv, fw
1112	IF ( simulated_time == 0.0_wp .AND. .NOT. velocity_seed_initialized ) THEN
1113	CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
1114	CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
1115	CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
1116
1117	IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent &
1118	.AND. .NOT. rans_mode ) ) THEN
1119	!
1120	!-- Generate turbulence at right boundary
1121	CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_right )
1122	CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_right )
1123	CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_right )
1124	!
1125	!-- Generate turbulence at north boundary
1126	CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_north )
1127	CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_north )
1128	CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_north )
1129	!
1130	!-- Generate turbulence at south boundary
1131	CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_south )
1132	CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_south )
1133	CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_south )
1134	ENDIF
1135	velocity_seed_initialized = .TRUE.
1136	ENDIF
1137	!
1138	!-- New set of fu, fv, fw
1139	CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
1140	CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
1141	CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
1142
1143	IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent &
1144	.AND. .NOT. rans_mode ) ) THEN
1145	!
1146	!-- Generate turbulence at right boundary
1147	CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
1148	CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
1149	CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
1150	!
1151	!-- Generate turbulence at north boundary
1152	CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
1153	CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
1154	CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
1155	!
1156	!-- Generate turbulence at south boundary
1157	CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
1158	CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
1159	CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
1160	ENDIF
1161	!
1162	!-- Turbulence generation at left and or right boundary
1163	IF ( myidx == id_stg_left .OR. myidx == id_stg_right ) THEN
1164
1165	DO j = nysg, nyng
1166	DO k = nzb, nzt + 1
1167	!
1168	!-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1169	IF ( tu(k) == 0.0_wp ) THEN
1170	fu_yz(k,j) = fuo_yz(k,j)
1171	ELSE
1172	fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
1173	fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
1174	ENDIF
1175
1176	IF ( tv(k) == 0.0_wp ) THEN
1177	fv_yz(k,j) = fvo_yz(k,j)
1178	ELSE
1179	fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
1180	fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
1181	ENDIF
1182
1183	IF ( tw(k) == 0.0_wp ) THEN
1184	fw_yz(k,j) = fwo_yz(k,j)
1185	ELSE
1186	fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
1187	fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
1188	ENDIF
1189	!
1190	!-- Lund rotation following Eq. 17 in Xie and Castro (2008).
1191	!-- Additional factors are added to improve the variance of v and w
1192	IF( k == 0 ) THEN
1193	dist_yz(k,j,1) = 0.0_wp
1194	dist_yz(k,j,2) = 0.0_wp
1195	dist_yz(k,j,3) = 0.0_wp
1196	! dist_yz(k,j,4) = 0.0_wp
1197	! dist_yz(k,j,5) = 0.0_wp
1198	ELSE
1199	dist_yz(k,j,1) = a11(k) * fu_yz(k,j)
1200	!experimental test of 1.2
1201	dist_yz(k,j,2) = ( SQRT( a22(k) / MAXVAL(a22) ) &
1202	* 1.2_wp ) &
1203	* ( a21(k) * fu_yz(k,j) &
1204	+ a22(k) * fv_yz(k,j) )
1205	dist_yz(k,j,3) = ( SQRT(a33(k) / MAXVAL(a33) ) &
1206	* 1.3_wp ) &
1207	* ( a31(k) * fu_yz(k,j) &
1208	+ a32(k) * fv_yz(k,j) &
1209	+ a33(k) * fw_yz(k,j) )
1210	! Calculation for pt and e not yet implemented
1211	! dist_yz(k,j,4) = 0.0_wp
1212	! dist_yz(k,j,5) = 0.0_wp
1213	ENDIF
1214
1215	ENDDO
1216	ENDDO
1217
1218	!
1219	!-- Mass flux correction following Kim et al. (2013)
1220	!-- This correction factor insures that the mass flux is preserved at the
1221	!-- inflow boundary
1222	IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN
1223	mc_factor_l = 0.0_wp
1224	mc_factor = 0.0_wp
1225	DO j = nys, nyn
1226	DO k = nzb+1, nzt
1227	mc_factor_l = mc_factor_l + dzw(k) * &
1228	( mean_inflow_profiles(k,1) + dist_yz(k,j,1) )
1229	ENDDO
1230	ENDDO
1231
1232	#if defined( __parallel )
1233	CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, &
1234	1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1235	#else
1236	mc_factor = mc_factor_l
1237	#endif
1238
1239	mc_factor = volume_flow_initial(1) / mc_factor
1240
1241	!
1242	!-- Add disturbance at the inflow
1243	DO j = nysg, nyng
1244	DO k = nzb, nzt+1
1245	u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) + &
1246	dist_yz(k,j,1) ) * mc_factor
1247	v(k,j,-nbgp:-1) = ( mean_inflow_profiles(k,2) + &
1248	dist_yz(k,j,2) ) * mc_factor
1249	w(k,j,-nbgp:-1) = dist_yz(k,j,3) * mc_factor
1250	ENDDO
1251	ENDDO
1252
1253	ELSE
1254	!
1255	!-- First, calculate volume flow at yz boundary
1256	IF ( myidx == id_stg_left ) i = nxl
1257	IF ( myidx == id_stg_right ) i = nxr+1
1258
1259	volume_flow_l = 0.0_wp
1260	volume_flow = 0.0_wp
1261	mc_factor_l = 0.0_wp
1262	mc_factor = 0.0_wp
1263	DO j = nys, nyn
1264	DO k = nzb+1, nzt
1265	volume_flow_l = volume_flow_l + u(k,j,i) * dzw(k) * dy &
1266	* MERGE( 1.0_wp, 0.0_wp, &
1267	BTEST( wall_flags_0(k,j,i), 1 ) )
1268
1269	mc_factor_l = mc_factor_l + ( u(k,j,i) + dist_yz(k,j,1) ) &
1270	* dzw(k) * dy &
1271	* MERGE( 1.0_wp, 0.0_wp, &
1272	BTEST( wall_flags_0(k,j,i), 1 ) )
1273	ENDDO
1274	ENDDO
1275	#if defined( __parallel )
1276	CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, &
1277	1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1278	CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, &
1279	1, MPI_REAL, MPI_SUM, comm1dy, ierr )
1280	#else
1281	volume_flow = volume_flow_l
1282	mc_factor = mc_factor_l
1283	#endif
1284
1285	mc_factor = volume_flow / mc_factor
1286
1287	!
1288	!-- Add disturbances
1289	IF ( myidx == id_stg_left ) THEN
1290
1291	DO j = nys, nyn
1292	DO k = nzb, nzt
1293	u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) ) &
1294	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1295	BTEST( wall_flags_0(k,j,0), 1 ) )
1296	v(k,j,-1) = ( v(k,j,-1) + dist_yz(k,j,2) ) &
1297	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1298	BTEST( wall_flags_0(k,j,-1), 2 ) )
1299	w(k,j,-1) = ( w(k,j,-1) + dist_yz(k,j,3) ) &
1300	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1301	BTEST( wall_flags_0(k,j,-1), 3 ) )
1302	ENDDO
1303	ENDDO
1304	ENDIF
1305	IF ( myidx == id_stg_right ) THEN
1306
1307	DO j = nys, nyn
1308	DO k = nzb, nzt
1309	u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) ) &
1310	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1311	BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
1312	v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) ) &
1313	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1314	BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
1315	w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) ) &
1316	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1317	BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
1318	ENDDO
1319	ENDDO
1320	ENDIF
1321	ENDIF
1322
1323	ENDIF
1324	!
1325	!-- Turbulence generation at north and south boundary
1326	IF ( myidy == id_stg_north .OR. myidy == id_stg_south ) THEN
1327
1328	DO i = nxlg, nxrg
1329	DO k = nzb, nzt + 1
1330	!
1331	!-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
1332	IF ( tu(k) == 0.0_wp ) THEN
1333	fu_xz(k,i) = fuo_xz(k,i)
1334	ELSE
1335	fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
1336	fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
1337	ENDIF
1338
1339	IF ( tv(k) == 0.0_wp ) THEN
1340	fv_xz(k,i) = fvo_xz(k,i)
1341	ELSE
1342	fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
1343	fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
1344	ENDIF
1345
1346	IF ( tw(k) == 0.0_wp ) THEN
1347	fw_xz(k,i) = fwo_xz(k,i)
1348	ELSE
1349	fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
1350	fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
1351	ENDIF
1352	!
1353	!-- Lund rotation following Eq. 17 in Xie and Castro (2008).
1354	!-- Additional factors are added to improve the variance of v and w
1355	IF( k == 0 ) THEN
1356	dist_xz(k,i,1) = 0.0_wp
1357	dist_xz(k,i,2) = 0.0_wp
1358	dist_xz(k,i,3) = 0.0_wp
1359
1360	ELSE
1361	dist_xz(k,i,1) = a11(k) * fu_xz(k,i)
1362	!experimental test of 1.2
1363	dist_xz(k,i,2) = ( SQRT( a22(k) / MAXVAL(a22) ) &
1364	* 1.2_wp ) &
1365	* ( a21(k) * fu_xz(k,i) &
1366	+ a22(k) * fv_xz(k,i) )
1367	dist_xz(k,i,3) = ( SQRT(a33(k) / MAXVAL(a33) ) &
1368	* 1.3_wp ) &
1369	* ( a31(k) * fu_xz(k,i) &
1370	+ a32(k) * fv_xz(k,i) &
1371	+ a33(k) * fw_xz(k,i) )
1372	ENDIF
1373
1374	ENDDO
1375	ENDDO
1376	!
1377	!-- Mass flux correction following Kim et al. (2013)
1378	!-- This correction factor insures that the mass flux is preserved at the
1379	!-- inflow boundary.
1380	!-- First, calculate volume flow at xz boundary
1381	IF ( myidy == id_stg_south ) j = nys
1382	IF ( myidy == id_stg_north ) j = nyn+1
1383
1384	volume_flow_l = 0.0_wp
1385	volume_flow = 0.0_wp
1386	mc_factor_l = 0.0_wp
1387	mc_factor = 0.0_wp
1388	DO i = nxl, nxr
1389	DO k = nzb+1, nzt
1390	volume_flow_l = volume_flow_l + v(k,j,i) * dzw(k) * dx &
1391	* MERGE( 1.0_wp, 0.0_wp, &
1392	BTEST( wall_flags_0(k,j,i), 2 ) )
1393
1394	mc_factor_l = mc_factor_l + ( v(k,j,i) + dist_xz(k,i,2) ) &
1395	* dzw(k) * dx &
1396	* MERGE( 1.0_wp, 0.0_wp, &
1397	BTEST( wall_flags_0(k,j,i), 2 ) )
1398	ENDDO
1399	ENDDO
1400	#if defined( __parallel )
1401	CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, &
1402	1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1403	CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, &
1404	1, MPI_REAL, MPI_SUM, comm1dx, ierr )
1405	#else
1406	volume_flow = volume_flow_l
1407	mc_factor = mc_factor_l
1408	#endif
1409
1410	mc_factor = volume_flow / mc_factor
1411
1412	!
1413	!-- Add disturbances
1414	IF ( myidy == id_stg_south ) THEN
1415
1416	DO i = nxl, nxr
1417	DO k = nzb, nzt
1418	u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) ) &
1419	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1420	BTEST( wall_flags_0(k,-1,i), 1 ) )
1421	v(k,0,i) = ( v(k,0,i) + dist_xz(k,i,2) ) &
1422	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1423	BTEST( wall_flags_0(k,0,i), 2 ) )
1424	w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3) ) &
1425	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1426	BTEST( wall_flags_0(k,-1,i), 3 ) )
1427	ENDDO
1428	ENDDO
1429	ENDIF
1430	IF ( myidy == id_stg_north ) THEN
1431
1432	DO i = nxl, nxr
1433	DO k = nzb, nzt
1434	u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) ) &
1435	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1436	BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
1437	v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,1) ) &
1438	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1439	BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
1440	w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,1) ) &
1441	* mc_factor * MERGE( 1.0_wp, 0.0_wp, &
1442	BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
1443	ENDDO
1444	ENDDO
1445	ENDIF
1446	ENDIF
1447
1448	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
1449
1450	END SUBROUTINE stg_main
1451
1452	!------------------------------------------------------------------------------!
1453	! Description:
1454	! ------------
1455	!> Generate a set of random number rand_it wich is equal on each PE
1456	!> and calculate the velocity seed f_n.
1457	!> f_n is splitted in vertical direction by the number of PEs in x-direction and
1458	!> and each PE calculates a vertical subsection of f_n. At the the end, all
1459	!> parts are collected to form the full array.
1460	!------------------------------------------------------------------------------!
1461	SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
1462
1463
1464	USE indices, &
1465	ONLY: ny
1466
1467	IMPLICIT NONE
1468
1469	INTEGER(iwp) :: id !< core ids at respective boundaries
1470	INTEGER(iwp) :: j !< loop index in y-direction
1471	INTEGER(iwp) :: jj !< loop index in y-direction
1472	INTEGER(iwp) :: k !< loop index in z-direction
1473	INTEGER(iwp) :: kk !< loop index in z-direction
1474	INTEGER(iwp) :: send_count !< send count for MPI_GATHERV
1475
1476	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y !< length scale in y-direction
1477	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z !< length scale in z-direction
1478	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y2 !< n_y*2
1479	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2 !< n_z*2
1480
1481	REAL(wp) :: nyz_inv !< inverse of number of grid points in yz-slice
1482	REAL(wp) :: rand_av !< average of random number
1483	REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number
1484
1485	REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1) :: b_y !< filter func in y-dir
1486	REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1) :: b_z !< filter func in z-dir
1487	REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l !< local velocity seed
1488	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng) :: f_n !< velocity seed
1489	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rand_it !< random number
1490
1491
1492	!
1493	!-- Generate random numbers using a seed generated in stg_init.
1494	!-- The set of random numbers are modified to have an average of 0 and
1495	!-- unit variance.
1496	ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:ny+mergp) )
1497
1498	rand_av = 0.0_wp
1499	rand_sigma_inv = 0.0_wp
1500	nyz_inv = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( ny+1 ), KIND=wp )
1501
1502	DO j = 0, ny
1503	DO k = nzb, nzt+1
1504	CALL RANDOM_NUMBER( rand_it(k,j) )
1505	rand_av = rand_av + rand_it(k,j)
1506	ENDDO
1507	ENDDO
1508
1509	rand_av = rand_av * nyz_inv
1510
1511	DO j = 0, ny
1512	DO k = nzb, nzt+1
1513	rand_it(k,j) = rand_it(k,j) - rand_av
1514	rand_sigma_inv = rand_sigma_inv + rand_it(k,j) ** 2
1515	ENDDO
1516	ENDDO
1517
1518	rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nyz_inv)
1519
1520	DO j = 0, ny
1521	DO k = nzb, nzt+1
1522	rand_it(k,j) = rand_it(k,j) * rand_sigma_inv
1523	ENDDO
1524	ENDDO
1525
1526	!
1527	!-- Periodic fill of random number in space
1528	DO j = 0, ny
1529	DO k = 1, mergp
1530	rand_it(nzb -k,j) = rand_it(nzt+2-k,j) ! bottom margin
1531	rand_it(nzt+1+k,j) = rand_it(nzb+k-1,j) ! top margin
1532	ENDDO
1533	ENDDO
1534	DO j = 1, mergp
1535	DO k = nzb-mergp, nzt+1+mergp
1536	rand_it(k, -j) = rand_it(k,ny-j+1) ! south margin
1537	rand_it(k,ny+j) = rand_it(k, j-1) ! north margin
1538	ENDDO
1539	ENDDO
1540
1541	!
1542	!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1543	n_y2 = n_y * 2
1544	n_z2 = n_z * 2
1545	f_n_l = 0.0_wp
1546
1547	DO j = nysg, nyng
1548	DO k = nzb_x_stg, nzt_x_stg+1
1549	DO jj = -n_y2(k), n_y2(k)
1550	DO kk = -n_z2(k), n_z2(k)
1551	f_n_l(k,j) = f_n_l(k,j) &
1552	+ b_y(jj,k) * b_z(kk,k) * rand_it(k+kk,j+jj)
1553	ENDDO
1554	ENDDO
1555	ENDDO
1556	ENDDO
1557
1558	DEALLOCATE( rand_it )
1559	!
1560	!-- Gather velocity seeds of full subdomain
1561	send_count = nzt_x_stg - nzb_x_stg + 1
1562	IF ( nzt_x_stg == nzt ) send_count = send_count + 1
1563
1564	#if defined( __parallel )
1565	CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small, &
1566	f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz, &
1567	id, comm1dx, ierr )
1568	#else
1569	f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x_stg:nzt_x_stg+1,nysg:nyng)
1570	#endif
1571
1572
1573	END SUBROUTINE stg_generate_seed_yz
1574
1575
1576
1577
1578	!------------------------------------------------------------------------------!
1579	! Description:
1580	! ------------
1581	!> Generate a set of random number rand_it wich is equal on each PE
1582	!> and calculate the velocity seed f_n.
1583	!> f_n is splitted in vertical direction by the number of PEs in y-direction and
1584	!> and each PE calculates a vertical subsection of f_n. At the the end, all
1585	!> parts are collected to form the full array.
1586	!------------------------------------------------------------------------------!
1587	SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
1588
1589
1590	USE indices, &
1591	ONLY: nx
1592
1593
1594	IMPLICIT NONE
1595
1596	INTEGER(iwp) :: id !< core ids at respective boundaries
1597	INTEGER(iwp) :: i !< loop index in x-direction
1598	INTEGER(iwp) :: ii !< loop index in x-direction
1599	INTEGER(iwp) :: k !< loop index in z-direction
1600	INTEGER(iwp) :: kk !< loop index in z-direction
1601	INTEGER(iwp) :: send_count !< send count for MPI_GATHERV
1602
1603	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x !< length scale in x-direction
1604	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z !< length scale in z-direction
1605	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2 !< n_y*2
1606	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2 !< n_z*2
1607
1608	REAL(wp) :: nxz_inv !< inverse of number of grid points in xz-slice
1609	REAL(wp) :: rand_av !< average of random number
1610	REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number
1611
1612	REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1) :: b_x !< filter func in y-dir
1613	REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1) :: b_z !< filter func in z-dir
1614	REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l !< local velocity seed
1615	REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg) :: f_n !< velocity seed
1616	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rand_it !< random number
1617
1618
1619	!
1620	!-- Generate random numbers using a seed generated in stg_init.
1621	!-- The set of random numbers are modified to have an average of 0 and
1622	!-- unit variance.
1623	ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:nx+mergp) )
1624
1625	rand_av = 0.0_wp
1626	rand_sigma_inv = 0.0_wp
1627	nxz_inv = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( nx+1 ), KIND=wp )
1628
1629	DO i = 0, nx
1630	DO k = nzb, nzt+1
1631	CALL RANDOM_NUMBER( rand_it(k,i) )
1632	rand_av = rand_av + rand_it(k,i)
1633	ENDDO
1634	ENDDO
1635
1636	rand_av = rand_av * nxz_inv
1637
1638	DO i = 0, nx
1639	DO k = nzb, nzt+1
1640	rand_it(k,i) = rand_it(k,i) - rand_av
1641	rand_sigma_inv = rand_sigma_inv + rand_it(k,i) ** 2
1642	ENDDO
1643	ENDDO
1644
1645	rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nxz_inv)
1646
1647	DO i = 0, nx
1648	DO k = nzb, nzt+1
1649	rand_it(k,i) = rand_it(k,i) * rand_sigma_inv
1650	ENDDO
1651	ENDDO
1652
1653	!
1654	!-- Periodic fill of random number in space
1655	DO i = 0, nx
1656	DO k = 1, mergp
1657	rand_it(nzb-k,i) = rand_it(nzt+2-k,i) ! bottom margin
1658	rand_it(nzt+1+k,i) = rand_it(nzb+k-1,i) ! top margin
1659	ENDDO
1660	ENDDO
1661	DO i = 1, mergp
1662	DO k = nzb-mergp, nzt+1+mergp
1663	rand_it(k,-i) = rand_it(k,nx-i+1) ! left margin
1664	rand_it(k,nx+i) = rand_it(k,i-1) ! right margin
1665	ENDDO
1666	ENDDO
1667
1668	!
1669	!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1670	n_x2 = n_x * 2
1671	n_z2 = n_z * 2
1672	f_n_l = 0.0_wp
1673
1674	DO i = nxlg, nxrg
1675	DO k = nzb_y_stg, nzt_y_stg+1
1676	DO ii = -n_x2(k), n_x2(k)
1677	DO kk = -n_z2(k), n_z2(k)
1678	f_n_l(k,i) = f_n_l(k,i) &
1679	+ b_x(ii,k) * b_z(kk,k) * rand_it(k+kk,i+ii)
1680	ENDDO
1681	ENDDO
1682	ENDDO
1683	ENDDO
1684
1685	DEALLOCATE( rand_it )
1686
1687	!
1688	!-- Gather velocity seeds of full subdomain
1689	send_count = nzt_y_stg - nzb_y_stg + 1
1690	IF ( nzt_y_stg == nzt ) send_count = send_count + 1
1691
1692
1693	#if defined( __parallel )
1694	CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small, &
1695	f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz, &
1696	id, comm1dy, ierr )
1697	#else
1698	f_n(nzb+1:nzt+1,nxlg:nxrg) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg)
1699	#endif
1700
1701
1702	END SUBROUTINE stg_generate_seed_xz
1703
1704	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 |