Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |