Home

Context Navigation

source: palm/trunk/SOURCE/lpm_advec.f90 @ 2101

Last change on this file since 2101 was 2101, checked in by suehring, 7 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 73.2 KB

Rev	Line
[1682]	1	!> @file lpm_advec.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[2101]	17	! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[849]	20	! Current revisions:
	21	! ------------------
[1930]	22	!
[2101]	23	!
[1930]	24	! Former revisions:
	25	! -----------------
	26	! $Id: lpm_advec.f90 2101 2017-01-05 16:42:31Z suehring $
	27	!
[2101]	28	! 2100 2017-01-05 16:40:16Z suehring
	29	! Prevent extremely large SGS-velocities in regions where TKE is zero, e.g.
	30	! at the begin of simulations and/or in non-turbulent regions.
	31	!
[2001]	32	! 2000 2016-08-20 18:09:15Z knoop
	33	! Forced header and separation lines into 80 columns
	34	!
[1937]	35	! 1936 2016-06-13 13:37:44Z suehring
	36	! Formatting adjustments
	37	!
[1930]	38	! 1929 2016-06-09 16:25:25Z suehring
[1929]	39	! Put stochastic equation in an extra subroutine.
	40	! Set flag for stochastic equation to communicate whether a particle is near
	41	! topography. This case, memory and drift term are disabled in the Weil equation.
[1889]	42	!
[1929]	43	! Enable vertical logarithmic interpolation also above topography. This case,
	44	! set a lower limit for the friction velocity, as it can become very small
[1930]	45	! in narrow street canyons, leading to too large particle speeds.
[1823]	46	!
[1889]	47	! 1888 2016-04-21 12:20:49Z suehring
	48	! Bugfix concerning logarithmic interpolation of particle speed
	49	!
[1823]	50	! 1822 2016-04-07 07:49:42Z hoffmann
[1822]	51	! Random velocity fluctuations for particles added. Terminal fall velocity
	52	! for droplets is calculated from a parameterization (which is better than
	53	! the previous, physically correct calculation, which demands a very short
	54	! time step that is not used in the model).
	55	!
	56	! Unused variables deleted.
[1321]	57	!
[1692]	58	! 1691 2015-10-26 16:17:44Z maronga
	59	! Renamed prandtl_layer to constant_flux_layer.
	60	!
[1686]	61	! 1685 2015-10-08 07:32:13Z raasch
	62	! TKE check for negative values (so far, only zero value was checked)
	63	! offset_ocean_nzt_m1 removed
	64	!
[1683]	65	! 1682 2015-10-07 23:56:08Z knoop
	66	! Code annotations made doxygen readable
	67	!
[1584]	68	! 1583 2015-04-15 12:16:27Z suehring
	69	! Bugfix: particle advection within Prandtl-layer in case of Galilei
	70	! transformation.
	71	!
[1370]	72	! 1369 2014-04-24 05:57:38Z raasch
	73	! usage of module interfaces removed
	74	!
[1360]	75	! 1359 2014-04-11 17:15:14Z hoffmann
	76	! New particle structure integrated.
	77	! Kind definition added to all floating point numbers.
	78	!
[1323]	79	! 1322 2014-03-20 16:38:49Z raasch
	80	! REAL constants defined as wp_kind
	81	!
[1321]	82	! 1320 2014-03-20 08:40:49Z raasch
[1320]	83	! ONLY-attribute added to USE-statements,
	84	! kind-parameters added to all INTEGER and REAL declaration statements,
	85	! kinds are defined in new module kinds,
	86	! revision history before 2012 removed,
	87	! comment fields (!:) to be used for variable explanations added to
	88	! all variable declaration statements
[849]	89	!
[1315]	90	! 1314 2014-03-14 18:25:17Z suehring
	91	! Vertical logarithmic interpolation of horizontal particle speed for particles
	92	! between roughness height and first vertical grid level.
	93	!
[1037]	94	! 1036 2012-10-22 13:43:42Z raasch
	95	! code put under GPL (PALM 3.9)
	96	!
[850]	97	! 849 2012-03-15 10:35:09Z raasch
	98	! initial revision (former part of advec_particles)
[849]	99	!
[850]	100	!
[849]	101	! Description:
	102	! ------------
[1682]	103	!> Calculation of new particle positions due to advection using a simple Euler
	104	!> scheme. Particles may feel inertia effects. SGS transport can be included
	105	!> using the stochastic model of Weil et al. (2004, JAS, 61, 2877-2887).
[849]	106	!------------------------------------------------------------------------------!
[1682]	107	SUBROUTINE lpm_advec (ip,jp,kp)
	108
[849]	109
[1320]	110	USE arrays_3d, &
[1822]	111	ONLY: de_dx, de_dy, de_dz, diss, e, km, u, us, usws, v, vsws, w, zu, zw
[849]	112
[1359]	113	USE cpulog
	114
	115	USE pegrid
	116
[1320]	117	USE control_parameters, &
[1691]	118	ONLY: atmos_ocean_sign, cloud_droplets, constant_flux_layer, dt_3d, &
[1822]	119	dt_3d_reached_l, dz, g, kappa, topography, u_gtrans, v_gtrans
[849]	120
[1320]	121	USE grid_variables, &
	122	ONLY: ddx, dx, ddy, dy
	123
	124	USE indices, &
	125	ONLY: nzb, nzb_s_inner, nzt
	126
	127	USE kinds
	128
	129	USE particle_attributes, &
[1822]	130	ONLY: block_offset, c_0, dt_min_part, grid_particles, &
[1359]	131	iran_part, log_z_z0, number_of_particles, number_of_sublayers, &
[1929]	132	particles, particle_groups, offset_ocean_nzt, sgs_wf_part, &
	133	use_sgs_for_particles, vertical_particle_advection, z0_av_global
[1320]	134
	135	USE statistics, &
	136	ONLY: hom
[849]	137
[1320]	138	IMPLICIT NONE
[849]	139
[1929]	140	INTEGER(iwp) :: agp !< loop variable
	141	INTEGER(iwp) :: gp_outside_of_building(1:8) !< number of grid points used for particle interpolation in case of topography
	142	INTEGER(iwp) :: i !< index variable along x
	143	INTEGER(iwp) :: ip !< index variable along x
	144	INTEGER(iwp) :: ilog !< index variable along x
	145	INTEGER(iwp) :: j !< index variable along y
	146	INTEGER(iwp) :: jp !< index variable along y
	147	INTEGER(iwp) :: jlog !< index variable along y
	148	INTEGER(iwp) :: k !< index variable along z
	149	INTEGER(iwp) :: kp !< index variable along z
	150	INTEGER(iwp) :: kw !< index variable along z
	151	INTEGER(iwp) :: n !< loop variable over all particles in a grid box
	152	INTEGER(iwp) :: nb !< block number particles are sorted in
	153	INTEGER(iwp) :: num_gp !< number of adjacent grid points inside topography
[849]	154
[1929]	155	INTEGER(iwp), DIMENSION(0:7) :: start_index !< start particle index for current block
	156	INTEGER(iwp), DIMENSION(0:7) :: end_index !< start particle index for current block
[1359]	157
[1929]	158	REAL(wp) :: aa !< dummy argument for horizontal particle interpolation
	159	REAL(wp) :: bb !< dummy argument for horizontal particle interpolation
	160	REAL(wp) :: cc !< dummy argument for horizontal particle interpolation
	161	REAL(wp) :: d_sum !< dummy argument for horizontal particle interpolation in case of topography
	162	REAL(wp) :: d_z_p_z0 !< inverse of interpolation length for logarithmic interpolation
	163	REAL(wp) :: dd !< dummy argument for horizontal particle interpolation
	164	REAL(wp) :: de_dx_int_l !< x/y-interpolated TKE gradient (x) at particle position at lower vertical level
	165	REAL(wp) :: de_dx_int_u !< x/y-interpolated TKE gradient (x) at particle position at upper vertical level
	166	REAL(wp) :: de_dy_int_l !< x/y-interpolated TKE gradient (y) at particle position at lower vertical level
	167	REAL(wp) :: de_dy_int_u !< x/y-interpolated TKE gradient (y) at particle position at upper vertical level
	168	REAL(wp) :: de_dt !< temporal derivative of TKE experienced by the particle
	169	REAL(wp) :: de_dt_min !< lower level for temporal TKE derivative
	170	REAL(wp) :: de_dz_int_l !< x/y-interpolated TKE gradient (z) at particle position at lower vertical level
	171	REAL(wp) :: de_dz_int_u !< x/y-interpolated TKE gradient (z) at particle position at upper vertical level
[1822]	172	REAL(wp) :: diameter !< diamter of droplet
[1929]	173	REAL(wp) :: diss_int_l !< x/y-interpolated dissipation at particle position at lower vertical level
	174	REAL(wp) :: diss_int_u !< x/y-interpolated dissipation at particle position at upper vertical level
	175	REAL(wp) :: dt_gap !< remaining time until particle time integration reaches LES time
	176	REAL(wp) :: dt_particle_m !< previous particle time step
	177	REAL(wp) :: e_int_l !< x/y-interpolated TKE at particle position at lower vertical level
	178	REAL(wp) :: e_int_u !< x/y-interpolated TKE at particle position at upper vertical level
	179	REAL(wp) :: e_mean_int !< horizontal mean TKE at particle height
[1682]	180	REAL(wp) :: exp_arg !<
	181	REAL(wp) :: exp_term !<
[1929]	182	REAL(wp) :: gg !< dummy argument for horizontal particle interpolation
	183	REAL(wp) :: height_p !< dummy argument for logarithmic interpolation
[1822]	184	REAL(wp) :: lagr_timescale !< Lagrangian timescale
[1929]	185	REAL(wp) :: location(1:30,1:3) !< wall locations
	186	REAL(wp) :: log_z_z0_int !< logarithmus used for surface_layer interpolation
[1682]	187	REAL(wp) :: random_gauss !<
[1822]	188	REAL(wp) :: RL !< Lagrangian autocorrelation coefficient
	189	REAL(wp) :: rg1 !< Gaussian distributed random number
	190	REAL(wp) :: rg2 !< Gaussian distributed random number
	191	REAL(wp) :: rg3 !< Gaussian distributed random number
	192	REAL(wp) :: sigma !< velocity standard deviation
[1929]	193	REAL(wp) :: u_int_l !< x/y-interpolated u-component at particle position at lower vertical level
	194	REAL(wp) :: u_int_u !< x/y-interpolated u-component at particle position at upper vertical level
	195	REAL(wp) :: us_int !< friction velocity at particle grid box
	196	REAL(wp) :: v_int_l !< x/y-interpolated v-component at particle position at lower vertical level
	197	REAL(wp) :: v_int_u !< x/y-interpolated v-component at particle position at upper vertical level
[1682]	198	REAL(wp) :: vv_int !<
[1929]	199	REAL(wp) :: w_int_l !< x/y-interpolated w-component at particle position at lower vertical level
	200	REAL(wp) :: w_int_u !< x/y-interpolated w-component at particle position at upper vertical level
[1822]	201	REAL(wp) :: w_s !< terminal velocity of droplets
[1929]	202	REAL(wp) :: x !< dummy argument for horizontal particle interpolation
	203	REAL(wp) :: y !< dummy argument for horizontal particle interpolation
	204	REAL(wp) :: z_p !< surface layer height (0.5 dz)
[849]	205
[1822]	206	REAL(wp), PARAMETER :: a_rog = 9.65_wp !< parameter for fall velocity
	207	REAL(wp), PARAMETER :: b_rog = 10.43_wp !< parameter for fall velocity
	208	REAL(wp), PARAMETER :: c_rog = 0.6_wp !< parameter for fall velocity
	209	REAL(wp), PARAMETER :: k_cap_rog = 4.0_wp !< parameter for fall velocity
	210	REAL(wp), PARAMETER :: k_low_rog = 12.0_wp !< parameter for fall velocity
	211	REAL(wp), PARAMETER :: d0_rog = 0.745_wp !< separation diameter
	212
[1929]	213	REAL(wp), DIMENSION(1:30) :: d_gp_pl !< dummy argument for particle interpolation scheme in case of topography
	214	REAL(wp), DIMENSION(1:30) :: de_dxi !< horizontal TKE gradient along x at adjacent wall
	215	REAL(wp), DIMENSION(1:30) :: de_dyi !< horizontal TKE gradient along y at adjacent wall
	216	REAL(wp), DIMENSION(1:30) :: de_dzi !< horizontal TKE gradient along z at adjacent wall
	217	REAL(wp), DIMENSION(1:30) :: dissi !< dissipation at adjacent wall
	218	REAL(wp), DIMENSION(1:30) :: ei !< TKE at adjacent wall
[849]	219
[1929]	220	REAL(wp), DIMENSION(number_of_particles) :: term_1_2 !< flag to communicate whether a particle is near topography or not
[1682]	221	REAL(wp), DIMENSION(number_of_particles) :: dens_ratio !<
[1929]	222	REAL(wp), DIMENSION(number_of_particles) :: de_dx_int !< horizontal TKE gradient along x at particle position
	223	REAL(wp), DIMENSION(number_of_particles) :: de_dy_int !< horizontal TKE gradient along y at particle position
	224	REAL(wp), DIMENSION(number_of_particles) :: de_dz_int !< horizontal TKE gradient along z at particle position
	225	REAL(wp), DIMENSION(number_of_particles) :: diss_int !< dissipation at particle position
	226	REAL(wp), DIMENSION(number_of_particles) :: dt_particle !< particle time step
	227	REAL(wp), DIMENSION(number_of_particles) :: e_int !< TKE at particle position
	228	REAL(wp), DIMENSION(number_of_particles) :: fs_int !< weighting factor for subgrid-scale particle speed
	229	REAL(wp), DIMENSION(number_of_particles) :: u_int !< u-component of particle speed
	230	REAL(wp), DIMENSION(number_of_particles) :: v_int !< v-component of particle speed
	231	REAL(wp), DIMENSION(number_of_particles) :: w_int !< w-component of particle speed
	232	REAL(wp), DIMENSION(number_of_particles) :: xv !< x-position
	233	REAL(wp), DIMENSION(number_of_particles) :: yv !< y-position
	234	REAL(wp), DIMENSION(number_of_particles) :: zv !< z-position
[1359]	235
[1929]	236	REAL(wp), DIMENSION(number_of_particles, 3) :: rg !< vector of Gaussian distributed random numbers
[1359]	237
	238	CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
	239
[1314]	240	!
	241	!-- Determine height of Prandtl layer and distance between Prandtl-layer
	242	!-- height and horizontal mean roughness height, which are required for
	243	!-- vertical logarithmic interpolation of horizontal particle speeds
	244	!-- (for particles below first vertical grid level).
	245	z_p = zu(nzb+1) - zw(nzb)
[1359]	246	d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
[849]	247
[1359]	248	start_index = grid_particles(kp,jp,ip)%start_index
	249	end_index = grid_particles(kp,jp,ip)%end_index
[849]	250
[1359]	251	xv = particles(1:number_of_particles)%x
	252	yv = particles(1:number_of_particles)%y
	253	zv = particles(1:number_of_particles)%z
[849]	254
[1359]	255	DO nb = 0, 7
[1314]	256
[1359]	257	i = ip
	258	j = jp + block_offset(nb)%j_off
	259	k = kp + block_offset(nb)%k_off
	260
[1929]	261
[849]	262	!
[1359]	263	!-- Interpolate u velocity-component
	264	DO n = start_index(nb), end_index(nb)
[1314]	265	!
[1359]	266	!-- Interpolation of the u velocity component onto particle position.
	267	!-- Particles are interpolation bi-linearly in the horizontal and a
	268	!-- linearly in the vertical. An exception is made for particles below
	269	!-- the first vertical grid level in case of a prandtl layer. In this
	270	!-- case the horizontal particle velocity components are determined using
	271	!-- Monin-Obukhov relations (if branch).
	272	!-- First, check if particle is located below first vertical grid level
	273	!-- (Prandtl-layer height)
[1929]	274	ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
	275	jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
	276
[1936]	277	IF ( constant_flux_layer .AND. &
	278	zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p ) THEN
[1314]	279	!
[1359]	280	!-- Resolved-scale horizontal particle velocity is zero below z0.
[1929]	281	IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global ) THEN
[1359]	282	u_int(n) = 0.0_wp
	283	ELSE
[1314]	284	!
[1359]	285	!-- Determine the sublayer. Further used as index.
[1936]	286	height_p = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global ) &
	287	* REAL( number_of_sublayers, KIND=wp ) &
[1359]	288	* d_z_p_z0
[1314]	289	!
[1359]	290	!-- Calculate LOG(z/z0) for exact particle height. Therefore,
	291	!-- interpolate linearly between precalculated logarithm.
[1929]	292	log_z_z0_int = log_z_z0(INT(height_p)) &
[1359]	293	+ ( height_p - INT(height_p) ) &
	294	* ( log_z_z0(INT(height_p)+1) &
	295	- log_z_z0(INT(height_p)) &
	296	)
[1314]	297	!
[1929]	298	!-- Limit friction velocity. In narrow canyons or holes the
	299	!-- friction velocity can become very small, resulting in a too
	300	!-- large particle speed.
	301	us_int = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog,ilog-1) ), &
	302	0.01_wp )
	303	!
[1359]	304	!-- Neutral solution is applied for all situations, e.g. also for
	305	!-- unstable and stable situations. Even though this is not exact
	306	!-- this saves a lot of CPU time since several calls of intrinsic
	307	!-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified
	308	!-- as sensitivity studies revealed no significant effect of
	309	!-- using the neutral solution also for un/stable situations.
[1929]	310	u_int(n) = -usws(jlog,ilog) / ( us_int * kappa + 1E-10_wp ) &
	311	* log_z_z0_int - u_gtrans
	312
[1359]	313	ENDIF
	314	!
	315	!-- Particle above the first grid level. Bi-linear interpolation in the
	316	!-- horizontal and linear interpolation in the vertical direction.
[1314]	317	ELSE
	318
[1359]	319	x = xv(n) + ( 0.5_wp - i ) * dx
	320	y = yv(n) - j * dy
	321	aa = x2 + y2
	322	bb = ( dx - x )2 + y2
	323	cc = x2 + ( dy - y )2
	324	dd = ( dx - x )2 + ( dy - y )2
	325	gg = aa + bb + cc + dd
[1314]	326
[1359]	327	u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) &
	328	+ ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * &
	329	u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
[1314]	330
[1359]	331	IF ( k == nzt ) THEN
	332	u_int(n) = u_int_l
	333	ELSE
	334	u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) &
	335	+ ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * &
	336	u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
	337	u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dz * &
	338	( u_int_u - u_int_l )
	339	ENDIF
[1929]	340
[1314]	341	ENDIF
	342
[1359]	343	ENDDO
[849]	344
[1359]	345	i = ip + block_offset(nb)%i_off
	346	j = jp
	347	k = kp + block_offset(nb)%k_off
[849]	348	!
[1359]	349	!-- Same procedure for interpolation of the v velocity-component
	350	DO n = start_index(nb), end_index(nb)
[1685]	351
[1929]	352	ilog = ( particles(n)%x + 0.5_wp * dx ) * ddx
	353	jlog = ( particles(n)%y + 0.5_wp * dy ) * ddy
[1936]	354	IF ( constant_flux_layer .AND. &
	355	zv(n) - zw(nzb_s_inner(jlog,ilog)) < z_p ) THEN
[849]	356
[1929]	357	IF ( zv(n) - zw(nzb_s_inner(jlog,ilog)) < z0_av_global ) THEN
[1314]	358	!
[1359]	359	!-- Resolved-scale horizontal particle velocity is zero below z0.
	360	v_int(n) = 0.0_wp
	361	ELSE
	362	!
[1929]	363	!-- Determine the sublayer. Further used as index. Please note,
	364	!-- logarithmus can not be reused from above, as in in case of
	365	!-- topography particle on u-grid can be above surface-layer height,
	366	!-- whereas it can be below on v-grid.
[1936]	367	height_p = ( zv(n) - zw(nzb_s_inner(jlog,ilog)) - z0_av_global ) &
	368	* REAL( number_of_sublayers, KIND=wp ) &
[1929]	369	* d_z_p_z0
	370	!
	371	!-- Calculate LOG(z/z0) for exact particle height. Therefore,
	372	!-- interpolate linearly between precalculated logarithm.
	373	log_z_z0_int = log_z_z0(INT(height_p)) &
	374	+ ( height_p - INT(height_p) ) &
	375	* ( log_z_z0(INT(height_p)+1) &
	376	- log_z_z0(INT(height_p)) &
	377	)
	378	!
	379	!-- Limit friction velocity. In narrow canyons or holes the
	380	!-- friction velocity can become very small, resulting in a too
	381	!-- large particle speed.
[1936]	382	us_int = MAX( 0.5_wp * ( us(jlog,ilog) + us(jlog-1,ilog) ), &
[1929]	383	0.01_wp )
	384	!
[1359]	385	!-- Neutral solution is applied for all situations, e.g. also for
	386	!-- unstable and stable situations. Even though this is not exact
	387	!-- this saves a lot of CPU time since several calls of intrinsic
	388	!-- FORTRAN procedures (LOG, ATAN) are avoided, This is justified
	389	!-- as sensitivity studies revealed no significant effect of
	390	!-- using the neutral solution also for un/stable situations.
[1936]	391	v_int(n) = -vsws(jlog,ilog) / ( us_int * kappa + 1E-10_wp ) &
[1929]	392	* log_z_z0_int - v_gtrans
[1314]	393
[1359]	394	ENDIF
[1929]	395
[1359]	396	ELSE
	397	x = xv(n) - i * dx
	398	y = yv(n) + ( 0.5_wp - j ) * dy
	399	aa = x2 + y2
	400	bb = ( dx - x )2 + y2
	401	cc = x2 + ( dy - y )2
	402	dd = ( dx - x )2 + ( dy - y )2
	403	gg = aa + bb + cc + dd
[1314]	404
[1359]	405	v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) &
	406	+ ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
	407	) / ( 3.0_wp * gg ) - v_gtrans
[1314]	408
[1359]	409	IF ( k == nzt ) THEN
	410	v_int(n) = v_int_l
	411	ELSE
	412	v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) &
	413	+ ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
	414	) / ( 3.0_wp * gg ) - v_gtrans
	415	v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dz * &
	416	( v_int_u - v_int_l )
	417	ENDIF
[1929]	418
[1314]	419	ENDIF
	420
[1359]	421	ENDDO
[1314]	422
[1359]	423	i = ip + block_offset(nb)%i_off
	424	j = jp + block_offset(nb)%j_off
[1929]	425	k = kp - 1
[849]	426	!
[1314]	427	!-- Same procedure for interpolation of the w velocity-component
[1359]	428	DO n = start_index(nb), end_index(nb)
[849]	429
[1359]	430	IF ( vertical_particle_advection(particles(n)%group) ) THEN
[849]	431
[1359]	432	x = xv(n) - i * dx
	433	y = yv(n) - j * dy
[849]	434	aa = x2 + y2
	435	bb = ( dx - x )2 + y2
	436	cc = x2 + ( dy - y )2
	437	dd = ( dx - x )2 + ( dy - y )2
	438	gg = aa + bb + cc + dd
	439
[1359]	440	w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) &
	441	+ ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
	442	) / ( 3.0_wp * gg )
[849]	443
[1359]	444	IF ( k == nzt ) THEN
	445	w_int(n) = w_int_l
[849]	446	ELSE
[1359]	447	w_int_u = ( ( gg-aa ) * w(k+1,j,i) + &
	448	( gg-bb ) * w(k+1,j,i+1) + &
	449	( gg-cc ) * w(k+1,j+1,i) + &
	450	( gg-dd ) * w(k+1,j+1,i+1) &
	451	) / ( 3.0_wp * gg )
	452	w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dz * &
	453	( w_int_u - w_int_l )
[849]	454	ENDIF
	455
[1359]	456	ELSE
[849]	457
[1359]	458	w_int(n) = 0.0_wp
[849]	459
[1359]	460	ENDIF
	461
	462	ENDDO
	463
	464	ENDDO
	465
	466	!-- Interpolate and calculate quantities needed for calculating the SGS
	467	!-- velocities
[1822]	468	IF ( use_sgs_for_particles .AND. .NOT. cloud_droplets ) THEN
[1359]	469
	470	IF ( topography == 'flat' ) THEN
	471
	472	DO nb = 0,7
	473
	474	i = ip + block_offset(nb)%i_off
	475	j = jp + block_offset(nb)%j_off
	476	k = kp + block_offset(nb)%k_off
	477
	478	DO n = start_index(nb), end_index(nb)
[849]	479	!
[1359]	480	!-- Interpolate TKE
	481	x = xv(n) - i * dx
	482	y = yv(n) - j * dy
	483	aa = x2 + y2
	484	bb = ( dx - x )2 + y2
	485	cc = x2 + ( dy - y )2
	486	dd = ( dx - x )2 + ( dy - y )2
	487	gg = aa + bb + cc + dd
[849]	488
[1359]	489	e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
	490	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
	491	) / ( 3.0_wp * gg )
	492
	493	IF ( k+1 == nzt+1 ) THEN
	494	e_int(n) = e_int_l
	495	ELSE
	496	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
	497	( gg - bb ) * e(k+1,j,i+1) + &
	498	( gg - cc ) * e(k+1,j+1,i) + &
	499	( gg - dd ) * e(k+1,j+1,i+1) &
	500	) / ( 3.0_wp * gg )
	501	e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
	502	( e_int_u - e_int_l )
	503	ENDIF
[849]	504	!
[1685]	505	!-- Needed to avoid NaN particle velocities (this might not be
	506	!-- required any more)
	507	IF ( e_int(n) <= 0.0_wp ) THEN
[1359]	508	e_int(n) = 1.0E-20_wp
	509	ENDIF
	510	!
	511	!-- Interpolate the TKE gradient along x (adopt incides i,j,k and
	512	!-- all position variables from above (TKE))
	513	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
	514	( gg - bb ) * de_dx(k,j,i+1) + &
	515	( gg - cc ) * de_dx(k,j+1,i) + &
	516	( gg - dd ) * de_dx(k,j+1,i+1) &
	517	) / ( 3.0_wp * gg )
[849]	518
	519	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
[1359]	520	de_dx_int(n) = de_dx_int_l
[849]	521	ELSE
[1359]	522	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
	523	( gg - bb ) * de_dx(k+1,j,i+1) + &
	524	( gg - cc ) * de_dx(k+1,j+1,i) + &
	525	( gg - dd ) * de_dx(k+1,j+1,i+1) &
	526	) / ( 3.0_wp * gg )
	527	de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / dz * &
	528	( de_dx_int_u - de_dx_int_l )
[849]	529	ENDIF
[1359]	530	!
	531	!-- Interpolate the TKE gradient along y
	532	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
	533	( gg - bb ) * de_dy(k,j,i+1) + &
	534	( gg - cc ) * de_dy(k,j+1,i) + &
	535	( gg - dd ) * de_dy(k,j+1,i+1) &
	536	) / ( 3.0_wp * gg )
	537	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	538	de_dy_int(n) = de_dy_int_l
	539	ELSE
	540	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
	541	( gg - bb ) * de_dy(k+1,j,i+1) + &
	542	( gg - cc ) * de_dy(k+1,j+1,i) + &
	543	( gg - dd ) * de_dy(k+1,j+1,i+1) &
	544	) / ( 3.0_wp * gg )
	545	de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
	546	( de_dy_int_u - de_dy_int_l )
	547	ENDIF
[849]	548
	549	!
[1359]	550	!-- Interpolate the TKE gradient along z
	551	IF ( zv(n) < 0.5_wp * dz ) THEN
	552	de_dz_int(n) = 0.0_wp
	553	ELSE
	554	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
	555	( gg - bb ) * de_dz(k,j,i+1) + &
	556	( gg - cc ) * de_dz(k,j+1,i) + &
	557	( gg - dd ) * de_dz(k,j+1,i+1) &
	558	) / ( 3.0_wp * gg )
[849]	559
[1359]	560	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	561	de_dz_int(n) = de_dz_int_l
	562	ELSE
	563	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
	564	( gg - bb ) * de_dz(k+1,j,i+1) + &
	565	( gg - cc ) * de_dz(k+1,j+1,i) + &
	566	( gg - dd ) * de_dz(k+1,j+1,i+1) &
	567	) / ( 3.0_wp * gg )
	568	de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) / dz * &
	569	( de_dz_int_u - de_dz_int_l )
	570	ENDIF
	571	ENDIF
[849]	572
[1359]	573	!
	574	!-- Interpolate the dissipation of TKE
	575	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
	576	( gg - bb ) * diss(k,j,i+1) + &
	577	( gg - cc ) * diss(k,j+1,i) + &
	578	( gg - dd ) * diss(k,j+1,i+1) &
	579	) / ( 3.0_wp * gg )
[849]	580
[1359]	581	IF ( k == nzt ) THEN
	582	diss_int(n) = diss_int_l
	583	ELSE
	584	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
	585	( gg - bb ) * diss(k+1,j,i+1) + &
	586	( gg - cc ) * diss(k+1,j+1,i) + &
	587	( gg - dd ) * diss(k+1,j+1,i+1) &
	588	) / ( 3.0_wp * gg )
	589	diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
	590	( diss_int_u - diss_int_l )
	591	ENDIF
	592
[1929]	593	!
	594	!-- Set flag for stochastic equation.
	595	term_1_2(n) = 1.0_wp
	596
[1359]	597	ENDDO
	598	ENDDO
	599
	600	ELSE ! non-flat topography, e.g., buildings
	601
	602	DO n = 1, number_of_particles
	603	i = particles(n)%x * ddx
	604	j = particles(n)%y * ddy
	605	k = ( zv(n) + 0.5_wp * dz * atmos_ocean_sign ) / dz &
	606	+ offset_ocean_nzt ! only exact if eq.dist
[849]	607	!
	608	!-- In case that there are buildings it has to be determined
	609	!-- how many of the gridpoints defining the particle box are
	610	!-- situated within a building
	611	!-- gp_outside_of_building(1): i,j,k
	612	!-- gp_outside_of_building(2): i,j+1,k
	613	!-- gp_outside_of_building(3): i,j,k+1
	614	!-- gp_outside_of_building(4): i,j+1,k+1
	615	!-- gp_outside_of_building(5): i+1,j,k
	616	!-- gp_outside_of_building(6): i+1,j+1,k
	617	!-- gp_outside_of_building(7): i+1,j,k+1
	618	!-- gp_outside_of_building(8): i+1,j+1,k+1
	619
	620	gp_outside_of_building = 0
[1359]	621	location = 0.0_wp
[849]	622	num_gp = 0
	623
	624	IF ( k > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
	625	num_gp = num_gp + 1
	626	gp_outside_of_building(1) = 1
	627	location(num_gp,1) = i * dx
	628	location(num_gp,2) = j * dy
[1359]	629	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	630	ei(num_gp) = e(k,j,i)
	631	dissi(num_gp) = diss(k,j,i)
	632	de_dxi(num_gp) = de_dx(k,j,i)
	633	de_dyi(num_gp) = de_dy(k,j,i)
	634	de_dzi(num_gp) = de_dz(k,j,i)
	635	ENDIF
[1929]	636	IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) THEN
[849]	637	num_gp = num_gp + 1
	638	gp_outside_of_building(2) = 1
	639	location(num_gp,1) = i * dx
	640	location(num_gp,2) = (j+1) * dy
[1359]	641	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	642	ei(num_gp) = e(k,j+1,i)
	643	dissi(num_gp) = diss(k,j+1,i)
	644	de_dxi(num_gp) = de_dx(k,j+1,i)
	645	de_dyi(num_gp) = de_dy(k,j+1,i)
	646	de_dzi(num_gp) = de_dz(k,j+1,i)
	647	ENDIF
	648
	649	IF ( k+1 > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
	650	num_gp = num_gp + 1
	651	gp_outside_of_building(3) = 1
	652	location(num_gp,1) = i * dx
	653	location(num_gp,2) = j * dy
[1359]	654	location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
[849]	655	ei(num_gp) = e(k+1,j,i)
	656	dissi(num_gp) = diss(k+1,j,i)
	657	de_dxi(num_gp) = de_dx(k+1,j,i)
	658	de_dyi(num_gp) = de_dy(k+1,j,i)
	659	de_dzi(num_gp) = de_dz(k+1,j,i)
	660	ENDIF
	661
[1929]	662	IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) THEN
[849]	663	num_gp = num_gp + 1
	664	gp_outside_of_building(4) = 1
	665	location(num_gp,1) = i * dx
	666	location(num_gp,2) = (j+1) * dy
[1359]	667	location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
[849]	668	ei(num_gp) = e(k+1,j+1,i)
	669	dissi(num_gp) = diss(k+1,j+1,i)
	670	de_dxi(num_gp) = de_dx(k+1,j+1,i)
	671	de_dyi(num_gp) = de_dy(k+1,j+1,i)
	672	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	673	ENDIF
	674
[1929]	675	IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) THEN
[849]	676	num_gp = num_gp + 1
	677	gp_outside_of_building(5) = 1
	678	location(num_gp,1) = (i+1) * dx
	679	location(num_gp,2) = j * dy
[1359]	680	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	681	ei(num_gp) = e(k,j,i+1)
	682	dissi(num_gp) = diss(k,j,i+1)
	683	de_dxi(num_gp) = de_dx(k,j,i+1)
	684	de_dyi(num_gp) = de_dy(k,j,i+1)
	685	de_dzi(num_gp) = de_dz(k,j,i+1)
	686	ENDIF
	687
[1929]	688	IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) THEN
[849]	689	num_gp = num_gp + 1
	690	gp_outside_of_building(6) = 1
	691	location(num_gp,1) = (i+1) * dx
	692	location(num_gp,2) = (j+1) * dy
[1359]	693	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	694	ei(num_gp) = e(k,j+1,i+1)
	695	dissi(num_gp) = diss(k,j+1,i+1)
	696	de_dxi(num_gp) = de_dx(k,j+1,i+1)
	697	de_dyi(num_gp) = de_dy(k,j+1,i+1)
	698	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	699	ENDIF
	700
[1929]	701	IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) THEN
[849]	702	num_gp = num_gp + 1
	703	gp_outside_of_building(7) = 1
	704	location(num_gp,1) = (i+1) * dx
	705	location(num_gp,2) = j * dy
[1359]	706	location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
[849]	707	ei(num_gp) = e(k+1,j,i+1)
	708	dissi(num_gp) = diss(k+1,j,i+1)
	709	de_dxi(num_gp) = de_dx(k+1,j,i+1)
	710	de_dyi(num_gp) = de_dy(k+1,j,i+1)
	711	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	712	ENDIF
	713
[1929]	714	IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0) THEN
[849]	715	num_gp = num_gp + 1
	716	gp_outside_of_building(8) = 1
	717	location(num_gp,1) = (i+1) * dx
	718	location(num_gp,2) = (j+1) * dy
[1359]	719	location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
[849]	720	ei(num_gp) = e(k+1,j+1,i+1)
	721	dissi(num_gp) = diss(k+1,j+1,i+1)
	722	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
	723	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
	724	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	725	ENDIF
	726	!
	727	!-- If all gridpoints are situated outside of a building, then the
	728	!-- ordinary interpolation scheme can be used.
	729	IF ( num_gp == 8 ) THEN
	730
	731	x = particles(n)%x - i * dx
	732	y = particles(n)%y - j * dy
	733	aa = x2 + y2
	734	bb = ( dx - x )2 + y2
	735	cc = x2 + ( dy - y )2
	736	dd = ( dx - x )2 + ( dy - y )2
	737	gg = aa + bb + cc + dd
[1929]	738
[1359]	739	e_int_l = ( ( gg - aa ) * e(k,j,i) + ( gg - bb ) * e(k,j,i+1) &
	740	+ ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
	741	) / ( 3.0_wp * gg )
[1929]	742
[1359]	743	IF ( k == nzt ) THEN
	744	e_int(n) = e_int_l
[849]	745	ELSE
	746	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
	747	( gg - bb ) * e(k+1,j,i+1) + &
	748	( gg - cc ) * e(k+1,j+1,i) + &
	749	( gg - dd ) * e(k+1,j+1,i+1) &
[1359]	750	) / ( 3.0_wp * gg )
	751	e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
[1929]	752	( e_int_u - e_int_l )
[849]	753	ENDIF
[1929]	754	!
[1685]	755	!-- Needed to avoid NaN particle velocities (this might not be
	756	!-- required any more)
	757	IF ( e_int(n) <= 0.0_wp ) THEN
[1359]	758	e_int(n) = 1.0E-20_wp
	759	ENDIF
	760	!
[849]	761	!-- Interpolate the TKE gradient along x (adopt incides i,j,k
	762	!-- and all position variables from above (TKE))
	763	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
	764	( gg - bb ) * de_dx(k,j,i+1) + &
	765	( gg - cc ) * de_dx(k,j+1,i) + &
	766	( gg - dd ) * de_dx(k,j+1,i+1) &
[1359]	767	) / ( 3.0_wp * gg )
[849]	768
[1359]	769	IF ( ( k == nzt ) .OR. ( k == nzb ) ) THEN
	770	de_dx_int(n) = de_dx_int_l
[849]	771	ELSE
	772	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
	773	( gg - bb ) * de_dx(k+1,j,i+1) + &
	774	( gg - cc ) * de_dx(k+1,j+1,i) + &
	775	( gg - dd ) * de_dx(k+1,j+1,i+1) &
[1359]	776	) / ( 3.0_wp * gg )
	777	de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
[849]	778	dz * ( de_dx_int_u - de_dx_int_l )
	779	ENDIF
	780
	781	!
	782	!-- Interpolate the TKE gradient along y
	783	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
	784	( gg - bb ) * de_dy(k,j,i+1) + &
	785	( gg - cc ) * de_dy(k,j+1,i) + &
	786	( gg - dd ) * de_dy(k,j+1,i+1) &
[1359]	787	) / ( 3.0_wp * gg )
[849]	788	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
[1359]	789	de_dy_int(n) = de_dy_int_l
[849]	790	ELSE
	791	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
	792	( gg - bb ) * de_dy(k+1,j,i+1) + &
	793	( gg - cc ) * de_dy(k+1,j+1,i) + &
	794	( gg - dd ) * de_dy(k+1,j+1,i+1) &
[1359]	795	) / ( 3.0_wp * gg )
	796	de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
[849]	797	dz * ( de_dy_int_u - de_dy_int_l )
	798	ENDIF
	799
	800	!
	801	!-- Interpolate the TKE gradient along z
[1359]	802	IF ( zv(n) < 0.5_wp * dz ) THEN
	803	de_dz_int(n) = 0.0_wp
[849]	804	ELSE
	805	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
	806	( gg - bb ) * de_dz(k,j,i+1) + &
	807	( gg - cc ) * de_dz(k,j+1,i) + &
	808	( gg - dd ) * de_dz(k,j+1,i+1) &
[1359]	809	) / ( 3.0_wp * gg )
[1929]	810
[849]	811	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
[1359]	812	de_dz_int(n) = de_dz_int_l
[849]	813	ELSE
	814	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
	815	( gg - bb ) * de_dz(k+1,j,i+1) + &
	816	( gg - cc ) * de_dz(k+1,j+1,i) + &
	817	( gg - dd ) * de_dz(k+1,j+1,i+1) &
[1359]	818	) / ( 3.0_wp * gg )
	819	de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
[849]	820	dz * ( de_dz_int_u - de_dz_int_l )
	821	ENDIF
	822	ENDIF
	823
	824	!
	825	!-- Interpolate the dissipation of TKE
	826	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
	827	( gg - bb ) * diss(k,j,i+1) + &
	828	( gg - cc ) * diss(k,j+1,i) + &
	829	( gg - dd ) * diss(k,j+1,i+1) &
[1359]	830	) / ( 3.0_wp * gg )
[1929]	831
[1359]	832	IF ( k == nzt ) THEN
	833	diss_int(n) = diss_int_l
[849]	834	ELSE
	835	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
	836	( gg - bb ) * diss(k+1,j,i+1) + &
	837	( gg - cc ) * diss(k+1,j+1,i) + &
	838	( gg - dd ) * diss(k+1,j+1,i+1) &
[1359]	839	) / ( 3.0_wp * gg )
	840	diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
[849]	841	( diss_int_u - diss_int_l )
	842	ENDIF
[1929]	843	!
	844	!-- Set flag for stochastic equation.
	845	term_1_2(n) = 1.0_wp
	846
[849]	847	ELSE
[1929]	848
[849]	849	!
	850	!-- If wall between gridpoint 1 and gridpoint 5, then
	851	!-- Neumann boundary condition has to be applied
	852	IF ( gp_outside_of_building(1) == 1 .AND. &
	853	gp_outside_of_building(5) == 0 ) THEN
	854	num_gp = num_gp + 1
[1359]	855	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	856	location(num_gp,2) = j * dy
[1359]	857	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	858	ei(num_gp) = e(k,j,i)
	859	dissi(num_gp) = diss(k,j,i)
[1359]	860	de_dxi(num_gp) = 0.0_wp
[849]	861	de_dyi(num_gp) = de_dy(k,j,i)
	862	de_dzi(num_gp) = de_dz(k,j,i)
	863	ENDIF
	864
	865	IF ( gp_outside_of_building(5) == 1 .AND. &
[1929]	866	gp_outside_of_building(1) == 0 ) THEN
[849]	867	num_gp = num_gp + 1
[1359]	868	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	869	location(num_gp,2) = j * dy
[1359]	870	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	871	ei(num_gp) = e(k,j,i+1)
	872	dissi(num_gp) = diss(k,j,i+1)
[1359]	873	de_dxi(num_gp) = 0.0_wp
[849]	874	de_dyi(num_gp) = de_dy(k,j,i+1)
	875	de_dzi(num_gp) = de_dz(k,j,i+1)
	876	ENDIF
	877
	878	!
	879	!-- If wall between gridpoint 5 and gridpoint 6, then
	880	!-- then Neumann boundary condition has to be applied
	881	IF ( gp_outside_of_building(5) == 1 .AND. &
	882	gp_outside_of_building(6) == 0 ) THEN
	883	num_gp = num_gp + 1
	884	location(num_gp,1) = (i+1) * dx
[1359]	885	location(num_gp,2) = j * dy + 0.5_wp * dy
	886	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	887	ei(num_gp) = e(k,j,i+1)
	888	dissi(num_gp) = diss(k,j,i+1)
	889	de_dxi(num_gp) = de_dx(k,j,i+1)
[1359]	890	de_dyi(num_gp) = 0.0_wp
[849]	891	de_dzi(num_gp) = de_dz(k,j,i+1)
	892	ENDIF
	893
	894	IF ( gp_outside_of_building(6) == 1 .AND. &
	895	gp_outside_of_building(5) == 0 ) THEN
	896	num_gp = num_gp + 1
	897	location(num_gp,1) = (i+1) * dx
[1359]	898	location(num_gp,2) = j * dy + 0.5_wp * dy
	899	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	900	ei(num_gp) = e(k,j+1,i+1)
	901	dissi(num_gp) = diss(k,j+1,i+1)
	902	de_dxi(num_gp) = de_dx(k,j+1,i+1)
[1359]	903	de_dyi(num_gp) = 0.0_wp
[849]	904	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	905	ENDIF
	906
	907	!
	908	!-- If wall between gridpoint 2 and gridpoint 6, then
	909	!-- Neumann boundary condition has to be applied
	910	IF ( gp_outside_of_building(2) == 1 .AND. &
	911	gp_outside_of_building(6) == 0 ) THEN
	912	num_gp = num_gp + 1
[1359]	913	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	914	location(num_gp,2) = (j+1) * dy
[1359]	915	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	916	ei(num_gp) = e(k,j+1,i)
	917	dissi(num_gp) = diss(k,j+1,i)
[1359]	918	de_dxi(num_gp) = 0.0_wp
[849]	919	de_dyi(num_gp) = de_dy(k,j+1,i)
	920	de_dzi(num_gp) = de_dz(k,j+1,i)
	921	ENDIF
	922
	923	IF ( gp_outside_of_building(6) == 1 .AND. &
	924	gp_outside_of_building(2) == 0 ) THEN
	925	num_gp = num_gp + 1
[1359]	926	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	927	location(num_gp,2) = (j+1) * dy
[1359]	928	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	929	ei(num_gp) = e(k,j+1,i+1)
	930	dissi(num_gp) = diss(k,j+1,i+1)
[1359]	931	de_dxi(num_gp) = 0.0_wp
[849]	932	de_dyi(num_gp) = de_dy(k,j+1,i+1)
	933	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	934	ENDIF
	935
	936	!
	937	!-- If wall between gridpoint 1 and gridpoint 2, then
	938	!-- Neumann boundary condition has to be applied
	939	IF ( gp_outside_of_building(1) == 1 .AND. &
	940	gp_outside_of_building(2) == 0 ) THEN
	941	num_gp = num_gp + 1
	942	location(num_gp,1) = i * dx
[1359]	943	location(num_gp,2) = j * dy + 0.5_wp * dy
	944	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	945	ei(num_gp) = e(k,j,i)
	946	dissi(num_gp) = diss(k,j,i)
	947	de_dxi(num_gp) = de_dx(k,j,i)
[1359]	948	de_dyi(num_gp) = 0.0_wp
[1929]	949	de_dzi(num_gp) = de_dz(k,j,i)
[849]	950	ENDIF
	951
	952	IF ( gp_outside_of_building(2) == 1 .AND. &
	953	gp_outside_of_building(1) == 0 ) THEN
	954	num_gp = num_gp + 1
	955	location(num_gp,1) = i * dx
[1359]	956	location(num_gp,2) = j * dy + 0.5_wp * dy
	957	location(num_gp,3) = k * dz - 0.5_wp * dz
[849]	958	ei(num_gp) = e(k,j+1,i)
	959	dissi(num_gp) = diss(k,j+1,i)
	960	de_dxi(num_gp) = de_dx(k,j+1,i)
[1359]	961	de_dyi(num_gp) = 0.0_wp
[849]	962	de_dzi(num_gp) = de_dz(k,j+1,i)
	963	ENDIF
	964
	965	!
	966	!-- If wall between gridpoint 3 and gridpoint 7, then
	967	!-- Neumann boundary condition has to be applied
	968	IF ( gp_outside_of_building(3) == 1 .AND. &
	969	gp_outside_of_building(7) == 0 ) THEN
	970	num_gp = num_gp + 1
[1359]	971	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	972	location(num_gp,2) = j * dy
[1359]	973	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	974	ei(num_gp) = e(k+1,j,i)
	975	dissi(num_gp) = diss(k+1,j,i)
[1359]	976	de_dxi(num_gp) = 0.0_wp
[849]	977	de_dyi(num_gp) = de_dy(k+1,j,i)
	978	de_dzi(num_gp) = de_dz(k+1,j,i)
	979	ENDIF
	980
	981	IF ( gp_outside_of_building(7) == 1 .AND. &
	982	gp_outside_of_building(3) == 0 ) THEN
	983	num_gp = num_gp + 1
[1359]	984	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	985	location(num_gp,2) = j * dy
[1359]	986	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	987	ei(num_gp) = e(k+1,j,i+1)
	988	dissi(num_gp) = diss(k+1,j,i+1)
[1359]	989	de_dxi(num_gp) = 0.0_wp
[849]	990	de_dyi(num_gp) = de_dy(k+1,j,i+1)
	991	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	992	ENDIF
	993
	994	!
	995	!-- If wall between gridpoint 7 and gridpoint 8, then
	996	!-- Neumann boundary condition has to be applied
	997	IF ( gp_outside_of_building(7) == 1 .AND. &
	998	gp_outside_of_building(8) == 0 ) THEN
	999	num_gp = num_gp + 1
	1000	location(num_gp,1) = (i+1) * dx
[1359]	1001	location(num_gp,2) = j * dy + 0.5_wp * dy
	1002	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	1003	ei(num_gp) = e(k+1,j,i+1)
	1004	dissi(num_gp) = diss(k+1,j,i+1)
	1005	de_dxi(num_gp) = de_dx(k+1,j,i+1)
[1359]	1006	de_dyi(num_gp) = 0.0_wp
[849]	1007	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	1008	ENDIF
	1009
	1010	IF ( gp_outside_of_building(8) == 1 .AND. &
	1011	gp_outside_of_building(7) == 0 ) THEN
	1012	num_gp = num_gp + 1
	1013	location(num_gp,1) = (i+1) * dx
[1359]	1014	location(num_gp,2) = j * dy + 0.5_wp * dy
	1015	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	1016	ei(num_gp) = e(k+1,j+1,i+1)
	1017	dissi(num_gp) = diss(k+1,j+1,i+1)
	1018	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
[1359]	1019	de_dyi(num_gp) = 0.0_wp
[849]	1020	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	1021	ENDIF
	1022
	1023	!
	1024	!-- If wall between gridpoint 4 and gridpoint 8, then
	1025	!-- Neumann boundary condition has to be applied
	1026	IF ( gp_outside_of_building(4) == 1 .AND. &
	1027	gp_outside_of_building(8) == 0 ) THEN
	1028	num_gp = num_gp + 1
[1359]	1029	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	1030	location(num_gp,2) = (j+1) * dy
[1359]	1031	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	1032	ei(num_gp) = e(k+1,j+1,i)
	1033	dissi(num_gp) = diss(k+1,j+1,i)
[1359]	1034	de_dxi(num_gp) = 0.0_wp
[849]	1035	de_dyi(num_gp) = de_dy(k+1,j+1,i)
	1036	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	1037	ENDIF
	1038
	1039	IF ( gp_outside_of_building(8) == 1 .AND. &
	1040	gp_outside_of_building(4) == 0 ) THEN
	1041	num_gp = num_gp + 1
[1359]	1042	location(num_gp,1) = i * dx + 0.5_wp * dx
[849]	1043	location(num_gp,2) = (j+1) * dy
[1359]	1044	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	1045	ei(num_gp) = e(k+1,j+1,i+1)
	1046	dissi(num_gp) = diss(k+1,j+1,i+1)
[1359]	1047	de_dxi(num_gp) = 0.0_wp
[849]	1048	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
	1049	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	1050	ENDIF
	1051
	1052	!
	1053	!-- If wall between gridpoint 3 and gridpoint 4, then
	1054	!-- Neumann boundary condition has to be applied
	1055	IF ( gp_outside_of_building(3) == 1 .AND. &
	1056	gp_outside_of_building(4) == 0 ) THEN
	1057	num_gp = num_gp + 1
	1058	location(num_gp,1) = i * dx
[1359]	1059	location(num_gp,2) = j * dy + 0.5_wp * dy
	1060	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	1061	ei(num_gp) = e(k+1,j,i)
	1062	dissi(num_gp) = diss(k+1,j,i)
	1063	de_dxi(num_gp) = de_dx(k+1,j,i)
[1359]	1064	de_dyi(num_gp) = 0.0_wp
[849]	1065	de_dzi(num_gp) = de_dz(k+1,j,i)
	1066	ENDIF
	1067
	1068	IF ( gp_outside_of_building(4) == 1 .AND. &
	1069	gp_outside_of_building(3) == 0 ) THEN
	1070	num_gp = num_gp + 1
	1071	location(num_gp,1) = i * dx
[1359]	1072	location(num_gp,2) = j * dy + 0.5_wp * dy
	1073	location(num_gp,3) = k * dz + 0.5_wp * dz
[849]	1074	ei(num_gp) = e(k+1,j+1,i)
	1075	dissi(num_gp) = diss(k+1,j+1,i)
	1076	de_dxi(num_gp) = de_dx(k+1,j+1,i)
[1359]	1077	de_dyi(num_gp) = 0.0_wp
[849]	1078	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	1079	ENDIF
	1080
	1081	!
	1082	!-- If wall between gridpoint 1 and gridpoint 3, then
	1083	!-- Neumann boundary condition has to be applied
	1084	!-- (only one case as only building beneath is possible)
	1085	IF ( gp_outside_of_building(1) == 0 .AND. &
	1086	gp_outside_of_building(3) == 1 ) THEN
	1087	num_gp = num_gp + 1
	1088	location(num_gp,1) = i * dx
	1089	location(num_gp,2) = j * dy
	1090	location(num_gp,3) = k * dz
	1091	ei(num_gp) = e(k+1,j,i)
	1092	dissi(num_gp) = diss(k+1,j,i)
	1093	de_dxi(num_gp) = de_dx(k+1,j,i)
	1094	de_dyi(num_gp) = de_dy(k+1,j,i)
[1359]	1095	de_dzi(num_gp) = 0.0_wp
[849]	1096	ENDIF
	1097
	1098	!
	1099	!-- If wall between gridpoint 5 and gridpoint 7, then
	1100	!-- Neumann boundary condition has to be applied
	1101	!-- (only one case as only building beneath is possible)
	1102	IF ( gp_outside_of_building(5) == 0 .AND. &
	1103	gp_outside_of_building(7) == 1 ) THEN
	1104	num_gp = num_gp + 1
	1105	location(num_gp,1) = (i+1) * dx
	1106	location(num_gp,2) = j * dy
	1107	location(num_gp,3) = k * dz
	1108	ei(num_gp) = e(k+1,j,i+1)
	1109	dissi(num_gp) = diss(k+1,j,i+1)
	1110	de_dxi(num_gp) = de_dx(k+1,j,i+1)
	1111	de_dyi(num_gp) = de_dy(k+1,j,i+1)
[1359]	1112	de_dzi(num_gp) = 0.0_wp
[849]	1113	ENDIF
	1114
	1115	!
	1116	!-- If wall between gridpoint 2 and gridpoint 4, then
	1117	!-- Neumann boundary condition has to be applied
	1118	!-- (only one case as only building beneath is possible)
	1119	IF ( gp_outside_of_building(2) == 0 .AND. &
	1120	gp_outside_of_building(4) == 1 ) THEN
	1121	num_gp = num_gp + 1
	1122	location(num_gp,1) = i * dx
	1123	location(num_gp,2) = (j+1) * dy
	1124	location(num_gp,3) = k * dz
	1125	ei(num_gp) = e(k+1,j+1,i)
	1126	dissi(num_gp) = diss(k+1,j+1,i)
	1127	de_dxi(num_gp) = de_dx(k+1,j+1,i)
	1128	de_dyi(num_gp) = de_dy(k+1,j+1,i)
[1359]	1129	de_dzi(num_gp) = 0.0_wp
[849]	1130	ENDIF
	1131
[1929]	1132	!
[849]	1133	!-- If wall between gridpoint 6 and gridpoint 8, then
	1134	!-- Neumann boundary condition has to be applied
	1135	!-- (only one case as only building beneath is possible)
	1136	IF ( gp_outside_of_building(6) == 0 .AND. &
	1137	gp_outside_of_building(8) == 1 ) THEN
	1138	num_gp = num_gp + 1
	1139	location(num_gp,1) = (i+1) * dx
	1140	location(num_gp,2) = (j+1) * dy
	1141	location(num_gp,3) = k * dz
	1142	ei(num_gp) = e(k+1,j+1,i+1)
	1143	dissi(num_gp) = diss(k+1,j+1,i+1)
	1144	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
	1145	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
[1359]	1146	de_dzi(num_gp) = 0.0_wp
[849]	1147	ENDIF
[1929]	1148
[849]	1149	!
	1150	!-- Carry out the interpolation
	1151	IF ( num_gp == 1 ) THEN
[1929]	1152	!
[849]	1153	!-- If only one of the gridpoints is situated outside of the
	1154	!-- building, it follows that the values at the particle
	1155	!-- location are the same as the gridpoint values
[1929]	1156	e_int(n) = ei(num_gp)
	1157	diss_int(n) = dissi(num_gp)
[1359]	1158	de_dx_int(n) = de_dxi(num_gp)
	1159	de_dy_int(n) = de_dyi(num_gp)
	1160	de_dz_int(n) = de_dzi(num_gp)
[1929]	1161	!
	1162	!-- Set flag for stochastic equation which disables calculation
	1163	!-- of drift and memory term near topography.
	1164	term_1_2(n) = 0.0_wp
[849]	1165	ELSE IF ( num_gp > 1 ) THEN
[1929]	1166
[1359]	1167	d_sum = 0.0_wp
[1929]	1168	!
[849]	1169	!-- Evaluation of the distances between the gridpoints
	1170	!-- contributing to the interpolated values, and the particle
	1171	!-- location
	1172	DO agp = 1, num_gp
	1173	d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 &
	1174	+ ( particles(n)%y-location(agp,2) )**2 &
[1359]	1175	+ ( zv(n)-location(agp,3) )**2
[849]	1176	d_sum = d_sum + d_gp_pl(agp)
	1177	ENDDO
[1929]	1178
[849]	1179	!
	1180	!-- Finally the interpolation can be carried out
[1359]	1181	e_int(n) = 0.0_wp
	1182	diss_int(n) = 0.0_wp
[1929]	1183	de_dx_int(n) = 0.0_wp
	1184	de_dy_int(n) = 0.0_wp
	1185	de_dz_int(n) = 0.0_wp
[849]	1186	DO agp = 1, num_gp
[1359]	1187	e_int(n) = e_int(n) + ( d_sum - d_gp_pl(agp) ) * &
[849]	1188	ei(agp) / ( (num_gp-1) * d_sum )
[1359]	1189	diss_int(n) = diss_int(n) + ( d_sum - d_gp_pl(agp) ) * &
[849]	1190	dissi(agp) / ( (num_gp-1) * d_sum )
[1359]	1191	de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
[849]	1192	de_dxi(agp) / ( (num_gp-1) * d_sum )
[1359]	1193	de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
[849]	1194	de_dyi(agp) / ( (num_gp-1) * d_sum )
[1359]	1195	de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
[849]	1196	de_dzi(agp) / ( (num_gp-1) * d_sum )
	1197	ENDDO
[1929]	1198
[849]	1199	ENDIF
[1929]	1200	e_int(n) = MAX( 1E-20_wp, e_int(n) )
	1201	diss_int(n) = MAX( 1E-20_wp, diss_int(n) )
	1202	de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) )
	1203	de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) )
	1204	de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) )
	1205	!
	1206	!-- Set flag for stochastic equation which disables calculation
	1207	!-- of drift and memory term near topography.
	1208	term_1_2(n) = 0.0_wp
[849]	1209	ENDIF
[1359]	1210	ENDDO
	1211	ENDIF
[849]	1212
[1359]	1213	DO nb = 0,7
	1214	i = ip + block_offset(nb)%i_off
	1215	j = jp + block_offset(nb)%j_off
	1216	k = kp + block_offset(nb)%k_off
[849]	1217
[1359]	1218	DO n = start_index(nb), end_index(nb)
[849]	1219	!
[1359]	1220	!-- Vertical interpolation of the horizontally averaged SGS TKE and
	1221	!-- resolved-scale velocity variances and use the interpolated values
	1222	!-- to calculate the coefficient fs, which is a measure of the ratio
	1223	!-- of the subgrid-scale turbulent kinetic energy to the total amount
	1224	!-- of turbulent kinetic energy.
	1225	IF ( k == 0 ) THEN
	1226	e_mean_int = hom(0,1,8,0)
	1227	ELSE
	1228	e_mean_int = hom(k,1,8,0) + &
	1229	( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
	1230	( zu(k+1) - zu(k) ) * &
	1231	( zv(n) - zu(k) )
	1232	ENDIF
[849]	1233
[1685]	1234	kw = kp - 1
[849]	1235
[1359]	1236	IF ( k == 0 ) THEN
	1237	aa = hom(k+1,1,30,0) * ( zv(n) / &
	1238	( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
	1239	bb = hom(k+1,1,31,0) * ( zv(n) / &
	1240	( 0.5_wp * ( zu(k+1) - zu(k) ) ) )
	1241	cc = hom(kw+1,1,32,0) * ( zv(n) / &
	1242	( 1.0_wp * ( zw(kw+1) - zw(kw) ) ) )
	1243	ELSE
	1244	aa = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
	1245	( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
	1246	bb = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
	1247	( ( zv(n) - zu(k) ) / ( zu(k+1) - zu(k) ) )
	1248	cc = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) * &
	1249	( ( zv(n) - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
	1250	ENDIF
[849]	1251
[1359]	1252	vv_int = ( 1.0_wp / 3.0_wp ) * ( aa + bb + cc )
	1253	!
	1254	!-- Needed to avoid NaN particle velocities. The value of 1.0 is just
	1255	!-- an educated guess for the given case.
	1256	IF ( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int == 0.0_wp ) THEN
	1257	fs_int(n) = 1.0_wp
	1258	ELSE
	1259	fs_int(n) = ( 2.0_wp / 3.0_wp ) * e_mean_int / &
	1260	( vv_int + ( 2.0_wp / 3.0_wp ) * e_mean_int )
	1261	ENDIF
[849]	1262
[1359]	1263	ENDDO
	1264	ENDDO
[849]	1265
[1359]	1266	DO n = 1, number_of_particles
	1267
	1268	rg(n,1) = random_gauss( iran_part, 5.0_wp )
	1269	rg(n,2) = random_gauss( iran_part, 5.0_wp )
	1270	rg(n,3) = random_gauss( iran_part, 5.0_wp )
	1271
	1272	ENDDO
	1273
	1274	DO n = 1, number_of_particles
[849]	1275	!
	1276	!-- Calculate the Lagrangian timescale according to Weil et al. (2004).
[1929]	1277	lagr_timescale = ( 4.0_wp * e_int(n) + 1E-20_wp ) / &
	1278	( 3.0_wp * fs_int(n) * c_0 * diss_int(n) + 1E-20_wp )
[849]	1279
	1280	!
	1281	!-- Calculate the next particle timestep. dt_gap is the time needed to
	1282	!-- complete the current LES timestep.
	1283	dt_gap = dt_3d - particles(n)%dt_sum
[1359]	1284	dt_particle(n) = MIN( dt_3d, 0.025_wp * lagr_timescale, dt_gap )
[849]	1285
	1286	!
	1287	!-- The particle timestep should not be too small in order to prevent
	1288	!-- the number of particle timesteps of getting too large
[1359]	1289	IF ( dt_particle(n) < dt_min_part .AND. dt_min_part < dt_gap ) THEN
	1290	dt_particle(n) = dt_min_part
[849]	1291	ENDIF
	1292
	1293	!
	1294	!-- Calculate the SGS velocity components
[1359]	1295	IF ( particles(n)%age == 0.0_wp ) THEN
[849]	1296	!
	1297	!-- For new particles the SGS components are derived from the SGS
	1298	!-- TKE. Limit the Gaussian random number to the interval
	1299	!-- [-5.0sigma, 5.0sigma] in order to prevent the SGS velocities
	1300	!-- from becoming unrealistically large.
[1929]	1301	particles(n)%rvar1 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * &
[1359]	1302	( rg(n,1) - 1.0_wp )
[1929]	1303	particles(n)%rvar2 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * &
[1359]	1304	( rg(n,2) - 1.0_wp )
[1929]	1305	particles(n)%rvar3 = SQRT( 2.0_wp * sgs_wf_part * e_int(n) + 1E-20_wp ) * &
[1359]	1306	( rg(n,3) - 1.0_wp )
[849]	1307
	1308	ELSE
	1309	!
	1310	!-- Restriction of the size of the new timestep: compared to the
	1311	!-- previous timestep the increase must not exceed 200%
	1312
	1313	dt_particle_m = particles(n)%age - particles(n)%age_m
[1359]	1314	IF ( dt_particle(n) > 2.0_wp * dt_particle_m ) THEN
	1315	dt_particle(n) = 2.0_wp * dt_particle_m
[849]	1316	ENDIF
	1317
	1318	!
	1319	!-- For old particles the SGS components are correlated with the
	1320	!-- values from the previous timestep. Random numbers have also to
	1321	!-- be limited (see above).
	1322	!-- As negative values for the subgrid TKE are not allowed, the
	1323	!-- change of the subgrid TKE with time cannot be smaller than
[1359]	1324	!-- -e_int(n)/dt_particle. This value is used as a lower boundary
[849]	1325	!-- value for the change of TKE
	1326
[1359]	1327	de_dt_min = - e_int(n) / dt_particle(n)
[849]	1328
[1359]	1329	de_dt = ( e_int(n) - particles(n)%e_m ) / dt_particle_m
[849]	1330
	1331	IF ( de_dt < de_dt_min ) THEN
	1332	de_dt = de_dt_min
	1333	ENDIF
	1334
[1929]	1335	CALL weil_stochastic_eq(particles(n)%rvar1, fs_int(n), e_int(n), &
	1336	de_dx_int(n), de_dt, diss_int(n), &
	1337	dt_particle(n), rg(n,1), term_1_2(n) )
[849]	1338
[1929]	1339	CALL weil_stochastic_eq(particles(n)%rvar2, fs_int(n), e_int(n), &
	1340	de_dy_int(n), de_dt, diss_int(n), &
	1341	dt_particle(n), rg(n,2), term_1_2(n) )
[849]	1342
[1929]	1343	CALL weil_stochastic_eq(particles(n)%rvar3, fs_int(n), e_int(n), &
	1344	de_dz_int(n), de_dt, diss_int(n), &
	1345	dt_particle(n), rg(n,3), term_1_2(n) )
[849]	1346
	1347	ENDIF
[1929]	1348
[1359]	1349	u_int(n) = u_int(n) + particles(n)%rvar1
	1350	v_int(n) = v_int(n) + particles(n)%rvar2
	1351	w_int(n) = w_int(n) + particles(n)%rvar3
[849]	1352	!
	1353	!-- Store the SGS TKE of the current timelevel which is needed for
	1354	!-- for calculating the SGS particle velocities at the next timestep
[1359]	1355	particles(n)%e_m = e_int(n)
	1356	ENDDO
[849]	1357
[1359]	1358	ELSE
[849]	1359	!
[1359]	1360	!-- If no SGS velocities are used, only the particle timestep has to
	1361	!-- be set
	1362	dt_particle = dt_3d
[849]	1363
[1359]	1364	ENDIF
[849]	1365	!
[1359]	1366	!-- Store the old age of the particle ( needed to prevent that a
	1367	!-- particle crosses several PEs during one timestep, and for the
	1368	!-- evaluation of the subgrid particle velocity fluctuations )
	1369	particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
[849]	1370
[1359]	1371	dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
[849]	1372
[1359]	1373	IF ( ANY( dens_ratio == 0.0_wp ) ) THEN
	1374	DO n = 1, number_of_particles
	1375
[849]	1376	!
[1359]	1377	!-- Particle advection
	1378	IF ( dens_ratio(n) == 0.0_wp ) THEN
[849]	1379	!
[1359]	1380	!-- Pure passive transport (without particle inertia)
	1381	particles(n)%x = xv(n) + u_int(n) * dt_particle(n)
	1382	particles(n)%y = yv(n) + v_int(n) * dt_particle(n)
	1383	particles(n)%z = zv(n) + w_int(n) * dt_particle(n)
[849]	1384
[1359]	1385	particles(n)%speed_x = u_int(n)
	1386	particles(n)%speed_y = v_int(n)
	1387	particles(n)%speed_z = w_int(n)
[849]	1388
[1359]	1389	ELSE
[849]	1390	!
[1359]	1391	!-- Transport of particles with inertia
	1392	particles(n)%x = particles(n)%x + particles(n)%speed_x * &
	1393	dt_particle(n)
	1394	particles(n)%y = particles(n)%y + particles(n)%speed_y * &
	1395	dt_particle(n)
	1396	particles(n)%z = particles(n)%z + particles(n)%speed_z * &
	1397	dt_particle(n)
[849]	1398
	1399	!
[1359]	1400	!-- Update of the particle velocity
	1401	IF ( cloud_droplets ) THEN
[1822]	1402	!
	1403	!-- Terminal velocity is computed for vertical direction (Rogers et
	1404	!-- al., 1993, J. Appl. Meteorol.)
	1405	diameter = particles(n)%radius * 2000.0_wp !diameter in mm
	1406	IF ( diameter <= d0_rog ) THEN
	1407	w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
	1408	ELSE
	1409	w_s = a_rog - b_rog * EXP( -c_rog * diameter )
	1410	ENDIF
[1359]	1411
[1822]	1412	!
	1413	!-- If selected, add random velocities following Soelch and Kaercher
	1414	!-- (2010, Q. J. R. Meteorol. Soc.)
	1415	IF ( use_sgs_for_particles ) THEN
	1416	lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
	1417	RL = EXP( -1.0_wp * dt_3d / lagr_timescale )
	1418	sigma = SQRT( e(kp,jp,ip) )
	1419
	1420	rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
	1421	rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
	1422	rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
	1423
	1424	particles(n)%rvar1 = RL * particles(n)%rvar1 + &
	1425	SQRT( 1.0_wp - RL*2 ) sigma * rg1
	1426	particles(n)%rvar2 = RL * particles(n)%rvar2 + &
	1427	SQRT( 1.0_wp - RL*2 ) sigma * rg2
	1428	particles(n)%rvar3 = RL * particles(n)%rvar3 + &
	1429	SQRT( 1.0_wp - RL*2 ) sigma * rg3
	1430
	1431	particles(n)%speed_x = u_int(n) + particles(n)%rvar1
	1432	particles(n)%speed_y = v_int(n) + particles(n)%rvar2
	1433	particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
	1434	ELSE
	1435	particles(n)%speed_x = u_int(n)
	1436	particles(n)%speed_y = v_int(n)
	1437	particles(n)%speed_z = w_int(n) - w_s
	1438	ENDIF
	1439
[1359]	1440	ELSE
[1822]	1441
	1442	IF ( use_sgs_for_particles ) THEN
	1443	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1444	exp_term = EXP( -exp_arg * dt_particle(n) )
	1445	ELSE
	1446	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1447	exp_term = particle_groups(particles(n)%group)%exp_term
	1448	ENDIF
	1449	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
	1450	u_int(n) * ( 1.0_wp - exp_term )
	1451	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
	1452	v_int(n) * ( 1.0_wp - exp_term )
	1453	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
	1454	( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * &
	1455	g / exp_arg ) * ( 1.0_wp - exp_term )
[1359]	1456	ENDIF
[1822]	1457
[1359]	1458	ENDIF
	1459
	1460	ENDDO
	1461
	1462	ELSE
	1463
	1464	DO n = 1, number_of_particles
	1465
	1466	!-- Transport of particles with inertia
	1467	particles(n)%x = xv(n) + particles(n)%speed_x * dt_particle(n)
	1468	particles(n)%y = yv(n) + particles(n)%speed_y * dt_particle(n)
	1469	particles(n)%z = zv(n) + particles(n)%speed_z * dt_particle(n)
	1470	!
[849]	1471	!-- Update of the particle velocity
	1472	IF ( cloud_droplets ) THEN
[1822]	1473	!
	1474	!-- Terminal velocity is computed for vertical direction (Rogers et al.,
	1475	!-- 1993, J. Appl. Meteorol.)
	1476	diameter = particles(n)%radius * 2000.0_wp !diameter in mm
	1477	IF ( diameter <= d0_rog ) THEN
	1478	w_s = k_cap_rog * diameter * ( 1.0_wp - EXP( -k_low_rog * diameter ) )
	1479	ELSE
	1480	w_s = a_rog - b_rog * EXP( -c_rog * diameter )
	1481	ENDIF
[1359]	1482
[1822]	1483	!
	1484	!-- If selected, add random velocities following Soelch and Kaercher
	1485	!-- (2010, Q. J. R. Meteorol. Soc.)
	1486	IF ( use_sgs_for_particles ) THEN
	1487	lagr_timescale = km(kp,jp,ip) / MAX( e(kp,jp,ip), 1.0E-20_wp )
	1488	RL = EXP( -1.0_wp * dt_3d / lagr_timescale )
	1489	sigma = SQRT( e(kp,jp,ip) )
[1359]	1490
[1822]	1491	rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
	1492	rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
	1493	rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
	1494
	1495	particles(n)%rvar1 = RL * particles(n)%rvar1 + &
	1496	SQRT( 1.0_wp - RL*2 ) sigma * rg1
	1497	particles(n)%rvar2 = RL * particles(n)%rvar2 + &
	1498	SQRT( 1.0_wp - RL*2 ) sigma * rg2
	1499	particles(n)%rvar3 = RL * particles(n)%rvar3 + &
	1500	SQRT( 1.0_wp - RL*2 ) sigma * rg3
	1501
	1502	particles(n)%speed_x = u_int(n) + particles(n)%rvar1
	1503	particles(n)%speed_y = v_int(n) + particles(n)%rvar2
	1504	particles(n)%speed_z = w_int(n) + particles(n)%rvar3 - w_s
	1505	ELSE
	1506	particles(n)%speed_x = u_int(n)
	1507	particles(n)%speed_y = v_int(n)
	1508	particles(n)%speed_z = w_int(n) - w_s
	1509	ENDIF
	1510
[849]	1511	ELSE
[1822]	1512
	1513	IF ( use_sgs_for_particles ) THEN
	1514	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1515	exp_term = EXP( -exp_arg * dt_particle(n) )
	1516	ELSE
	1517	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1518	exp_term = particle_groups(particles(n)%group)%exp_term
	1519	ENDIF
	1520	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
	1521	u_int(n) * ( 1.0_wp - exp_term )
	1522	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
	1523	v_int(n) * ( 1.0_wp - exp_term )
	1524	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
	1525	( w_int(n) - ( 1.0_wp - dens_ratio(n) ) * g / &
	1526	exp_arg ) * ( 1.0_wp - exp_term )
[849]	1527	ENDIF
[1822]	1528
[1359]	1529	ENDDO
[849]	1530
[1359]	1531	ENDIF
	1532
	1533	DO n = 1, number_of_particles
[849]	1534	!
	1535	!-- Increment the particle age and the total time that the particle
	1536	!-- has advanced within the particle timestep procedure
[1359]	1537	particles(n)%age = particles(n)%age + dt_particle(n)
	1538	particles(n)%dt_sum = particles(n)%dt_sum + dt_particle(n)
[849]	1539
	1540	!
	1541	!-- Check whether there is still a particle that has not yet completed
	1542	!-- the total LES timestep
[1359]	1543	IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp ) THEN
[849]	1544	dt_3d_reached_l = .FALSE.
	1545	ENDIF
	1546
	1547	ENDDO
	1548
[1359]	1549	CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
[849]	1550
[1929]	1551
[849]	1552	END SUBROUTINE lpm_advec
[1929]	1553
	1554	! Description:
	1555	! ------------
	1556	!> Calculation of subgrid-scale particle speed using the stochastic model
	1557	!> of Weil et al. (2004, JAS, 61, 2877-2887).
	1558	!------------------------------------------------------------------------------!
	1559	SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n, &
	1560	dt_n, rg_n, fac )
	1561
	1562	USE kinds
	1563
	1564	USE particle_attributes, &
	1565	ONLY: c_0, sgs_wf_part
	1566
	1567	IMPLICIT NONE
	1568
	1569	REAL(wp) :: a1 !< dummy argument
	1570	REAL(wp) :: dedt_n !< time derivative of TKE at particle position
	1571	REAL(wp) :: dedxi_n !< horizontal derivative of TKE at particle position
	1572	REAL(wp) :: diss_n !< dissipation at particle position
	1573	REAL(wp) :: dt_n !< particle timestep
	1574	REAL(wp) :: e_n !< TKE at particle position
	1575	REAL(wp) :: fac !< flag to identify adjacent topography
	1576	REAL(wp) :: fs_n !< weighting factor to prevent that subgrid-scale particle speed becomes too large
	1577	REAL(wp) :: sgs_w !< constant (1/3)
	1578	REAL(wp) :: rg_n !< random number
	1579	REAL(wp) :: term1 !< memory term
	1580	REAL(wp) :: term2 !< drift correction term
	1581	REAL(wp) :: term3 !< random term
	1582	REAL(wp) :: v_sgs !< subgrid-scale velocity component
	1583
[2100]	1584	!-- At first, limit TKE to a small non-zero number, in order to prevent
	1585	!-- the occurrence of extremely large SGS-velocities in case TKE is zero,
	1586	!-- (could occur at the simulation begin).
	1587	e_n = MAX( e_n, 1E-20_wp )
[1929]	1588	!
	1589	!-- Please note, terms 1 and 2 (drift and memory term, respectively) are
	1590	!-- multiplied by a flag to switch of both terms near topography.
	1591	!-- This is necessary, as both terms may cause a subgrid-scale velocity build up
	1592	!-- if particles are trapped in regions with very small TKE, e.g. in narrow street
	1593	!-- canyons resolved by only a few grid points. Hence, term 1 and term 2 are
	1594	!-- disabled if one of the adjacent grid points belongs to topography.
	1595	!-- Moreover, in this case, the previous subgrid-scale component is also set
	1596	!-- to zero.
	1597
	1598	a1 = fs_n * c_0 * diss_n
	1599	!
	1600	!-- Memory term
	1601	term1 = - a1 * v_sgs * dt_n / ( 4.0_wp * sgs_wf_part * e_n + 1E-20_wp ) &
	1602	* fac
	1603	!
	1604	!-- Drift correction term
	1605	term2 = ( ( dedt_n * v_sgs / e_n ) + dedxi_n ) * 0.5_wp * dt_n &
	1606	* fac
	1607	!
	1608	!-- Random term
	1609	term3 = SQRT( MAX( a1, 1E-20 ) ) * ( rg_n - 1.0_wp ) * SQRT( dt_n )
	1610	!
	1611	!-- In cese one of the adjacent grid-boxes belongs to topograhy, the previous
	1612	!-- subgrid-scale velocity component is set to zero, in order to prevent a
	1613	!-- velocity build-up.
	1614
	1615	!-- This case, set also previous subgrid-scale component to zero.
	1616	v_sgs = v_sgs * fac + term1 + term2 + term3
	1617
	1618	END SUBROUTINE weil_stochastic_eq

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |