Home

Context Navigation

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

Last change on this file since 1217 was 1037, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 47.7 KB

Rev	Line
[849]	1	SUBROUTINE lpm_advec
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
	17	! Copyright 1997-2012 Leibniz University Hannover
	18	!--------------------------------------------------------------------------------!
	19	!
[849]	20	! Current revisions:
	21	! ------------------
	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: lpm_advec.f90 1037 2012-10-22 14:10:22Z raasch $
	27	!
[1037]	28	! 1036 2012-10-22 13:43:42Z raasch
	29	! code put under GPL (PALM 3.9)
	30	!
[850]	31	! 849 2012-03-15 10:35:09Z raasch
	32	! initial revision (former part of advec_particles)
[849]	33	!
[850]	34	!
[849]	35	! Description:
	36	! ------------
	37	! Calculation of new particle positions due to advection using a simple Euler
	38	! scheme. Particles may feel inertia effects. SGS transport can be included
	39	! using the stochastic model of Weil et al. (2004, JAS, 61, 2877-2887).
	40	!------------------------------------------------------------------------------!
	41
	42	USE arrays_3d
	43	USE control_parameters
	44	USE grid_variables
	45	USE indices
	46	USE particle_attributes
	47	USE statistics
	48
	49	IMPLICIT NONE
	50
	51	INTEGER :: i, j, k, n
	52
	53	REAL :: aa, bb, cc, dd, dens_ratio, exp_arg, exp_term, gg, u_int, &
	54	u_int_l, u_int_u, v_int, v_int_l, v_int_u, w_int, w_int_l, &
	55	w_int_u, x, y
	56
	57
	58	INTEGER :: agp, kw, num_gp
	59	INTEGER :: gp_outside_of_building(1:8)
	60
	61	REAL :: d_sum, de_dx_int, de_dx_int_l, de_dx_int_u, de_dy_int, &
	62	de_dy_int_l, de_dy_int_u, de_dt, de_dt_min, de_dz_int, &
	63	de_dz_int_l, de_dz_int_u, diss_int, diss_int_l, diss_int_u, &
	64	dt_gap, dt_particle, dt_particle_m, e_int, e_int_l, e_int_u, &
	65	e_mean_int, fs_int, lagr_timescale, random_gauss, vv_int
	66
	67	REAL :: location(1:30,1:3)
	68	REAL, DIMENSION(1:30) :: de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
	69
	70
	71	DO n = 1, number_of_particles
	72
	73	!
	74	!-- Move particle only if the LES timestep has not (approximately) been
	75	!-- reached
	76	IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 ) CYCLE
	77
	78	!
	79	!-- Interpolate u velocity-component, determine left, front, bottom
	80	!-- index of u-array
	81	i = ( particles(n)%x + 0.5 * dx ) * ddx
	82	j = particles(n)%y * ddy
	83	k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz &
	84	+ offset_ocean_nzt ! only exact if equidistant
	85
	86	!
	87	!-- Interpolation of the velocity components in the xy-plane
	88	x = particles(n)%x + ( 0.5 - i ) * dx
	89	y = particles(n)%y - j * dy
	90	aa = x2 + y2
	91	bb = ( dx - x )2 + y2
	92	cc = x2 + ( dy - y )2
	93	dd = ( dx - x )2 + ( dy - y )2
	94	gg = aa + bb + cc + dd
	95
	96	u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) &
	97	+ ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
	98	) / ( 3.0 * gg ) - u_gtrans
	99	IF ( k+1 == nzt+1 ) THEN
	100	u_int = u_int_l
	101	ELSE
	102	u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) &
	103	+ ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
	104	) / ( 3.0 * gg ) - u_gtrans
	105	u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
	106	( u_int_u - u_int_l )
	107	ENDIF
	108
	109	!
	110	!-- Same procedure for interpolation of the v velocity-component (adopt
	111	!-- index k from u velocity-component)
	112	i = particles(n)%x * ddx
	113	j = ( particles(n)%y + 0.5 * dy ) * ddy
	114
	115	x = particles(n)%x - i * dx
	116	y = particles(n)%y + ( 0.5 - j ) * dy
	117	aa = x2 + y2
	118	bb = ( dx - x )2 + y2
	119	cc = x2 + ( dy - y )2
	120	dd = ( dx - x )2 + ( dy - y )2
	121	gg = aa + bb + cc + dd
	122
	123	v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) &
	124	+ ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
	125	) / ( 3.0 * gg ) - v_gtrans
	126	IF ( k+1 == nzt+1 ) THEN
	127	v_int = v_int_l
	128	ELSE
	129	v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) &
	130	+ ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
	131	) / ( 3.0 * gg ) - v_gtrans
	132	v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
	133	( v_int_u - v_int_l )
	134	ENDIF
	135
	136	!
	137	!-- Same procedure for interpolation of the w velocity-component (adopt
	138	!-- index i from v velocity-component)
	139	IF ( vertical_particle_advection(particles(n)%group) ) THEN
	140	j = particles(n)%y * ddy
	141	k = particles(n)%z / dz + offset_ocean_nzt_m1
	142
	143	x = particles(n)%x - i * dx
	144	y = particles(n)%y - j * dy
	145	aa = x2 + y2
	146	bb = ( dx - x )2 + y2
	147	cc = x2 + ( dy - y )2
	148	dd = ( dx - x )2 + ( dy - y )2
	149	gg = aa + bb + cc + dd
	150
	151	w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) &
	152	+ ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
	153	) / ( 3.0 * gg )
	154	IF ( k+1 == nzt+1 ) THEN
	155	w_int = w_int_l
	156	ELSE
	157	w_int_u = ( ( gg-aa ) * w(k+1,j,i) + &
	158	( gg-bb ) * w(k+1,j,i+1) + &
	159	( gg-cc ) * w(k+1,j+1,i) + &
	160	( gg-dd ) * w(k+1,j+1,i+1) &
	161	) / ( 3.0 * gg )
	162	w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
	163	( w_int_u - w_int_l )
	164	ENDIF
	165	ELSE
	166	w_int = 0.0
	167	ENDIF
	168
	169	!
	170	!-- Interpolate and calculate quantities needed for calculating the SGS
	171	!-- velocities
	172	IF ( use_sgs_for_particles ) THEN
	173	!
	174	!-- Interpolate TKE
	175	i = particles(n)%x * ddx
	176	j = particles(n)%y * ddy
	177	k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz &
	178	+ offset_ocean_nzt ! only exact if eq.dist
	179
	180	IF ( topography == 'flat' ) THEN
	181
	182	x = particles(n)%x - i * dx
	183	y = particles(n)%y - j * dy
	184	aa = x2 + y2
	185	bb = ( dx - x )2 + y2
	186	cc = x2 + ( dy - y )2
	187	dd = ( dx - x )2 + ( dy - y )2
	188	gg = aa + bb + cc + dd
	189
	190	e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
	191	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
	192	) / ( 3.0 * gg )
	193
	194	IF ( k+1 == nzt+1 ) THEN
	195	e_int = e_int_l
	196	ELSE
	197	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
	198	( gg - bb ) * e(k+1,j,i+1) + &
	199	( gg - cc ) * e(k+1,j+1,i) + &
	200	( gg - dd ) * e(k+1,j+1,i+1) &
	201	) / ( 3.0 * gg )
	202	e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
	203	( e_int_u - e_int_l )
	204	ENDIF
	205
	206	!
	207	!-- Interpolate the TKE gradient along x (adopt incides i,j,k and
	208	!-- all position variables from above (TKE))
	209	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
	210	( gg - bb ) * de_dx(k,j,i+1) + &
	211	( gg - cc ) * de_dx(k,j+1,i) + &
	212	( gg - dd ) * de_dx(k,j+1,i+1) &
	213	) / ( 3.0 * gg )
	214
	215	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	216	de_dx_int = de_dx_int_l
	217	ELSE
	218	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
	219	( gg - bb ) * de_dx(k+1,j,i+1) + &
	220	( gg - cc ) * de_dx(k+1,j+1,i) + &
	221	( gg - dd ) * de_dx(k+1,j+1,i+1) &
	222	) / ( 3.0 * gg )
	223	de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
	224	( de_dx_int_u - de_dx_int_l )
	225	ENDIF
	226
	227	!
	228	!-- Interpolate the TKE gradient along y
	229	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
	230	( gg - bb ) * de_dy(k,j,i+1) + &
	231	( gg - cc ) * de_dy(k,j+1,i) + &
	232	( gg - dd ) * de_dy(k,j+1,i+1) &
	233	) / ( 3.0 * gg )
	234	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	235	de_dy_int = de_dy_int_l
	236	ELSE
	237	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
	238	( gg - bb ) * de_dy(k+1,j,i+1) + &
	239	( gg - cc ) * de_dy(k+1,j+1,i) + &
	240	( gg - dd ) * de_dy(k+1,j+1,i+1) &
	241	) / ( 3.0 * gg )
	242	de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
	243	( de_dy_int_u - de_dy_int_l )
	244	ENDIF
	245
	246	!
	247	!-- Interpolate the TKE gradient along z
	248	IF ( particles(n)%z < 0.5 * dz ) THEN
	249	de_dz_int = 0.0
	250	ELSE
	251	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
	252	( gg - bb ) * de_dz(k,j,i+1) + &
	253	( gg - cc ) * de_dz(k,j+1,i) + &
	254	( gg - dd ) * de_dz(k,j+1,i+1) &
	255	) / ( 3.0 * gg )
	256
	257	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	258	de_dz_int = de_dz_int_l
	259	ELSE
	260	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
	261	( gg - bb ) * de_dz(k+1,j,i+1) + &
	262	( gg - cc ) * de_dz(k+1,j+1,i) + &
	263	( gg - dd ) * de_dz(k+1,j+1,i+1) &
	264	) / ( 3.0 * gg )
	265	de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
	266	( de_dz_int_u - de_dz_int_l )
	267	ENDIF
	268	ENDIF
	269
	270	!
	271	!-- Interpolate the dissipation of TKE
	272	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
	273	( gg - bb ) * diss(k,j,i+1) + &
	274	( gg - cc ) * diss(k,j+1,i) + &
	275	( gg - dd ) * diss(k,j+1,i+1) &
	276	) / ( 3.0 * gg )
	277
	278	IF ( k+1 == nzt+1 ) THEN
	279	diss_int = diss_int_l
	280	ELSE
	281	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
	282	( gg - bb ) * diss(k+1,j,i+1) + &
	283	( gg - cc ) * diss(k+1,j+1,i) + &
	284	( gg - dd ) * diss(k+1,j+1,i+1) &
	285	) / ( 3.0 * gg )
	286	diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
	287	( diss_int_u - diss_int_l )
	288	ENDIF
	289
	290	ELSE
	291
	292	!
	293	!-- In case that there are buildings it has to be determined
	294	!-- how many of the gridpoints defining the particle box are
	295	!-- situated within a building
	296	!-- gp_outside_of_building(1): i,j,k
	297	!-- gp_outside_of_building(2): i,j+1,k
	298	!-- gp_outside_of_building(3): i,j,k+1
	299	!-- gp_outside_of_building(4): i,j+1,k+1
	300	!-- gp_outside_of_building(5): i+1,j,k
	301	!-- gp_outside_of_building(6): i+1,j+1,k
	302	!-- gp_outside_of_building(7): i+1,j,k+1
	303	!-- gp_outside_of_building(8): i+1,j+1,k+1
	304
	305	gp_outside_of_building = 0
	306	location = 0.0
	307	num_gp = 0
	308
	309	IF ( k > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
	310	num_gp = num_gp + 1
	311	gp_outside_of_building(1) = 1
	312	location(num_gp,1) = i * dx
	313	location(num_gp,2) = j * dy
	314	location(num_gp,3) = k * dz - 0.5 * dz
	315	ei(num_gp) = e(k,j,i)
	316	dissi(num_gp) = diss(k,j,i)
	317	de_dxi(num_gp) = de_dx(k,j,i)
	318	de_dyi(num_gp) = de_dy(k,j,i)
	319	de_dzi(num_gp) = de_dz(k,j,i)
	320	ENDIF
	321
	322	IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) &
	323	THEN
	324	num_gp = num_gp + 1
	325	gp_outside_of_building(2) = 1
	326	location(num_gp,1) = i * dx
	327	location(num_gp,2) = (j+1) * dy
	328	location(num_gp,3) = k * dz - 0.5 * dz
	329	ei(num_gp) = e(k,j+1,i)
	330	dissi(num_gp) = diss(k,j+1,i)
	331	de_dxi(num_gp) = de_dx(k,j+1,i)
	332	de_dyi(num_gp) = de_dy(k,j+1,i)
	333	de_dzi(num_gp) = de_dz(k,j+1,i)
	334	ENDIF
	335
	336	IF ( k+1 > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
	337	num_gp = num_gp + 1
	338	gp_outside_of_building(3) = 1
	339	location(num_gp,1) = i * dx
	340	location(num_gp,2) = j * dy
	341	location(num_gp,3) = (k+1) * dz - 0.5 * dz
	342	ei(num_gp) = e(k+1,j,i)
	343	dissi(num_gp) = diss(k+1,j,i)
	344	de_dxi(num_gp) = de_dx(k+1,j,i)
	345	de_dyi(num_gp) = de_dy(k+1,j,i)
	346	de_dzi(num_gp) = de_dz(k+1,j,i)
	347	ENDIF
	348
	349	IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) &
	350	THEN
	351	num_gp = num_gp + 1
	352	gp_outside_of_building(4) = 1
	353	location(num_gp,1) = i * dx
	354	location(num_gp,2) = (j+1) * dy
	355	location(num_gp,3) = (k+1) * dz - 0.5 * dz
	356	ei(num_gp) = e(k+1,j+1,i)
	357	dissi(num_gp) = diss(k+1,j+1,i)
	358	de_dxi(num_gp) = de_dx(k+1,j+1,i)
	359	de_dyi(num_gp) = de_dy(k+1,j+1,i)
	360	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	361	ENDIF
	362
	363	IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) &
	364	THEN
	365	num_gp = num_gp + 1
	366	gp_outside_of_building(5) = 1
	367	location(num_gp,1) = (i+1) * dx
	368	location(num_gp,2) = j * dy
	369	location(num_gp,3) = k * dz - 0.5 * dz
	370	ei(num_gp) = e(k,j,i+1)
	371	dissi(num_gp) = diss(k,j,i+1)
	372	de_dxi(num_gp) = de_dx(k,j,i+1)
	373	de_dyi(num_gp) = de_dy(k,j,i+1)
	374	de_dzi(num_gp) = de_dz(k,j,i+1)
	375	ENDIF
	376
	377	IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) &
	378	THEN
	379	num_gp = num_gp + 1
	380	gp_outside_of_building(6) = 1
	381	location(num_gp,1) = (i+1) * dx
	382	location(num_gp,2) = (j+1) * dy
	383	location(num_gp,3) = k * dz - 0.5 * dz
	384	ei(num_gp) = e(k,j+1,i+1)
	385	dissi(num_gp) = diss(k,j+1,i+1)
	386	de_dxi(num_gp) = de_dx(k,j+1,i+1)
	387	de_dyi(num_gp) = de_dy(k,j+1,i+1)
	388	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	389	ENDIF
	390
	391	IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) &
	392	THEN
	393	num_gp = num_gp + 1
	394	gp_outside_of_building(7) = 1
	395	location(num_gp,1) = (i+1) * dx
	396	location(num_gp,2) = j * dy
	397	location(num_gp,3) = (k+1) * dz - 0.5 * dz
	398	ei(num_gp) = e(k+1,j,i+1)
	399	dissi(num_gp) = diss(k+1,j,i+1)
	400	de_dxi(num_gp) = de_dx(k+1,j,i+1)
	401	de_dyi(num_gp) = de_dy(k+1,j,i+1)
	402	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	403	ENDIF
	404
	405	IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)&
	406	THEN
	407	num_gp = num_gp + 1
	408	gp_outside_of_building(8) = 1
	409	location(num_gp,1) = (i+1) * dx
	410	location(num_gp,2) = (j+1) * dy
	411	location(num_gp,3) = (k+1) * dz - 0.5 * dz
	412	ei(num_gp) = e(k+1,j+1,i+1)
	413	dissi(num_gp) = diss(k+1,j+1,i+1)
	414	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
	415	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
	416	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	417	ENDIF
	418
	419	!
	420	!-- If all gridpoints are situated outside of a building, then the
	421	!-- ordinary interpolation scheme can be used.
	422	IF ( num_gp == 8 ) THEN
	423
	424	x = particles(n)%x - i * dx
	425	y = particles(n)%y - j * dy
	426	aa = x2 + y2
	427	bb = ( dx - x )2 + y2
	428	cc = x2 + ( dy - y )2
	429	dd = ( dx - x )2 + ( dy - y )2
	430	gg = aa + bb + cc + dd
	431
	432	e_int_l = (( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
	433	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)&
	434	) / ( 3.0 * gg )
	435
	436	IF ( k+1 == nzt+1 ) THEN
	437	e_int = e_int_l
	438	ELSE
	439	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
	440	( gg - bb ) * e(k+1,j,i+1) + &
	441	( gg - cc ) * e(k+1,j+1,i) + &
	442	( gg - dd ) * e(k+1,j+1,i+1) &
	443	) / ( 3.0 * gg )
	444	e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
	445	( e_int_u - e_int_l )
	446	ENDIF
	447
	448	!
	449	!-- Interpolate the TKE gradient along x (adopt incides i,j,k
	450	!-- and all position variables from above (TKE))
	451	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
	452	( gg - bb ) * de_dx(k,j,i+1) + &
	453	( gg - cc ) * de_dx(k,j+1,i) + &
	454	( gg - dd ) * de_dx(k,j+1,i+1) &
	455	) / ( 3.0 * gg )
	456
	457	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	458	de_dx_int = de_dx_int_l
	459	ELSE
	460	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
	461	( gg - bb ) * de_dx(k+1,j,i+1) + &
	462	( gg - cc ) * de_dx(k+1,j+1,i) + &
	463	( gg - dd ) * de_dx(k+1,j+1,i+1) &
	464	) / ( 3.0 * gg )
	465	de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / &
	466	dz * ( de_dx_int_u - de_dx_int_l )
	467	ENDIF
	468
	469	!
	470	!-- Interpolate the TKE gradient along y
	471	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
	472	( gg - bb ) * de_dy(k,j,i+1) + &
	473	( gg - cc ) * de_dy(k,j+1,i) + &
	474	( gg - dd ) * de_dy(k,j+1,i+1) &
	475	) / ( 3.0 * gg )
	476	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	477	de_dy_int = de_dy_int_l
	478	ELSE
	479	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
	480	( gg - bb ) * de_dy(k+1,j,i+1) + &
	481	( gg - cc ) * de_dy(k+1,j+1,i) + &
	482	( gg - dd ) * de_dy(k+1,j+1,i+1) &
	483	) / ( 3.0 * gg )
	484	de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / &
	485	dz * ( de_dy_int_u - de_dy_int_l )
	486	ENDIF
	487
	488	!
	489	!-- Interpolate the TKE gradient along z
	490	IF ( particles(n)%z < 0.5 * dz ) THEN
	491	de_dz_int = 0.0
	492	ELSE
	493	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
	494	( gg - bb ) * de_dz(k,j,i+1) + &
	495	( gg - cc ) * de_dz(k,j+1,i) + &
	496	( gg - dd ) * de_dz(k,j+1,i+1) &
	497	) / ( 3.0 * gg )
	498
	499	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
	500	de_dz_int = de_dz_int_l
	501	ELSE
	502	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
	503	( gg - bb ) * de_dz(k+1,j,i+1) + &
	504	( gg - cc ) * de_dz(k+1,j+1,i) + &
	505	( gg - dd ) * de_dz(k+1,j+1,i+1) &
	506	) / ( 3.0 * gg )
	507	de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /&
	508	dz * ( de_dz_int_u - de_dz_int_l )
	509	ENDIF
	510	ENDIF
	511
	512	!
	513	!-- Interpolate the dissipation of TKE
	514	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
	515	( gg - bb ) * diss(k,j,i+1) + &
	516	( gg - cc ) * diss(k,j+1,i) + &
	517	( gg - dd ) * diss(k,j+1,i+1) &
	518	) / ( 3.0 * gg )
	519
	520	IF ( k+1 == nzt+1 ) THEN
	521	diss_int = diss_int_l
	522	ELSE
	523	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
	524	( gg - bb ) * diss(k+1,j,i+1) + &
	525	( gg - cc ) * diss(k+1,j+1,i) + &
	526	( gg - dd ) * diss(k+1,j+1,i+1) &
	527	) / ( 3.0 * gg )
	528	diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *&
	529	( diss_int_u - diss_int_l )
	530	ENDIF
	531
	532	ELSE
	533
	534	!
	535	!-- If wall between gridpoint 1 and gridpoint 5, then
	536	!-- Neumann boundary condition has to be applied
	537	IF ( gp_outside_of_building(1) == 1 .AND. &
	538	gp_outside_of_building(5) == 0 ) THEN
	539	num_gp = num_gp + 1
	540	location(num_gp,1) = i * dx + 0.5 * dx
	541	location(num_gp,2) = j * dy
	542	location(num_gp,3) = k * dz - 0.5 * dz
	543	ei(num_gp) = e(k,j,i)
	544	dissi(num_gp) = diss(k,j,i)
	545	de_dxi(num_gp) = 0.0
	546	de_dyi(num_gp) = de_dy(k,j,i)
	547	de_dzi(num_gp) = de_dz(k,j,i)
	548	ENDIF
	549
	550	IF ( gp_outside_of_building(5) == 1 .AND. &
	551	gp_outside_of_building(1) == 0 ) THEN
	552	num_gp = num_gp + 1
	553	location(num_gp,1) = i * dx + 0.5 * dx
	554	location(num_gp,2) = j * dy
	555	location(num_gp,3) = k * dz - 0.5 * dz
	556	ei(num_gp) = e(k,j,i+1)
	557	dissi(num_gp) = diss(k,j,i+1)
	558	de_dxi(num_gp) = 0.0
	559	de_dyi(num_gp) = de_dy(k,j,i+1)
	560	de_dzi(num_gp) = de_dz(k,j,i+1)
	561	ENDIF
	562
	563	!
	564	!-- If wall between gridpoint 5 and gridpoint 6, then
	565	!-- then Neumann boundary condition has to be applied
	566	IF ( gp_outside_of_building(5) == 1 .AND. &
	567	gp_outside_of_building(6) == 0 ) THEN
	568	num_gp = num_gp + 1
	569	location(num_gp,1) = (i+1) * dx
	570	location(num_gp,2) = j * dy + 0.5 * dy
	571	location(num_gp,3) = k * dz - 0.5 * dz
	572	ei(num_gp) = e(k,j,i+1)
	573	dissi(num_gp) = diss(k,j,i+1)
	574	de_dxi(num_gp) = de_dx(k,j,i+1)
	575	de_dyi(num_gp) = 0.0
	576	de_dzi(num_gp) = de_dz(k,j,i+1)
	577	ENDIF
	578
	579	IF ( gp_outside_of_building(6) == 1 .AND. &
	580	gp_outside_of_building(5) == 0 ) THEN
	581	num_gp = num_gp + 1
	582	location(num_gp,1) = (i+1) * dx
	583	location(num_gp,2) = j * dy + 0.5 * dy
	584	location(num_gp,3) = k * dz - 0.5 * dz
	585	ei(num_gp) = e(k,j+1,i+1)
	586	dissi(num_gp) = diss(k,j+1,i+1)
	587	de_dxi(num_gp) = de_dx(k,j+1,i+1)
	588	de_dyi(num_gp) = 0.0
	589	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	590	ENDIF
	591
	592	!
	593	!-- If wall between gridpoint 2 and gridpoint 6, then
	594	!-- Neumann boundary condition has to be applied
	595	IF ( gp_outside_of_building(2) == 1 .AND. &
	596	gp_outside_of_building(6) == 0 ) THEN
	597	num_gp = num_gp + 1
	598	location(num_gp,1) = i * dx + 0.5 * dx
	599	location(num_gp,2) = (j+1) * dy
	600	location(num_gp,3) = k * dz - 0.5 * dz
	601	ei(num_gp) = e(k,j+1,i)
	602	dissi(num_gp) = diss(k,j+1,i)
	603	de_dxi(num_gp) = 0.0
	604	de_dyi(num_gp) = de_dy(k,j+1,i)
	605	de_dzi(num_gp) = de_dz(k,j+1,i)
	606	ENDIF
	607
	608	IF ( gp_outside_of_building(6) == 1 .AND. &
	609	gp_outside_of_building(2) == 0 ) THEN
	610	num_gp = num_gp + 1
	611	location(num_gp,1) = i * dx + 0.5 * dx
	612	location(num_gp,2) = (j+1) * dy
	613	location(num_gp,3) = k * dz - 0.5 * dz
	614	ei(num_gp) = e(k,j+1,i+1)
	615	dissi(num_gp) = diss(k,j+1,i+1)
	616	de_dxi(num_gp) = 0.0
	617	de_dyi(num_gp) = de_dy(k,j+1,i+1)
	618	de_dzi(num_gp) = de_dz(k,j+1,i+1)
	619	ENDIF
	620
	621	!
	622	!-- If wall between gridpoint 1 and gridpoint 2, then
	623	!-- Neumann boundary condition has to be applied
	624	IF ( gp_outside_of_building(1) == 1 .AND. &
	625	gp_outside_of_building(2) == 0 ) THEN
	626	num_gp = num_gp + 1
	627	location(num_gp,1) = i * dx
	628	location(num_gp,2) = j * dy + 0.5 * dy
	629	location(num_gp,3) = k * dz - 0.5 * dz
	630	ei(num_gp) = e(k,j,i)
	631	dissi(num_gp) = diss(k,j,i)
	632	de_dxi(num_gp) = de_dx(k,j,i)
	633	de_dyi(num_gp) = 0.0
	634	de_dzi(num_gp) = de_dz(k,j,i)
	635	ENDIF
	636
	637	IF ( gp_outside_of_building(2) == 1 .AND. &
	638	gp_outside_of_building(1) == 0 ) THEN
	639	num_gp = num_gp + 1
	640	location(num_gp,1) = i * dx
	641	location(num_gp,2) = j * dy + 0.5 * dy
	642	location(num_gp,3) = k * dz - 0.5 * dz
	643	ei(num_gp) = e(k,j+1,i)
	644	dissi(num_gp) = diss(k,j+1,i)
	645	de_dxi(num_gp) = de_dx(k,j+1,i)
	646	de_dyi(num_gp) = 0.0
	647	de_dzi(num_gp) = de_dz(k,j+1,i)
	648	ENDIF
	649
	650	!
	651	!-- If wall between gridpoint 3 and gridpoint 7, then
	652	!-- Neumann boundary condition has to be applied
	653	IF ( gp_outside_of_building(3) == 1 .AND. &
	654	gp_outside_of_building(7) == 0 ) THEN
	655	num_gp = num_gp + 1
	656	location(num_gp,1) = i * dx + 0.5 * dx
	657	location(num_gp,2) = j * dy
	658	location(num_gp,3) = k * dz + 0.5 * dz
	659	ei(num_gp) = e(k+1,j,i)
	660	dissi(num_gp) = diss(k+1,j,i)
	661	de_dxi(num_gp) = 0.0
	662	de_dyi(num_gp) = de_dy(k+1,j,i)
	663	de_dzi(num_gp) = de_dz(k+1,j,i)
	664	ENDIF
	665
	666	IF ( gp_outside_of_building(7) == 1 .AND. &
	667	gp_outside_of_building(3) == 0 ) THEN
	668	num_gp = num_gp + 1
	669	location(num_gp,1) = i * dx + 0.5 * dx
	670	location(num_gp,2) = j * dy
	671	location(num_gp,3) = k * dz + 0.5 * dz
	672	ei(num_gp) = e(k+1,j,i+1)
	673	dissi(num_gp) = diss(k+1,j,i+1)
	674	de_dxi(num_gp) = 0.0
	675	de_dyi(num_gp) = de_dy(k+1,j,i+1)
	676	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	677	ENDIF
	678
	679	!
	680	!-- If wall between gridpoint 7 and gridpoint 8, then
	681	!-- Neumann boundary condition has to be applied
	682	IF ( gp_outside_of_building(7) == 1 .AND. &
	683	gp_outside_of_building(8) == 0 ) THEN
	684	num_gp = num_gp + 1
	685	location(num_gp,1) = (i+1) * dx
	686	location(num_gp,2) = j * dy + 0.5 * dy
	687	location(num_gp,3) = k * dz + 0.5 * dz
	688	ei(num_gp) = e(k+1,j,i+1)
	689	dissi(num_gp) = diss(k+1,j,i+1)
	690	de_dxi(num_gp) = de_dx(k+1,j,i+1)
	691	de_dyi(num_gp) = 0.0
	692	de_dzi(num_gp) = de_dz(k+1,j,i+1)
	693	ENDIF
	694
	695	IF ( gp_outside_of_building(8) == 1 .AND. &
	696	gp_outside_of_building(7) == 0 ) THEN
	697	num_gp = num_gp + 1
	698	location(num_gp,1) = (i+1) * dx
	699	location(num_gp,2) = j * dy + 0.5 * dy
	700	location(num_gp,3) = k * dz + 0.5 * dz
	701	ei(num_gp) = e(k+1,j+1,i+1)
	702	dissi(num_gp) = diss(k+1,j+1,i+1)
	703	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
	704	de_dyi(num_gp) = 0.0
	705	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	706	ENDIF
	707
	708	!
	709	!-- If wall between gridpoint 4 and gridpoint 8, then
	710	!-- Neumann boundary condition has to be applied
	711	IF ( gp_outside_of_building(4) == 1 .AND. &
	712	gp_outside_of_building(8) == 0 ) THEN
	713	num_gp = num_gp + 1
	714	location(num_gp,1) = i * dx + 0.5 * dx
	715	location(num_gp,2) = (j+1) * dy
	716	location(num_gp,3) = k * dz + 0.5 * dz
	717	ei(num_gp) = e(k+1,j+1,i)
	718	dissi(num_gp) = diss(k+1,j+1,i)
	719	de_dxi(num_gp) = 0.0
	720	de_dyi(num_gp) = de_dy(k+1,j+1,i)
	721	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	722	ENDIF
	723
	724	IF ( gp_outside_of_building(8) == 1 .AND. &
	725	gp_outside_of_building(4) == 0 ) THEN
	726	num_gp = num_gp + 1
	727	location(num_gp,1) = i * dx + 0.5 * dx
	728	location(num_gp,2) = (j+1) * dy
	729	location(num_gp,3) = k * dz + 0.5 * dz
	730	ei(num_gp) = e(k+1,j+1,i+1)
	731	dissi(num_gp) = diss(k+1,j+1,i+1)
	732	de_dxi(num_gp) = 0.0
	733	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
	734	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
	735	ENDIF
	736
	737	!
	738	!-- If wall between gridpoint 3 and gridpoint 4, then
	739	!-- Neumann boundary condition has to be applied
	740	IF ( gp_outside_of_building(3) == 1 .AND. &
	741	gp_outside_of_building(4) == 0 ) THEN
	742	num_gp = num_gp + 1
	743	location(num_gp,1) = i * dx
	744	location(num_gp,2) = j * dy + 0.5 * dy
	745	location(num_gp,3) = k * dz + 0.5 * dz
	746	ei(num_gp) = e(k+1,j,i)
	747	dissi(num_gp) = diss(k+1,j,i)
	748	de_dxi(num_gp) = de_dx(k+1,j,i)
	749	de_dyi(num_gp) = 0.0
	750	de_dzi(num_gp) = de_dz(k+1,j,i)
	751	ENDIF
	752
	753	IF ( gp_outside_of_building(4) == 1 .AND. &
	754	gp_outside_of_building(3) == 0 ) THEN
	755	num_gp = num_gp + 1
	756	location(num_gp,1) = i * dx
	757	location(num_gp,2) = j * dy + 0.5 * dy
	758	location(num_gp,3) = k * dz + 0.5 * dz
	759	ei(num_gp) = e(k+1,j+1,i)
	760	dissi(num_gp) = diss(k+1,j+1,i)
	761	de_dxi(num_gp) = de_dx(k+1,j+1,i)
	762	de_dyi(num_gp) = 0.0
	763	de_dzi(num_gp) = de_dz(k+1,j+1,i)
	764	ENDIF
	765
	766	!
	767	!-- If wall between gridpoint 1 and gridpoint 3, then
	768	!-- Neumann boundary condition has to be applied
	769	!-- (only one case as only building beneath is possible)
	770	IF ( gp_outside_of_building(1) == 0 .AND. &
	771	gp_outside_of_building(3) == 1 ) THEN
	772	num_gp = num_gp + 1
	773	location(num_gp,1) = i * dx
	774	location(num_gp,2) = j * dy
	775	location(num_gp,3) = k * dz
	776	ei(num_gp) = e(k+1,j,i)
	777	dissi(num_gp) = diss(k+1,j,i)
	778	de_dxi(num_gp) = de_dx(k+1,j,i)
	779	de_dyi(num_gp) = de_dy(k+1,j,i)
	780	de_dzi(num_gp) = 0.0
	781	ENDIF
	782
	783	!
	784	!-- If wall between gridpoint 5 and gridpoint 7, then
	785	!-- Neumann boundary condition has to be applied
	786	!-- (only one case as only building beneath is possible)
	787	IF ( gp_outside_of_building(5) == 0 .AND. &
	788	gp_outside_of_building(7) == 1 ) THEN
	789	num_gp = num_gp + 1
	790	location(num_gp,1) = (i+1) * dx
	791	location(num_gp,2) = j * dy
	792	location(num_gp,3) = k * dz
	793	ei(num_gp) = e(k+1,j,i+1)
	794	dissi(num_gp) = diss(k+1,j,i+1)
	795	de_dxi(num_gp) = de_dx(k+1,j,i+1)
	796	de_dyi(num_gp) = de_dy(k+1,j,i+1)
	797	de_dzi(num_gp) = 0.0
	798	ENDIF
	799
	800	!
	801	!-- If wall between gridpoint 2 and gridpoint 4, then
	802	!-- Neumann boundary condition has to be applied
	803	!-- (only one case as only building beneath is possible)
	804	IF ( gp_outside_of_building(2) == 0 .AND. &
	805	gp_outside_of_building(4) == 1 ) THEN
	806	num_gp = num_gp + 1
	807	location(num_gp,1) = i * dx
	808	location(num_gp,2) = (j+1) * dy
	809	location(num_gp,3) = k * dz
	810	ei(num_gp) = e(k+1,j+1,i)
	811	dissi(num_gp) = diss(k+1,j+1,i)
	812	de_dxi(num_gp) = de_dx(k+1,j+1,i)
	813	de_dyi(num_gp) = de_dy(k+1,j+1,i)
	814	de_dzi(num_gp) = 0.0
	815	ENDIF
	816
	817	!
	818	!-- If wall between gridpoint 6 and gridpoint 8, then
	819	!-- Neumann boundary condition has to be applied
	820	!-- (only one case as only building beneath is possible)
	821	IF ( gp_outside_of_building(6) == 0 .AND. &
	822	gp_outside_of_building(8) == 1 ) THEN
	823	num_gp = num_gp + 1
	824	location(num_gp,1) = (i+1) * dx
	825	location(num_gp,2) = (j+1) * dy
	826	location(num_gp,3) = k * dz
	827	ei(num_gp) = e(k+1,j+1,i+1)
	828	dissi(num_gp) = diss(k+1,j+1,i+1)
	829	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
	830	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
	831	de_dzi(num_gp) = 0.0
	832	ENDIF
	833
	834	!
	835	!-- Carry out the interpolation
	836	IF ( num_gp == 1 ) THEN
	837	!
	838	!-- If only one of the gridpoints is situated outside of the
	839	!-- building, it follows that the values at the particle
	840	!-- location are the same as the gridpoint values
	841	e_int = ei(num_gp)
	842	diss_int = dissi(num_gp)
	843	de_dx_int = de_dxi(num_gp)
	844	de_dy_int = de_dyi(num_gp)
	845	de_dz_int = de_dzi(num_gp)
	846	ELSE IF ( num_gp > 1 ) THEN
	847
	848	d_sum = 0.0
	849	!
	850	!-- Evaluation of the distances between the gridpoints
	851	!-- contributing to the interpolated values, and the particle
	852	!-- location
	853	DO agp = 1, num_gp
	854	d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 &
	855	+ ( particles(n)%y-location(agp,2) )**2 &
	856	+ ( particles(n)%z-location(agp,3) )**2
	857	d_sum = d_sum + d_gp_pl(agp)
	858	ENDDO
	859
	860	!
	861	!-- Finally the interpolation can be carried out
	862	e_int = 0.0
	863	diss_int = 0.0
	864	de_dx_int = 0.0
	865	de_dy_int = 0.0
	866	de_dz_int = 0.0
	867	DO agp = 1, num_gp
	868	e_int = e_int + ( d_sum - d_gp_pl(agp) ) * &
	869	ei(agp) / ( (num_gp-1) * d_sum )
	870	diss_int = diss_int + ( d_sum - d_gp_pl(agp) ) * &
	871	dissi(agp) / ( (num_gp-1) * d_sum )
	872	de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * &
	873	de_dxi(agp) / ( (num_gp-1) * d_sum )
	874	de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * &
	875	de_dyi(agp) / ( (num_gp-1) * d_sum )
	876	de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * &
	877	de_dzi(agp) / ( (num_gp-1) * d_sum )
	878	ENDDO
	879
	880	ENDIF
	881
	882	ENDIF
	883
	884	ENDIF
	885
	886	!
	887	!-- Vertically interpolate the horizontally averaged SGS TKE and
	888	!-- resolved-scale velocity variances and use the interpolated values
	889	!-- to calculate the coefficient fs, which is a measure of the ratio
	890	!-- of the subgrid-scale turbulent kinetic energy to the total amount
	891	!-- of turbulent kinetic energy.
	892	IF ( k == 0 ) THEN
	893	e_mean_int = hom(0,1,8,0)
	894	ELSE
	895	e_mean_int = hom(k,1,8,0) + &
	896	( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
	897	( zu(k+1) - zu(k) ) * &
	898	( particles(n)%z - zu(k) )
	899	ENDIF
	900
	901	kw = particles(n)%z / dz
	902
	903	IF ( k == 0 ) THEN
	904	aa = hom(k+1,1,30,0) * ( particles(n)%z / &
	905	( 0.5 * ( zu(k+1) - zu(k) ) ) )
	906	bb = hom(k+1,1,31,0) * ( particles(n)%z / &
	907	( 0.5 * ( zu(k+1) - zu(k) ) ) )
	908	cc = hom(kw+1,1,32,0) * ( particles(n)%z / &
	909	( 1.0 * ( zw(kw+1) - zw(kw) ) ) )
	910	ELSE
	911	aa = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
	912	( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
	913	bb = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
	914	( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
	915	cc = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *&
	916	( ( particles(n)%z - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
	917	ENDIF
	918
	919	vv_int = ( 1.0 / 3.0 ) * ( aa + bb + cc )
	920
	921	fs_int = ( 2.0 / 3.0 ) * e_mean_int / &
	922	( vv_int + ( 2.0 / 3.0 ) * e_mean_int )
	923
	924	!
	925	!-- Calculate the Lagrangian timescale according to Weil et al. (2004).
	926	lagr_timescale = ( 4.0 * e_int ) / &
	927	( 3.0 * fs_int * c_0 * diss_int )
	928
	929	!
	930	!-- Calculate the next particle timestep. dt_gap is the time needed to
	931	!-- complete the current LES timestep.
	932	dt_gap = dt_3d - particles(n)%dt_sum
	933	dt_particle = MIN( dt_3d, 0.025 * lagr_timescale, dt_gap )
	934
	935	!
	936	!-- The particle timestep should not be too small in order to prevent
	937	!-- the number of particle timesteps of getting too large
	938	IF ( dt_particle < dt_min_part .AND. dt_min_part < dt_gap ) &
	939	THEN
	940	dt_particle = dt_min_part
	941	ENDIF
	942
	943	!
	944	!-- Calculate the SGS velocity components
	945	IF ( particles(n)%age == 0.0 ) THEN
	946	!
	947	!-- For new particles the SGS components are derived from the SGS
	948	!-- TKE. Limit the Gaussian random number to the interval
	949	!-- [-5.0sigma, 5.0sigma] in order to prevent the SGS velocities
	950	!-- from becoming unrealistically large.
	951	particles(n)%rvar1 = SQRT( 2.0 * sgs_wfu_part * e_int ) * &
	952	( random_gauss( iran_part, 5.0 ) - 1.0 )
	953	particles(n)%rvar2 = SQRT( 2.0 * sgs_wfv_part * e_int ) * &
	954	( random_gauss( iran_part, 5.0 ) - 1.0 )
	955	particles(n)%rvar3 = SQRT( 2.0 * sgs_wfw_part * e_int ) * &
	956	( random_gauss( iran_part, 5.0 ) - 1.0 )
	957
	958	ELSE
	959
	960	!
	961	!-- Restriction of the size of the new timestep: compared to the
	962	!-- previous timestep the increase must not exceed 200%
	963
	964	dt_particle_m = particles(n)%age - particles(n)%age_m
	965	IF ( dt_particle > 2.0 * dt_particle_m ) THEN
	966	dt_particle = 2.0 * dt_particle_m
	967	ENDIF
	968
	969	!
	970	!-- For old particles the SGS components are correlated with the
	971	!-- values from the previous timestep. Random numbers have also to
	972	!-- be limited (see above).
	973	!-- As negative values for the subgrid TKE are not allowed, the
	974	!-- change of the subgrid TKE with time cannot be smaller than
	975	!-- -e_int/dt_particle. This value is used as a lower boundary
	976	!-- value for the change of TKE
	977
	978	de_dt_min = - e_int / dt_particle
	979
	980	de_dt = ( e_int - particles(n)%e_m ) / dt_particle_m
	981
	982	IF ( de_dt < de_dt_min ) THEN
	983	de_dt = de_dt_min
	984	ENDIF
	985
	986	particles(n)%rvar1 = particles(n)%rvar1 - fs_int * c_0 * &
	987	diss_int * particles(n)%rvar1 * dt_particle /&
	988	( 4.0 * sgs_wfu_part * e_int ) + &
	989	( 2.0 * sgs_wfu_part * de_dt * &
	990	particles(n)%rvar1 / &
	991	( 2.0 * sgs_wfu_part * e_int ) + de_dx_int &
	992	) * dt_particle / 2.0 + &
	993	SQRT( fs_int * c_0 * diss_int ) * &
	994	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
	995	SQRT( dt_particle )
	996
	997	particles(n)%rvar2 = particles(n)%rvar2 - fs_int * c_0 * &
	998	diss_int * particles(n)%rvar2 * dt_particle /&
	999	( 4.0 * sgs_wfv_part * e_int ) + &
	1000	( 2.0 * sgs_wfv_part * de_dt * &
	1001	particles(n)%rvar2 / &
	1002	( 2.0 * sgs_wfv_part * e_int ) + de_dy_int &
	1003	) * dt_particle / 2.0 + &
	1004	SQRT( fs_int * c_0 * diss_int ) * &
	1005	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
	1006	SQRT( dt_particle )
	1007
	1008	particles(n)%rvar3 = particles(n)%rvar3 - fs_int * c_0 * &
	1009	diss_int * particles(n)%rvar3 * dt_particle /&
	1010	( 4.0 * sgs_wfw_part * e_int ) + &
	1011	( 2.0 * sgs_wfw_part * de_dt * &
	1012	particles(n)%rvar3 / &
	1013	( 2.0 * sgs_wfw_part * e_int ) + de_dz_int &
	1014	) * dt_particle / 2.0 + &
	1015	SQRT( fs_int * c_0 * diss_int ) * &
	1016	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
	1017	SQRT( dt_particle )
	1018
	1019	ENDIF
	1020
	1021	u_int = u_int + particles(n)%rvar1
	1022	v_int = v_int + particles(n)%rvar2
	1023	w_int = w_int + particles(n)%rvar3
	1024
	1025	!
	1026	!-- Store the SGS TKE of the current timelevel which is needed for
	1027	!-- for calculating the SGS particle velocities at the next timestep
	1028	particles(n)%e_m = e_int
	1029
	1030	ELSE
	1031	!
	1032	!-- If no SGS velocities are used, only the particle timestep has to
	1033	!-- be set
	1034	dt_particle = dt_3d
	1035
	1036	ENDIF
	1037
	1038	!
	1039	!-- Store the old age of the particle ( needed to prevent that a
	1040	!-- particle crosses several PEs during one timestep, and for the
	1041	!-- evaluation of the subgrid particle velocity fluctuations )
	1042	particles(n)%age_m = particles(n)%age
	1043
	1044
	1045	!
	1046	!-- Particle advection
	1047	IF ( particle_groups(particles(n)%group)%density_ratio == 0.0 ) THEN
	1048	!
	1049	!-- Pure passive transport (without particle inertia)
	1050	particles(n)%x = particles(n)%x + u_int * dt_particle
	1051	particles(n)%y = particles(n)%y + v_int * dt_particle
	1052	particles(n)%z = particles(n)%z + w_int * dt_particle
	1053
	1054	particles(n)%speed_x = u_int
	1055	particles(n)%speed_y = v_int
	1056	particles(n)%speed_z = w_int
	1057
	1058	ELSE
	1059	!
	1060	!-- Transport of particles with inertia
	1061	particles(n)%x = particles(n)%x + particles(n)%speed_x * &
	1062	dt_particle
	1063	particles(n)%y = particles(n)%y + particles(n)%speed_y * &
	1064	dt_particle
	1065	particles(n)%z = particles(n)%z + particles(n)%speed_z * &
	1066	dt_particle
	1067
	1068	!
	1069	!-- Update of the particle velocity
	1070	dens_ratio = particle_groups(particles(n)%group)%density_ratio
	1071	IF ( cloud_droplets ) THEN
	1072	exp_arg = 4.5 * dens_ratio * molecular_viscosity / &
	1073	( particles(n)%radius )*2 &
	1074	( 1.0 + 0.15 * ( 2.0 * particles(n)%radius * &
	1075	SQRT( ( u_int - particles(n)%speed_x )**2 + &
	1076	( v_int - particles(n)%speed_y )**2 + &
	1077	( w_int - particles(n)%speed_z )**2 ) / &
	1078	molecular_viscosity )**0.687 &
	1079	)
	1080	exp_term = EXP( -exp_arg * dt_particle )
	1081	ELSEIF ( use_sgs_for_particles ) THEN
	1082	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1083	exp_term = EXP( -exp_arg * dt_particle )
	1084	ELSE
	1085	exp_arg = particle_groups(particles(n)%group)%exp_arg
	1086	exp_term = particle_groups(particles(n)%group)%exp_term
	1087	ENDIF
	1088	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
	1089	u_int * ( 1.0 - exp_term )
	1090	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
	1091	v_int * ( 1.0 - exp_term )
	1092	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
	1093	( w_int - ( 1.0 - dens_ratio ) * g / exp_arg )&
	1094	* ( 1.0 - exp_term )
	1095	ENDIF
	1096
	1097	!
	1098	!-- Increment the particle age and the total time that the particle
	1099	!-- has advanced within the particle timestep procedure
	1100	particles(n)%age = particles(n)%age + dt_particle
	1101	particles(n)%dt_sum = particles(n)%dt_sum + dt_particle
	1102
	1103	!
	1104	!-- Check whether there is still a particle that has not yet completed
	1105	!-- the total LES timestep
	1106	IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8 ) THEN
	1107	dt_3d_reached_l = .FALSE.
	1108	ENDIF
	1109
	1110	ENDDO
	1111
	1112
	1113	END SUBROUTINE lpm_advec

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |