Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |