Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |