Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |