Home

Context Navigation

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

Last change on this file since 957 was 850, checked in by raasch, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 46.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |