Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |