Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |