Home

Context Navigation

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

Last change on this file since 2317 was 2317, checked in by suehring, 7 years ago
get topograpyh top index via function call
Property svn:keywords set to `Id`
File size: 76.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |