Home

Context Navigation

lpm_advec.f90 @ 849

Last change on this file since 849 was 849, checked in by raasch, 12 years ago

Changed:

Original routine advec_particles split into several new subroutines and renamed
lpm.
init_particles renamed lpm_init
user_advec_particles renamed user_lpm_advec,
particle_boundary_conds renamed lpm_boundary_conds,
set_particle_attributes renamed lpm_set_attributes,
user_init_particles renamed user_lpm_init,
user_particle_attributes renamed user_lpm_set_attributes
(Makefile, lpm_droplet_collision, lpm_droplet_condensation, init_3d_model, modules, palm, read_var_list, time_integration, write_var_list, deleted: advec_particles, init_particles, particle_boundary_conds, set_particle_attributes, user_advec_particles, user_init_particles, user_particle_attributes, new: lpm, lpm_advec, lpm_boundary_conds, lpm_calc_liquid_water_content, lpm_data_output_particles, lpm_droplet_collision, lpm_drollet_condensation, lpm_exchange_horiz, lpm_extend_particle_array, lpm_extend_tails, lpm_extend_tail_array, lpm_init, lpm_init_sgs_tke, lpm_pack_arrays, lpm_read_restart_file, lpm_release_set, lpm_set_attributes, lpm_sort_arrays, lpm_write_exchange_statistics, lpm_write_restart_file, user_lpm_advec, user_lpm_init, user_lpm_set_attributes

Property svn:keywords set to Id

File size: 46.7 KB

Line
1	SUBROUTINE lpm_advec
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! ------------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: lpm_advec.f90 849 2012-03-15 10:35:09Z raasch $
11	!
12	!
13	! Description:
14	! ------------
15	! Calculation of new particle positions due to advection using a simple Euler
16	! scheme. Particles may feel inertia effects. SGS transport can be included
17	! using the stochastic model of Weil et al. (2004, JAS, 61, 2877-2887).
18	!------------------------------------------------------------------------------!
19
20	USE arrays_3d
21	USE control_parameters
22	USE grid_variables
23	USE indices
24	USE particle_attributes
25	USE statistics
26
27	IMPLICIT NONE
28
29	INTEGER :: i, j, k, n
30
31	REAL :: aa, bb, cc, dd, dens_ratio, exp_arg, exp_term, gg, u_int, &
32	u_int_l, u_int_u, v_int, v_int_l, v_int_u, w_int, w_int_l, &
33	w_int_u, x, y
34
35
36	INTEGER :: agp, kw, num_gp
37	INTEGER :: gp_outside_of_building(1:8)
38
39	REAL :: d_sum, de_dx_int, de_dx_int_l, de_dx_int_u, de_dy_int, &
40	de_dy_int_l, de_dy_int_u, de_dt, de_dt_min, de_dz_int, &
41	de_dz_int_l, de_dz_int_u, diss_int, diss_int_l, diss_int_u, &
42	dt_gap, dt_particle, dt_particle_m, e_int, e_int_l, e_int_u, &
43	e_mean_int, fs_int, lagr_timescale, random_gauss, vv_int
44
45	REAL :: location(1:30,1:3)
46	REAL, DIMENSION(1:30) :: de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
47
48
49	DO n = 1, number_of_particles
50
51	!
52	!-- Move particle only if the LES timestep has not (approximately) been
53	!-- reached
54	IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 ) CYCLE
55
56	!
57	!-- Interpolate u velocity-component, determine left, front, bottom
58	!-- index of u-array
59	i = ( particles(n)%x + 0.5 * dx ) * ddx
60	j = particles(n)%y * ddy
61	k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz &
62	+ offset_ocean_nzt ! only exact if equidistant
63
64	!
65	!-- Interpolation of the velocity components in the xy-plane
66	x = particles(n)%x + ( 0.5 - i ) * dx
67	y = particles(n)%y - j * dy
68	aa = x2 + y2
69	bb = ( dx - x )2 + y2
70	cc = x2 + ( dy - y )2
71	dd = ( dx - x )2 + ( dy - y )2
72	gg = aa + bb + cc + dd
73
74	u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) &
75	+ ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
76	) / ( 3.0 * gg ) - u_gtrans
77	IF ( k+1 == nzt+1 ) THEN
78	u_int = u_int_l
79	ELSE
80	u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) &
81	+ ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
82	) / ( 3.0 * gg ) - u_gtrans
83	u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
84	( u_int_u - u_int_l )
85	ENDIF
86
87	!
88	!-- Same procedure for interpolation of the v velocity-component (adopt
89	!-- index k from u velocity-component)
90	i = particles(n)%x * ddx
91	j = ( particles(n)%y + 0.5 * dy ) * ddy
92
93	x = particles(n)%x - i * dx
94	y = particles(n)%y + ( 0.5 - j ) * dy
95	aa = x2 + y2
96	bb = ( dx - x )2 + y2
97	cc = x2 + ( dy - y )2
98	dd = ( dx - x )2 + ( dy - y )2
99	gg = aa + bb + cc + dd
100
101	v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) &
102	+ ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
103	) / ( 3.0 * gg ) - v_gtrans
104	IF ( k+1 == nzt+1 ) THEN
105	v_int = v_int_l
106	ELSE
107	v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) &
108	+ ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
109	) / ( 3.0 * gg ) - v_gtrans
110	v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
111	( v_int_u - v_int_l )
112	ENDIF
113
114	!
115	!-- Same procedure for interpolation of the w velocity-component (adopt
116	!-- index i from v velocity-component)
117	IF ( vertical_particle_advection(particles(n)%group) ) THEN
118	j = particles(n)%y * ddy
119	k = particles(n)%z / dz + offset_ocean_nzt_m1
120
121	x = particles(n)%x - i * dx
122	y = particles(n)%y - j * dy
123	aa = x2 + y2
124	bb = ( dx - x )2 + y2
125	cc = x2 + ( dy - y )2
126	dd = ( dx - x )2 + ( dy - y )2
127	gg = aa + bb + cc + dd
128
129	w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) &
130	+ ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
131	) / ( 3.0 * gg )
132	IF ( k+1 == nzt+1 ) THEN
133	w_int = w_int_l
134	ELSE
135	w_int_u = ( ( gg-aa ) * w(k+1,j,i) + &
136	( gg-bb ) * w(k+1,j,i+1) + &
137	( gg-cc ) * w(k+1,j+1,i) + &
138	( gg-dd ) * w(k+1,j+1,i+1) &
139	) / ( 3.0 * gg )
140	w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
141	( w_int_u - w_int_l )
142	ENDIF
143	ELSE
144	w_int = 0.0
145	ENDIF
146
147	!
148	!-- Interpolate and calculate quantities needed for calculating the SGS
149	!-- velocities
150	IF ( use_sgs_for_particles ) THEN
151	!
152	!-- Interpolate TKE
153	i = particles(n)%x * ddx
154	j = particles(n)%y * ddy
155	k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz &
156	+ offset_ocean_nzt ! only exact if eq.dist
157
158	IF ( topography == 'flat' ) THEN
159
160	x = particles(n)%x - i * dx
161	y = particles(n)%y - j * dy
162	aa = x2 + y2
163	bb = ( dx - x )2 + y2
164	cc = x2 + ( dy - y )2
165	dd = ( dx - x )2 + ( dy - y )2
166	gg = aa + bb + cc + dd
167
168	e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
169	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
170	) / ( 3.0 * gg )
171
172	IF ( k+1 == nzt+1 ) THEN
173	e_int = e_int_l
174	ELSE
175	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
176	( gg - bb ) * e(k+1,j,i+1) + &
177	( gg - cc ) * e(k+1,j+1,i) + &
178	( gg - dd ) * e(k+1,j+1,i+1) &
179	) / ( 3.0 * gg )
180	e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
181	( e_int_u - e_int_l )
182	ENDIF
183
184	!
185	!-- Interpolate the TKE gradient along x (adopt incides i,j,k and
186	!-- all position variables from above (TKE))
187	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
188	( gg - bb ) * de_dx(k,j,i+1) + &
189	( gg - cc ) * de_dx(k,j+1,i) + &
190	( gg - dd ) * de_dx(k,j+1,i+1) &
191	) / ( 3.0 * gg )
192
193	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
194	de_dx_int = de_dx_int_l
195	ELSE
196	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
197	( gg - bb ) * de_dx(k+1,j,i+1) + &
198	( gg - cc ) * de_dx(k+1,j+1,i) + &
199	( gg - dd ) * de_dx(k+1,j+1,i+1) &
200	) / ( 3.0 * gg )
201	de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
202	( de_dx_int_u - de_dx_int_l )
203	ENDIF
204
205	!
206	!-- Interpolate the TKE gradient along y
207	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
208	( gg - bb ) * de_dy(k,j,i+1) + &
209	( gg - cc ) * de_dy(k,j+1,i) + &
210	( gg - dd ) * de_dy(k,j+1,i+1) &
211	) / ( 3.0 * gg )
212	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
213	de_dy_int = de_dy_int_l
214	ELSE
215	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
216	( gg - bb ) * de_dy(k+1,j,i+1) + &
217	( gg - cc ) * de_dy(k+1,j+1,i) + &
218	( gg - dd ) * de_dy(k+1,j+1,i+1) &
219	) / ( 3.0 * gg )
220	de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
221	( de_dy_int_u - de_dy_int_l )
222	ENDIF
223
224	!
225	!-- Interpolate the TKE gradient along z
226	IF ( particles(n)%z < 0.5 * dz ) THEN
227	de_dz_int = 0.0
228	ELSE
229	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
230	( gg - bb ) * de_dz(k,j,i+1) + &
231	( gg - cc ) * de_dz(k,j+1,i) + &
232	( gg - dd ) * de_dz(k,j+1,i+1) &
233	) / ( 3.0 * gg )
234
235	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
236	de_dz_int = de_dz_int_l
237	ELSE
238	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
239	( gg - bb ) * de_dz(k+1,j,i+1) + &
240	( gg - cc ) * de_dz(k+1,j+1,i) + &
241	( gg - dd ) * de_dz(k+1,j+1,i+1) &
242	) / ( 3.0 * gg )
243	de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
244	( de_dz_int_u - de_dz_int_l )
245	ENDIF
246	ENDIF
247
248	!
249	!-- Interpolate the dissipation of TKE
250	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
251	( gg - bb ) * diss(k,j,i+1) + &
252	( gg - cc ) * diss(k,j+1,i) + &
253	( gg - dd ) * diss(k,j+1,i+1) &
254	) / ( 3.0 * gg )
255
256	IF ( k+1 == nzt+1 ) THEN
257	diss_int = diss_int_l
258	ELSE
259	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
260	( gg - bb ) * diss(k+1,j,i+1) + &
261	( gg - cc ) * diss(k+1,j+1,i) + &
262	( gg - dd ) * diss(k+1,j+1,i+1) &
263	) / ( 3.0 * gg )
264	diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
265	( diss_int_u - diss_int_l )
266	ENDIF
267
268	ELSE
269
270	!
271	!-- In case that there are buildings it has to be determined
272	!-- how many of the gridpoints defining the particle box are
273	!-- situated within a building
274	!-- gp_outside_of_building(1): i,j,k
275	!-- gp_outside_of_building(2): i,j+1,k
276	!-- gp_outside_of_building(3): i,j,k+1
277	!-- gp_outside_of_building(4): i,j+1,k+1
278	!-- gp_outside_of_building(5): i+1,j,k
279	!-- gp_outside_of_building(6): i+1,j+1,k
280	!-- gp_outside_of_building(7): i+1,j,k+1
281	!-- gp_outside_of_building(8): i+1,j+1,k+1
282
283	gp_outside_of_building = 0
284	location = 0.0
285	num_gp = 0
286
287	IF ( k > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
288	num_gp = num_gp + 1
289	gp_outside_of_building(1) = 1
290	location(num_gp,1) = i * dx
291	location(num_gp,2) = j * dy
292	location(num_gp,3) = k * dz - 0.5 * dz
293	ei(num_gp) = e(k,j,i)
294	dissi(num_gp) = diss(k,j,i)
295	de_dxi(num_gp) = de_dx(k,j,i)
296	de_dyi(num_gp) = de_dy(k,j,i)
297	de_dzi(num_gp) = de_dz(k,j,i)
298	ENDIF
299
300	IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) &
301	THEN
302	num_gp = num_gp + 1
303	gp_outside_of_building(2) = 1
304	location(num_gp,1) = i * dx
305	location(num_gp,2) = (j+1) * dy
306	location(num_gp,3) = k * dz - 0.5 * dz
307	ei(num_gp) = e(k,j+1,i)
308	dissi(num_gp) = diss(k,j+1,i)
309	de_dxi(num_gp) = de_dx(k,j+1,i)
310	de_dyi(num_gp) = de_dy(k,j+1,i)
311	de_dzi(num_gp) = de_dz(k,j+1,i)
312	ENDIF
313
314	IF ( k+1 > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
315	num_gp = num_gp + 1
316	gp_outside_of_building(3) = 1
317	location(num_gp,1) = i * dx
318	location(num_gp,2) = j * dy
319	location(num_gp,3) = (k+1) * dz - 0.5 * dz
320	ei(num_gp) = e(k+1,j,i)
321	dissi(num_gp) = diss(k+1,j,i)
322	de_dxi(num_gp) = de_dx(k+1,j,i)
323	de_dyi(num_gp) = de_dy(k+1,j,i)
324	de_dzi(num_gp) = de_dz(k+1,j,i)
325	ENDIF
326
327	IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) &
328	THEN
329	num_gp = num_gp + 1
330	gp_outside_of_building(4) = 1
331	location(num_gp,1) = i * dx
332	location(num_gp,2) = (j+1) * dy
333	location(num_gp,3) = (k+1) * dz - 0.5 * dz
334	ei(num_gp) = e(k+1,j+1,i)
335	dissi(num_gp) = diss(k+1,j+1,i)
336	de_dxi(num_gp) = de_dx(k+1,j+1,i)
337	de_dyi(num_gp) = de_dy(k+1,j+1,i)
338	de_dzi(num_gp) = de_dz(k+1,j+1,i)
339	ENDIF
340
341	IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) &
342	THEN
343	num_gp = num_gp + 1
344	gp_outside_of_building(5) = 1
345	location(num_gp,1) = (i+1) * dx
346	location(num_gp,2) = j * dy
347	location(num_gp,3) = k * dz - 0.5 * dz
348	ei(num_gp) = e(k,j,i+1)
349	dissi(num_gp) = diss(k,j,i+1)
350	de_dxi(num_gp) = de_dx(k,j,i+1)
351	de_dyi(num_gp) = de_dy(k,j,i+1)
352	de_dzi(num_gp) = de_dz(k,j,i+1)
353	ENDIF
354
355	IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) &
356	THEN
357	num_gp = num_gp + 1
358	gp_outside_of_building(6) = 1
359	location(num_gp,1) = (i+1) * dx
360	location(num_gp,2) = (j+1) * dy
361	location(num_gp,3) = k * dz - 0.5 * dz
362	ei(num_gp) = e(k,j+1,i+1)
363	dissi(num_gp) = diss(k,j+1,i+1)
364	de_dxi(num_gp) = de_dx(k,j+1,i+1)
365	de_dyi(num_gp) = de_dy(k,j+1,i+1)
366	de_dzi(num_gp) = de_dz(k,j+1,i+1)
367	ENDIF
368
369	IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) &
370	THEN
371	num_gp = num_gp + 1
372	gp_outside_of_building(7) = 1
373	location(num_gp,1) = (i+1) * dx
374	location(num_gp,2) = j * dy
375	location(num_gp,3) = (k+1) * dz - 0.5 * dz
376	ei(num_gp) = e(k+1,j,i+1)
377	dissi(num_gp) = diss(k+1,j,i+1)
378	de_dxi(num_gp) = de_dx(k+1,j,i+1)
379	de_dyi(num_gp) = de_dy(k+1,j,i+1)
380	de_dzi(num_gp) = de_dz(k+1,j,i+1)
381	ENDIF
382
383	IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)&
384	THEN
385	num_gp = num_gp + 1
386	gp_outside_of_building(8) = 1
387	location(num_gp,1) = (i+1) * dx
388	location(num_gp,2) = (j+1) * dy
389	location(num_gp,3) = (k+1) * dz - 0.5 * dz
390	ei(num_gp) = e(k+1,j+1,i+1)
391	dissi(num_gp) = diss(k+1,j+1,i+1)
392	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
393	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
394	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
395	ENDIF
396
397	!
398	!-- If all gridpoints are situated outside of a building, then the
399	!-- ordinary interpolation scheme can be used.
400	IF ( num_gp == 8 ) THEN
401
402	x = particles(n)%x - i * dx
403	y = particles(n)%y - j * dy
404	aa = x2 + y2
405	bb = ( dx - x )2 + y2
406	cc = x2 + ( dy - y )2
407	dd = ( dx - x )2 + ( dy - y )2
408	gg = aa + bb + cc + dd
409
410	e_int_l = (( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
411	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)&
412	) / ( 3.0 * gg )
413
414	IF ( k+1 == nzt+1 ) THEN
415	e_int = e_int_l
416	ELSE
417	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
418	( gg - bb ) * e(k+1,j,i+1) + &
419	( gg - cc ) * e(k+1,j+1,i) + &
420	( gg - dd ) * e(k+1,j+1,i+1) &
421	) / ( 3.0 * gg )
422	e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
423	( e_int_u - e_int_l )
424	ENDIF
425
426	!
427	!-- Interpolate the TKE gradient along x (adopt incides i,j,k
428	!-- and all position variables from above (TKE))
429	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
430	( gg - bb ) * de_dx(k,j,i+1) + &
431	( gg - cc ) * de_dx(k,j+1,i) + &
432	( gg - dd ) * de_dx(k,j+1,i+1) &
433	) / ( 3.0 * gg )
434
435	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
436	de_dx_int = de_dx_int_l
437	ELSE
438	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
439	( gg - bb ) * de_dx(k+1,j,i+1) + &
440	( gg - cc ) * de_dx(k+1,j+1,i) + &
441	( gg - dd ) * de_dx(k+1,j+1,i+1) &
442	) / ( 3.0 * gg )
443	de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / &
444	dz * ( de_dx_int_u - de_dx_int_l )
445	ENDIF
446
447	!
448	!-- Interpolate the TKE gradient along y
449	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
450	( gg - bb ) * de_dy(k,j,i+1) + &
451	( gg - cc ) * de_dy(k,j+1,i) + &
452	( gg - dd ) * de_dy(k,j+1,i+1) &
453	) / ( 3.0 * gg )
454	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
455	de_dy_int = de_dy_int_l
456	ELSE
457	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
458	( gg - bb ) * de_dy(k+1,j,i+1) + &
459	( gg - cc ) * de_dy(k+1,j+1,i) + &
460	( gg - dd ) * de_dy(k+1,j+1,i+1) &
461	) / ( 3.0 * gg )
462	de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / &
463	dz * ( de_dy_int_u - de_dy_int_l )
464	ENDIF
465
466	!
467	!-- Interpolate the TKE gradient along z
468	IF ( particles(n)%z < 0.5 * dz ) THEN
469	de_dz_int = 0.0
470	ELSE
471	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
472	( gg - bb ) * de_dz(k,j,i+1) + &
473	( gg - cc ) * de_dz(k,j+1,i) + &
474	( gg - dd ) * de_dz(k,j+1,i+1) &
475	) / ( 3.0 * gg )
476
477	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
478	de_dz_int = de_dz_int_l
479	ELSE
480	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
481	( gg - bb ) * de_dz(k+1,j,i+1) + &
482	( gg - cc ) * de_dz(k+1,j+1,i) + &
483	( gg - dd ) * de_dz(k+1,j+1,i+1) &
484	) / ( 3.0 * gg )
485	de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /&
486	dz * ( de_dz_int_u - de_dz_int_l )
487	ENDIF
488	ENDIF
489
490	!
491	!-- Interpolate the dissipation of TKE
492	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
493	( gg - bb ) * diss(k,j,i+1) + &
494	( gg - cc ) * diss(k,j+1,i) + &
495	( gg - dd ) * diss(k,j+1,i+1) &
496	) / ( 3.0 * gg )
497
498	IF ( k+1 == nzt+1 ) THEN
499	diss_int = diss_int_l
500	ELSE
501	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
502	( gg - bb ) * diss(k+1,j,i+1) + &
503	( gg - cc ) * diss(k+1,j+1,i) + &
504	( gg - dd ) * diss(k+1,j+1,i+1) &
505	) / ( 3.0 * gg )
506	diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *&
507	( diss_int_u - diss_int_l )
508	ENDIF
509
510	ELSE
511
512	!
513	!-- If wall between gridpoint 1 and gridpoint 5, then
514	!-- Neumann boundary condition has to be applied
515	IF ( gp_outside_of_building(1) == 1 .AND. &
516	gp_outside_of_building(5) == 0 ) THEN
517	num_gp = num_gp + 1
518	location(num_gp,1) = i * dx + 0.5 * dx
519	location(num_gp,2) = j * dy
520	location(num_gp,3) = k * dz - 0.5 * dz
521	ei(num_gp) = e(k,j,i)
522	dissi(num_gp) = diss(k,j,i)
523	de_dxi(num_gp) = 0.0
524	de_dyi(num_gp) = de_dy(k,j,i)
525	de_dzi(num_gp) = de_dz(k,j,i)
526	ENDIF
527
528	IF ( gp_outside_of_building(5) == 1 .AND. &
529	gp_outside_of_building(1) == 0 ) THEN
530	num_gp = num_gp + 1
531	location(num_gp,1) = i * dx + 0.5 * dx
532	location(num_gp,2) = j * dy
533	location(num_gp,3) = k * dz - 0.5 * dz
534	ei(num_gp) = e(k,j,i+1)
535	dissi(num_gp) = diss(k,j,i+1)
536	de_dxi(num_gp) = 0.0
537	de_dyi(num_gp) = de_dy(k,j,i+1)
538	de_dzi(num_gp) = de_dz(k,j,i+1)
539	ENDIF
540
541	!
542	!-- If wall between gridpoint 5 and gridpoint 6, then
543	!-- then Neumann boundary condition has to be applied
544	IF ( gp_outside_of_building(5) == 1 .AND. &
545	gp_outside_of_building(6) == 0 ) THEN
546	num_gp = num_gp + 1
547	location(num_gp,1) = (i+1) * dx
548	location(num_gp,2) = j * dy + 0.5 * dy
549	location(num_gp,3) = k * dz - 0.5 * dz
550	ei(num_gp) = e(k,j,i+1)
551	dissi(num_gp) = diss(k,j,i+1)
552	de_dxi(num_gp) = de_dx(k,j,i+1)
553	de_dyi(num_gp) = 0.0
554	de_dzi(num_gp) = de_dz(k,j,i+1)
555	ENDIF
556
557	IF ( gp_outside_of_building(6) == 1 .AND. &
558	gp_outside_of_building(5) == 0 ) THEN
559	num_gp = num_gp + 1
560	location(num_gp,1) = (i+1) * dx
561	location(num_gp,2) = j * dy + 0.5 * dy
562	location(num_gp,3) = k * dz - 0.5 * dz
563	ei(num_gp) = e(k,j+1,i+1)
564	dissi(num_gp) = diss(k,j+1,i+1)
565	de_dxi(num_gp) = de_dx(k,j+1,i+1)
566	de_dyi(num_gp) = 0.0
567	de_dzi(num_gp) = de_dz(k,j+1,i+1)
568	ENDIF
569
570	!
571	!-- If wall between gridpoint 2 and gridpoint 6, then
572	!-- Neumann boundary condition has to be applied
573	IF ( gp_outside_of_building(2) == 1 .AND. &
574	gp_outside_of_building(6) == 0 ) THEN
575	num_gp = num_gp + 1
576	location(num_gp,1) = i * dx + 0.5 * dx
577	location(num_gp,2) = (j+1) * dy
578	location(num_gp,3) = k * dz - 0.5 * dz
579	ei(num_gp) = e(k,j+1,i)
580	dissi(num_gp) = diss(k,j+1,i)
581	de_dxi(num_gp) = 0.0
582	de_dyi(num_gp) = de_dy(k,j+1,i)
583	de_dzi(num_gp) = de_dz(k,j+1,i)
584	ENDIF
585
586	IF ( gp_outside_of_building(6) == 1 .AND. &
587	gp_outside_of_building(2) == 0 ) THEN
588	num_gp = num_gp + 1
589	location(num_gp,1) = i * dx + 0.5 * dx
590	location(num_gp,2) = (j+1) * dy
591	location(num_gp,3) = k * dz - 0.5 * dz
592	ei(num_gp) = e(k,j+1,i+1)
593	dissi(num_gp) = diss(k,j+1,i+1)
594	de_dxi(num_gp) = 0.0
595	de_dyi(num_gp) = de_dy(k,j+1,i+1)
596	de_dzi(num_gp) = de_dz(k,j+1,i+1)
597	ENDIF
598
599	!
600	!-- If wall between gridpoint 1 and gridpoint 2, then
601	!-- Neumann boundary condition has to be applied
602	IF ( gp_outside_of_building(1) == 1 .AND. &
603	gp_outside_of_building(2) == 0 ) THEN
604	num_gp = num_gp + 1
605	location(num_gp,1) = i * dx
606	location(num_gp,2) = j * dy + 0.5 * dy
607	location(num_gp,3) = k * dz - 0.5 * dz
608	ei(num_gp) = e(k,j,i)
609	dissi(num_gp) = diss(k,j,i)
610	de_dxi(num_gp) = de_dx(k,j,i)
611	de_dyi(num_gp) = 0.0
612	de_dzi(num_gp) = de_dz(k,j,i)
613	ENDIF
614
615	IF ( gp_outside_of_building(2) == 1 .AND. &
616	gp_outside_of_building(1) == 0 ) THEN
617	num_gp = num_gp + 1
618	location(num_gp,1) = i * dx
619	location(num_gp,2) = j * dy + 0.5 * dy
620	location(num_gp,3) = k * dz - 0.5 * dz
621	ei(num_gp) = e(k,j+1,i)
622	dissi(num_gp) = diss(k,j+1,i)
623	de_dxi(num_gp) = de_dx(k,j+1,i)
624	de_dyi(num_gp) = 0.0
625	de_dzi(num_gp) = de_dz(k,j+1,i)
626	ENDIF
627
628	!
629	!-- If wall between gridpoint 3 and gridpoint 7, then
630	!-- Neumann boundary condition has to be applied
631	IF ( gp_outside_of_building(3) == 1 .AND. &
632	gp_outside_of_building(7) == 0 ) THEN
633	num_gp = num_gp + 1
634	location(num_gp,1) = i * dx + 0.5 * dx
635	location(num_gp,2) = j * dy
636	location(num_gp,3) = k * dz + 0.5 * dz
637	ei(num_gp) = e(k+1,j,i)
638	dissi(num_gp) = diss(k+1,j,i)
639	de_dxi(num_gp) = 0.0
640	de_dyi(num_gp) = de_dy(k+1,j,i)
641	de_dzi(num_gp) = de_dz(k+1,j,i)
642	ENDIF
643
644	IF ( gp_outside_of_building(7) == 1 .AND. &
645	gp_outside_of_building(3) == 0 ) THEN
646	num_gp = num_gp + 1
647	location(num_gp,1) = i * dx + 0.5 * dx
648	location(num_gp,2) = j * dy
649	location(num_gp,3) = k * dz + 0.5 * dz
650	ei(num_gp) = e(k+1,j,i+1)
651	dissi(num_gp) = diss(k+1,j,i+1)
652	de_dxi(num_gp) = 0.0
653	de_dyi(num_gp) = de_dy(k+1,j,i+1)
654	de_dzi(num_gp) = de_dz(k+1,j,i+1)
655	ENDIF
656
657	!
658	!-- If wall between gridpoint 7 and gridpoint 8, then
659	!-- Neumann boundary condition has to be applied
660	IF ( gp_outside_of_building(7) == 1 .AND. &
661	gp_outside_of_building(8) == 0 ) THEN
662	num_gp = num_gp + 1
663	location(num_gp,1) = (i+1) * dx
664	location(num_gp,2) = j * dy + 0.5 * dy
665	location(num_gp,3) = k * dz + 0.5 * dz
666	ei(num_gp) = e(k+1,j,i+1)
667	dissi(num_gp) = diss(k+1,j,i+1)
668	de_dxi(num_gp) = de_dx(k+1,j,i+1)
669	de_dyi(num_gp) = 0.0
670	de_dzi(num_gp) = de_dz(k+1,j,i+1)
671	ENDIF
672
673	IF ( gp_outside_of_building(8) == 1 .AND. &
674	gp_outside_of_building(7) == 0 ) THEN
675	num_gp = num_gp + 1
676	location(num_gp,1) = (i+1) * dx
677	location(num_gp,2) = j * dy + 0.5 * dy
678	location(num_gp,3) = k * dz + 0.5 * dz
679	ei(num_gp) = e(k+1,j+1,i+1)
680	dissi(num_gp) = diss(k+1,j+1,i+1)
681	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
682	de_dyi(num_gp) = 0.0
683	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
684	ENDIF
685
686	!
687	!-- If wall between gridpoint 4 and gridpoint 8, then
688	!-- Neumann boundary condition has to be applied
689	IF ( gp_outside_of_building(4) == 1 .AND. &
690	gp_outside_of_building(8) == 0 ) THEN
691	num_gp = num_gp + 1
692	location(num_gp,1) = i * dx + 0.5 * dx
693	location(num_gp,2) = (j+1) * dy
694	location(num_gp,3) = k * dz + 0.5 * dz
695	ei(num_gp) = e(k+1,j+1,i)
696	dissi(num_gp) = diss(k+1,j+1,i)
697	de_dxi(num_gp) = 0.0
698	de_dyi(num_gp) = de_dy(k+1,j+1,i)
699	de_dzi(num_gp) = de_dz(k+1,j+1,i)
700	ENDIF
701
702	IF ( gp_outside_of_building(8) == 1 .AND. &
703	gp_outside_of_building(4) == 0 ) THEN
704	num_gp = num_gp + 1
705	location(num_gp,1) = i * dx + 0.5 * dx
706	location(num_gp,2) = (j+1) * dy
707	location(num_gp,3) = k * dz + 0.5 * dz
708	ei(num_gp) = e(k+1,j+1,i+1)
709	dissi(num_gp) = diss(k+1,j+1,i+1)
710	de_dxi(num_gp) = 0.0
711	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
712	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
713	ENDIF
714
715	!
716	!-- If wall between gridpoint 3 and gridpoint 4, then
717	!-- Neumann boundary condition has to be applied
718	IF ( gp_outside_of_building(3) == 1 .AND. &
719	gp_outside_of_building(4) == 0 ) THEN
720	num_gp = num_gp + 1
721	location(num_gp,1) = i * dx
722	location(num_gp,2) = j * dy + 0.5 * dy
723	location(num_gp,3) = k * dz + 0.5 * dz
724	ei(num_gp) = e(k+1,j,i)
725	dissi(num_gp) = diss(k+1,j,i)
726	de_dxi(num_gp) = de_dx(k+1,j,i)
727	de_dyi(num_gp) = 0.0
728	de_dzi(num_gp) = de_dz(k+1,j,i)
729	ENDIF
730
731	IF ( gp_outside_of_building(4) == 1 .AND. &
732	gp_outside_of_building(3) == 0 ) THEN
733	num_gp = num_gp + 1
734	location(num_gp,1) = i * dx
735	location(num_gp,2) = j * dy + 0.5 * dy
736	location(num_gp,3) = k * dz + 0.5 * dz
737	ei(num_gp) = e(k+1,j+1,i)
738	dissi(num_gp) = diss(k+1,j+1,i)
739	de_dxi(num_gp) = de_dx(k+1,j+1,i)
740	de_dyi(num_gp) = 0.0
741	de_dzi(num_gp) = de_dz(k+1,j+1,i)
742	ENDIF
743
744	!
745	!-- If wall between gridpoint 1 and gridpoint 3, then
746	!-- Neumann boundary condition has to be applied
747	!-- (only one case as only building beneath is possible)
748	IF ( gp_outside_of_building(1) == 0 .AND. &
749	gp_outside_of_building(3) == 1 ) THEN
750	num_gp = num_gp + 1
751	location(num_gp,1) = i * dx
752	location(num_gp,2) = j * dy
753	location(num_gp,3) = k * dz
754	ei(num_gp) = e(k+1,j,i)
755	dissi(num_gp) = diss(k+1,j,i)
756	de_dxi(num_gp) = de_dx(k+1,j,i)
757	de_dyi(num_gp) = de_dy(k+1,j,i)
758	de_dzi(num_gp) = 0.0
759	ENDIF
760
761	!
762	!-- If wall between gridpoint 5 and gridpoint 7, then
763	!-- Neumann boundary condition has to be applied
764	!-- (only one case as only building beneath is possible)
765	IF ( gp_outside_of_building(5) == 0 .AND. &
766	gp_outside_of_building(7) == 1 ) THEN
767	num_gp = num_gp + 1
768	location(num_gp,1) = (i+1) * dx
769	location(num_gp,2) = j * dy
770	location(num_gp,3) = k * dz
771	ei(num_gp) = e(k+1,j,i+1)
772	dissi(num_gp) = diss(k+1,j,i+1)
773	de_dxi(num_gp) = de_dx(k+1,j,i+1)
774	de_dyi(num_gp) = de_dy(k+1,j,i+1)
775	de_dzi(num_gp) = 0.0
776	ENDIF
777
778	!
779	!-- If wall between gridpoint 2 and gridpoint 4, then
780	!-- Neumann boundary condition has to be applied
781	!-- (only one case as only building beneath is possible)
782	IF ( gp_outside_of_building(2) == 0 .AND. &
783	gp_outside_of_building(4) == 1 ) THEN
784	num_gp = num_gp + 1
785	location(num_gp,1) = i * dx
786	location(num_gp,2) = (j+1) * dy
787	location(num_gp,3) = k * dz
788	ei(num_gp) = e(k+1,j+1,i)
789	dissi(num_gp) = diss(k+1,j+1,i)
790	de_dxi(num_gp) = de_dx(k+1,j+1,i)
791	de_dyi(num_gp) = de_dy(k+1,j+1,i)
792	de_dzi(num_gp) = 0.0
793	ENDIF
794
795	!
796	!-- If wall between gridpoint 6 and gridpoint 8, then
797	!-- Neumann boundary condition has to be applied
798	!-- (only one case as only building beneath is possible)
799	IF ( gp_outside_of_building(6) == 0 .AND. &
800	gp_outside_of_building(8) == 1 ) THEN
801	num_gp = num_gp + 1
802	location(num_gp,1) = (i+1) * dx
803	location(num_gp,2) = (j+1) * dy
804	location(num_gp,3) = k * dz
805	ei(num_gp) = e(k+1,j+1,i+1)
806	dissi(num_gp) = diss(k+1,j+1,i+1)
807	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
808	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
809	de_dzi(num_gp) = 0.0
810	ENDIF
811
812	!
813	!-- Carry out the interpolation
814	IF ( num_gp == 1 ) THEN
815	!
816	!-- If only one of the gridpoints is situated outside of the
817	!-- building, it follows that the values at the particle
818	!-- location are the same as the gridpoint values
819	e_int = ei(num_gp)
820	diss_int = dissi(num_gp)
821	de_dx_int = de_dxi(num_gp)
822	de_dy_int = de_dyi(num_gp)
823	de_dz_int = de_dzi(num_gp)
824	ELSE IF ( num_gp > 1 ) THEN
825
826	d_sum = 0.0
827	!
828	!-- Evaluation of the distances between the gridpoints
829	!-- contributing to the interpolated values, and the particle
830	!-- location
831	DO agp = 1, num_gp
832	d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 &
833	+ ( particles(n)%y-location(agp,2) )**2 &
834	+ ( particles(n)%z-location(agp,3) )**2
835	d_sum = d_sum + d_gp_pl(agp)
836	ENDDO
837
838	!
839	!-- Finally the interpolation can be carried out
840	e_int = 0.0
841	diss_int = 0.0
842	de_dx_int = 0.0
843	de_dy_int = 0.0
844	de_dz_int = 0.0
845	DO agp = 1, num_gp
846	e_int = e_int + ( d_sum - d_gp_pl(agp) ) * &
847	ei(agp) / ( (num_gp-1) * d_sum )
848	diss_int = diss_int + ( d_sum - d_gp_pl(agp) ) * &
849	dissi(agp) / ( (num_gp-1) * d_sum )
850	de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * &
851	de_dxi(agp) / ( (num_gp-1) * d_sum )
852	de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * &
853	de_dyi(agp) / ( (num_gp-1) * d_sum )
854	de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * &
855	de_dzi(agp) / ( (num_gp-1) * d_sum )
856	ENDDO
857
858	ENDIF
859
860	ENDIF
861
862	ENDIF
863
864	!
865	!-- Vertically interpolate the horizontally averaged SGS TKE and
866	!-- resolved-scale velocity variances and use the interpolated values
867	!-- to calculate the coefficient fs, which is a measure of the ratio
868	!-- of the subgrid-scale turbulent kinetic energy to the total amount
869	!-- of turbulent kinetic energy.
870	IF ( k == 0 ) THEN
871	e_mean_int = hom(0,1,8,0)
872	ELSE
873	e_mean_int = hom(k,1,8,0) + &
874	( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
875	( zu(k+1) - zu(k) ) * &
876	( particles(n)%z - zu(k) )
877	ENDIF
878
879	kw = particles(n)%z / dz
880
881	IF ( k == 0 ) THEN
882	aa = hom(k+1,1,30,0) * ( particles(n)%z / &
883	( 0.5 * ( zu(k+1) - zu(k) ) ) )
884	bb = hom(k+1,1,31,0) * ( particles(n)%z / &
885	( 0.5 * ( zu(k+1) - zu(k) ) ) )
886	cc = hom(kw+1,1,32,0) * ( particles(n)%z / &
887	( 1.0 * ( zw(kw+1) - zw(kw) ) ) )
888	ELSE
889	aa = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
890	( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
891	bb = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
892	( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
893	cc = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *&
894	( ( particles(n)%z - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
895	ENDIF
896
897	vv_int = ( 1.0 / 3.0 ) * ( aa + bb + cc )
898
899	fs_int = ( 2.0 / 3.0 ) * e_mean_int / &
900	( vv_int + ( 2.0 / 3.0 ) * e_mean_int )
901
902	!
903	!-- Calculate the Lagrangian timescale according to Weil et al. (2004).
904	lagr_timescale = ( 4.0 * e_int ) / &
905	( 3.0 * fs_int * c_0 * diss_int )
906
907	!
908	!-- Calculate the next particle timestep. dt_gap is the time needed to
909	!-- complete the current LES timestep.
910	dt_gap = dt_3d - particles(n)%dt_sum
911	dt_particle = MIN( dt_3d, 0.025 * lagr_timescale, dt_gap )
912
913	!
914	!-- The particle timestep should not be too small in order to prevent
915	!-- the number of particle timesteps of getting too large
916	IF ( dt_particle < dt_min_part .AND. dt_min_part < dt_gap ) &
917	THEN
918	dt_particle = dt_min_part
919	ENDIF
920
921	!
922	!-- Calculate the SGS velocity components
923	IF ( particles(n)%age == 0.0 ) THEN
924	!
925	!-- For new particles the SGS components are derived from the SGS
926	!-- TKE. Limit the Gaussian random number to the interval
927	!-- [-5.0sigma, 5.0sigma] in order to prevent the SGS velocities
928	!-- from becoming unrealistically large.
929	particles(n)%rvar1 = SQRT( 2.0 * sgs_wfu_part * e_int ) * &
930	( random_gauss( iran_part, 5.0 ) - 1.0 )
931	particles(n)%rvar2 = SQRT( 2.0 * sgs_wfv_part * e_int ) * &
932	( random_gauss( iran_part, 5.0 ) - 1.0 )
933	particles(n)%rvar3 = SQRT( 2.0 * sgs_wfw_part * e_int ) * &
934	( random_gauss( iran_part, 5.0 ) - 1.0 )
935
936	ELSE
937
938	!
939	!-- Restriction of the size of the new timestep: compared to the
940	!-- previous timestep the increase must not exceed 200%
941
942	dt_particle_m = particles(n)%age - particles(n)%age_m
943	IF ( dt_particle > 2.0 * dt_particle_m ) THEN
944	dt_particle = 2.0 * dt_particle_m
945	ENDIF
946
947	!
948	!-- For old particles the SGS components are correlated with the
949	!-- values from the previous timestep. Random numbers have also to
950	!-- be limited (see above).
951	!-- As negative values for the subgrid TKE are not allowed, the
952	!-- change of the subgrid TKE with time cannot be smaller than
953	!-- -e_int/dt_particle. This value is used as a lower boundary
954	!-- value for the change of TKE
955
956	de_dt_min = - e_int / dt_particle
957
958	de_dt = ( e_int - particles(n)%e_m ) / dt_particle_m
959
960	IF ( de_dt < de_dt_min ) THEN
961	de_dt = de_dt_min
962	ENDIF
963
964	particles(n)%rvar1 = particles(n)%rvar1 - fs_int * c_0 * &
965	diss_int * particles(n)%rvar1 * dt_particle /&
966	( 4.0 * sgs_wfu_part * e_int ) + &
967	( 2.0 * sgs_wfu_part * de_dt * &
968	particles(n)%rvar1 / &
969	( 2.0 * sgs_wfu_part * e_int ) + de_dx_int &
970	) * dt_particle / 2.0 + &
971	SQRT( fs_int * c_0 * diss_int ) * &
972	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
973	SQRT( dt_particle )
974
975	particles(n)%rvar2 = particles(n)%rvar2 - fs_int * c_0 * &
976	diss_int * particles(n)%rvar2 * dt_particle /&
977	( 4.0 * sgs_wfv_part * e_int ) + &
978	( 2.0 * sgs_wfv_part * de_dt * &
979	particles(n)%rvar2 / &
980	( 2.0 * sgs_wfv_part * e_int ) + de_dy_int &
981	) * dt_particle / 2.0 + &
982	SQRT( fs_int * c_0 * diss_int ) * &
983	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
984	SQRT( dt_particle )
985
986	particles(n)%rvar3 = particles(n)%rvar3 - fs_int * c_0 * &
987	diss_int * particles(n)%rvar3 * dt_particle /&
988	( 4.0 * sgs_wfw_part * e_int ) + &
989	( 2.0 * sgs_wfw_part * de_dt * &
990	particles(n)%rvar3 / &
991	( 2.0 * sgs_wfw_part * e_int ) + de_dz_int &
992	) * dt_particle / 2.0 + &
993	SQRT( fs_int * c_0 * diss_int ) * &
994	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
995	SQRT( dt_particle )
996
997	ENDIF
998
999	u_int = u_int + particles(n)%rvar1
1000	v_int = v_int + particles(n)%rvar2
1001	w_int = w_int + particles(n)%rvar3
1002
1003	!
1004	!-- Store the SGS TKE of the current timelevel which is needed for
1005	!-- for calculating the SGS particle velocities at the next timestep
1006	particles(n)%e_m = e_int
1007
1008	ELSE
1009	!
1010	!-- If no SGS velocities are used, only the particle timestep has to
1011	!-- be set
1012	dt_particle = dt_3d
1013
1014	ENDIF
1015
1016	!
1017	!-- Store the old age of the particle ( needed to prevent that a
1018	!-- particle crosses several PEs during one timestep, and for the
1019	!-- evaluation of the subgrid particle velocity fluctuations )
1020	particles(n)%age_m = particles(n)%age
1021
1022
1023	!
1024	!-- Particle advection
1025	IF ( particle_groups(particles(n)%group)%density_ratio == 0.0 ) THEN
1026	!
1027	!-- Pure passive transport (without particle inertia)
1028	particles(n)%x = particles(n)%x + u_int * dt_particle
1029	particles(n)%y = particles(n)%y + v_int * dt_particle
1030	particles(n)%z = particles(n)%z + w_int * dt_particle
1031
1032	particles(n)%speed_x = u_int
1033	particles(n)%speed_y = v_int
1034	particles(n)%speed_z = w_int
1035
1036	ELSE
1037	!
1038	!-- Transport of particles with inertia
1039	particles(n)%x = particles(n)%x + particles(n)%speed_x * &
1040	dt_particle
1041	particles(n)%y = particles(n)%y + particles(n)%speed_y * &
1042	dt_particle
1043	particles(n)%z = particles(n)%z + particles(n)%speed_z * &
1044	dt_particle
1045
1046	!
1047	!-- Update of the particle velocity
1048	dens_ratio = particle_groups(particles(n)%group)%density_ratio
1049	IF ( cloud_droplets ) THEN
1050	exp_arg = 4.5 * dens_ratio * molecular_viscosity / &
1051	( particles(n)%radius )*2 &
1052	( 1.0 + 0.15 * ( 2.0 * particles(n)%radius * &
1053	SQRT( ( u_int - particles(n)%speed_x )**2 + &
1054	( v_int - particles(n)%speed_y )**2 + &
1055	( w_int - particles(n)%speed_z )**2 ) / &
1056	molecular_viscosity )**0.687 &
1057	)
1058	exp_term = EXP( -exp_arg * dt_particle )
1059	ELSEIF ( use_sgs_for_particles ) THEN
1060	exp_arg = particle_groups(particles(n)%group)%exp_arg
1061	exp_term = EXP( -exp_arg * dt_particle )
1062	ELSE
1063	exp_arg = particle_groups(particles(n)%group)%exp_arg
1064	exp_term = particle_groups(particles(n)%group)%exp_term
1065	ENDIF
1066	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
1067	u_int * ( 1.0 - exp_term )
1068	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
1069	v_int * ( 1.0 - exp_term )
1070	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
1071	( w_int - ( 1.0 - dens_ratio ) * g / exp_arg )&
1072	* ( 1.0 - exp_term )
1073	ENDIF
1074
1075	!
1076	!-- Increment the particle age and the total time that the particle
1077	!-- has advanced within the particle timestep procedure
1078	particles(n)%age = particles(n)%age + dt_particle
1079	particles(n)%dt_sum = particles(n)%dt_sum + dt_particle
1080
1081	!
1082	!-- Check whether there is still a particle that has not yet completed
1083	!-- the total LES timestep
1084	IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8 ) THEN
1085	dt_3d_reached_l = .FALSE.
1086	ENDIF
1087
1088	ENDDO
1089
1090
1091	END SUBROUTINE lpm_advec

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

Context Navigation

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

Download in other formats: