Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |