Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |