Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |