Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |