Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |