Home

Context Navigation

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

Last change on this file since 2100 was 2100, checked in by suehring, 7 years ago
Prevent extremely large SGS-velocities by limiting TKE in stochastic terms.
Property svn:keywords set to `Id`
File size: 73.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |