Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |