Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |