Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |