Home

Context Navigation

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

Last change on this file since 2000 was 2000, checked in by knoop, 8 years ago
Forced header and separation lines into 80 columns
Property svn:keywords set to `Id`
File size: 72.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |