Home

Context Navigation

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

Last change on this file since 2549 was 2417, checked in by suehring, 7 years ago
Major bugfix in modeling SGS particle speeds.
Property svn:keywords set to `Id`
File size: 80.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |