Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |