Home

Context Navigation

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

Last change on this file since 2345 was 2318, checked in by suehring, 7 years ago
get topograpyhy top index via function call
Property svn:keywords set to `Id`
File size: 77.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |