Home

Context Navigation

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

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |