Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |