Home

Context Navigation

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

Last change on this file since 1385 was 1370, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 61.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |