Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |