Home

Context Navigation

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

Last change on this file since 1889 was 1889, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 65.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |