Home

Context Navigation

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

Last change on this file since 1936 was 1936, checked in by suehring, 8 years ago
deallocation of unused particle memory, formatting adjustments
Property svn:keywords set to `Id`
File size: 72.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |