Home

Context Navigation

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

Last change on this file since 2177 was 2101, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 73.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |