Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |