Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |