Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |