Home

Context Navigation

source: palm/trunk/SOURCE/advec_particles.f90 @ 130

Last change on this file since 130 was 119, checked in by raasch, 17 years ago
small bugfixes in SGS part of adved_particles
Property svn:keywords set to `Id`
File size: 164.4 KB

Line
1	SUBROUTINE advec_particles
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! Sorting of particles is controlled by dt_sort_particles and moved from
7	! the SGS timestep loop after the end of this loop.
8	! Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
9	! numbers along y
10	! Small bugfixes in the SGS part
11	! TEST: PRINT statements on unit 9 (commented out)
12	!
13	! Former revisions:
14	! -----------------
15	! $Id: advec_particles.f90 119 2007-10-17 10:27:13Z letzel $
16	!
17	! 106 2007-08-16 14:30:26Z raasch
18	! remaining variables iran changed to iran_part
19	!
20	! 95 2007-06-02 16:48:38Z raasch
21	! hydro_press renamed hyp
22	!
23	! 75 2007-03-22 09:54:05Z raasch
24	! Particle reflection at vertical walls implemented in new subroutine
25	! particle_boundary_conds,
26	! vertical walls are regarded in the SGS model,
27	! + user_advec_particles, particles-package is now part of the defaut code,
28	! array arguments in sendrecv calls have to refer to first element (1) due to
29	! mpich (mpiI) interface requirements,
30	! 2nd+3rd argument removed from exchange horiz
31	!
32	! 16 2007-02-15 13:16:47Z raasch
33	! Bugfix: wrong if-clause from revision 1.32
34	!
35	! r4 \| raasch \| 2007-02-13 12:33:16 +0100 (Tue, 13 Feb 2007)
36	! RCS Log replace by Id keyword, revision history cleaned up
37	!
38	! Revision 1.32 2007/02/11 12:48:20 raasch
39	! Allways the lower level k is used for interpolation
40	! Bugfix: new particles are released only if end_time_prel > simulated_time
41	! Bugfix: transfer of particles when x < -0.5*dx (0.0 before), etc.,
42	! index i,j used instead of cartesian (x,y) coordinate to check for
43	! transfer because this failed under very rare conditions
44	! Bugfix: calculation of number of particles with same radius as the current
45	! particle (cloud droplet code)
46	!
47	! Revision 1.31 2006/08/17 09:21:01 raasch
48	! Two more compilation errors removed from the last revision
49	!
50	! Revision 1.30 2006/08/17 09:11:17 raasch
51	! Two compilation errors removed from the last revision
52	!
53	! Revision 1.29 2006/08/04 14:05:01 raasch
54	! Subgrid scale velocities are (optionally) included for calculating the
55	! particle advection, new counters trlp_count_sum, etc. for accumulating
56	! the number of particles exchanged between the subdomains during all
57	! sub-timesteps (if sgs velocities are included), +3d-arrays de_dx/y/z,
58	! izuf renamed iran, output of particle time series
59	!
60	! Revision 1.1 1999/11/25 16:16:06 raasch
61	! Initial revision
62	!
63	!
64	! Description:
65	! ------------
66	! Particle advection
67	!------------------------------------------------------------------------------!
68
69	USE arrays_3d
70	USE cloud_parameters
71	USE constants
72	USE control_parameters
73	USE cpulog
74	USE grid_variables
75	USE indices
76	USE interfaces
77	USE netcdf_control
78	USE particle_attributes
79	USE pegrid
80	USE random_function_mod
81	USE statistics
82
83	IMPLICIT NONE
84
85	INTEGER :: agp, deleted_particles, deleted_tails, i, ie, ii, inc, is, j, &
86	jj, js, k, kk, kw, m, n, nc, nn, num_gp, psi, tlength, &
87	trlp_count, trlp_count_sum, trlp_count_recv, &
88	trlp_count_recv_sum, trlpt_count, trlpt_count_recv, &
89	trnp_count, trnp_count_sum, trnp_count_recv, &
90	trnp_count_recv_sum, trnpt_count, trnpt_count_recv, &
91	trrp_count, trrp_count_sum, trrp_count_recv, &
92	trrp_count_recv_sum, trrpt_count, trrpt_count_recv, &
93	trsp_count, trsp_count_sum, trsp_count_recv, &
94	trsp_count_recv_sum, trspt_count, trspt_count_recv, nd
95
96	INTEGER :: gp_outside_of_building(1:8)
97
98	LOGICAL :: dt_3d_reached, dt_3d_reached_l, prt_position
99
100	REAL :: aa, arg, bb, cc, dd, delta_r, dens_ratio, de_dt, de_dt_min, &
101	de_dx_int, de_dx_int_l, de_dx_int_u, de_dy_int, de_dy_int_l, &
102	de_dy_int_u, de_dz_int, de_dz_int_l, de_dz_int_u, diss_int, &
103	diss_int_l, diss_int_u, distance, dt_gap, dt_particle, &
104	dt_particle_m, d_radius, d_sum, e_a, e_int, e_int_l, e_int_u, &
105	e_mean_int, e_s, exp_arg, exp_term, fs_int, gg, &
106	lagr_timescale, mean_r, new_r, p_int, pt_int, pt_int_l, &
107	pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, ql_int_u, &
108	random_gauss, sl_r3, sl_r4, s_r3, s_r4, t_int, u_int, u_int_l, &
109	u_int_u, vv_int, v_int, v_int_l, v_int_u, w_int, w_int_l, &
110	w_int_u, x, y
111
112	REAL, DIMENSION(1:30) :: de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
113
114	REAL :: location(1:30,1:3)
115
116	REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: de_dx, de_dy, de_dz
117
118	REAL, DIMENSION(:,:,:), ALLOCATABLE :: trlpt, trnpt, trrpt, trspt
119
120	TYPE(particle_type) :: tmp_particle
121
122	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: trlp, trnp, trrp, trsp
123
124
125	CALL cpu_log( log_point(25), 'advec_particles', 'start' )
126
127	! IF ( number_of_particles /= number_of_tails ) THEN
128	! WRITE (9,*) '--- advec_particles: #1'
129	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
130	! CALL local_flush( 9 )
131	! ENDIF
132	!
133	!-- Write particle data on file for later analysis.
134	!-- This has to be done here (before particles are advected) in order
135	!-- to allow correct output in case of dt_write_particle_data = dt_prel =
136	!-- particle_maximum_age. Otherwise (if output is done at the end of this
137	!-- subroutine), the relevant particles would have been already deleted.
138	!-- The MOD function allows for changes in the output interval with restart
139	!-- runs.
140	!-- Attention: change version number for unit 85 (in routine check_open)
141	!-- whenever the output format for this unit is changed!
142	time_write_particle_data = time_write_particle_data + dt_3d
143	IF ( time_write_particle_data >= dt_write_particle_data ) THEN
144
145	CALL cpu_log( log_point_s(40), 'advec_part_io', 'start' )
146	CALL check_open( 85 )
147	WRITE ( 85 ) simulated_time, maximum_number_of_particles, &
148	number_of_particles
149	WRITE ( 85 ) particles
150	WRITE ( 85 ) maximum_number_of_tailpoints, maximum_number_of_tails, &
151	number_of_tails
152	WRITE ( 85 ) particle_tail_coordinates
153	CALL close_file( 85 )
154
155	IF ( netcdf_output ) CALL output_particles_netcdf
156
157	time_write_particle_data = MOD( time_write_particle_data, &
158	MAX( dt_write_particle_data, dt_3d ) )
159	CALL cpu_log( log_point_s(40), 'advec_part_io', 'stop' )
160	ENDIF
161
162	!
163	!-- Calculate exponential term used in case of particle inertia for each
164	!-- of the particle groups
165	CALL cpu_log( log_point_s(41), 'advec_part_exp', 'start' )
166	DO m = 1, number_of_particle_groups
167	IF ( particle_groups(m)%density_ratio /= 0.0 ) THEN
168	particle_groups(m)%exp_arg = &
169	4.5 * particle_groups(m)%density_ratio * &
170	molecular_viscosity / ( particle_groups(m)%radius )**2
171	particle_groups(m)%exp_term = EXP( -particle_groups(m)%exp_arg * &
172	dt_3d )
173	ENDIF
174	ENDDO
175	CALL cpu_log( log_point_s(41), 'advec_part_exp', 'stop' )
176
177	! WRITE ( 9, * ) '*** advec_particles: ##0.3'
178	! CALL local_flush( 9 )
179	! nd = 0
180	! DO n = 1, number_of_particles
181	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
182	! ENDDO
183	! IF ( nd /= deleted_particles ) THEN
184	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
185	! CALL local_flush( 9 )
186	! CALL MPI_ABORT( comm2d, 9999, ierr )
187	! ENDIF
188
189	!
190	!-- Particle (droplet) growth by condensation/evaporation and collision
191	IF ( cloud_droplets ) THEN
192
193	!
194	!-- Reset summation arrays
195	ql_c = 0.0; ql_v = 0.0; ql_vp = 0.0
196
197	!
198	!-- Particle growth by condensation/evaporation
199	CALL cpu_log( log_point_s(42), 'advec_part_cond', 'start' )
200	DO n = 1, number_of_particles
201	!
202	!-- Interpolate temperature and humidity.
203	!-- First determine left, south, and bottom index of the arrays.
204	i = particles(n)%x * ddx
205	j = particles(n)%y * ddy
206	k = ( particles(n)%z + 0.5 * dz ) / dz ! only exact if equidistant
207
208	x = particles(n)%x - i * dx
209	y = particles(n)%y - j * dy
210	aa = x2 + y2
211	bb = ( dx - x )2 + y2
212	cc = x2 + ( dy - y )2
213	dd = ( dx - x )2 + ( dy - y )2
214	gg = aa + bb + cc + dd
215
216	pt_int_l = ( ( gg - aa ) * pt(k,j,i) + ( gg - bb ) * pt(k,j,i+1) &
217	+ ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
218	) / ( 3.0 * gg )
219
220	pt_int_u = ( ( gg-aa ) * pt(k+1,j,i) + ( gg-bb ) * pt(k+1,j,i+1) &
221	+ ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
222	) / ( 3.0 * gg )
223
224	pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
225	( pt_int_u - pt_int_l )
226
227	q_int_l = ( ( gg - aa ) * q(k,j,i) + ( gg - bb ) * q(k,j,i+1) &
228	+ ( gg - cc ) * q(k,j+1,i) + ( gg - dd ) * q(k,j+1,i+1) &
229	) / ( 3.0 * gg )
230
231	q_int_u = ( ( gg-aa ) * q(k+1,j,i) + ( gg-bb ) * q(k+1,j,i+1) &
232	+ ( gg-cc ) * q(k+1,j+1,i) + ( gg-dd ) * q(k+1,j+1,i+1) &
233	) / ( 3.0 * gg )
234
235	q_int = q_int_l + ( particles(n)%z - zu(k) ) / dz * &
236	( q_int_u - q_int_l )
237
238	ql_int_l = ( ( gg - aa ) * ql(k,j,i) + ( gg - bb ) * ql(k,j,i+1) &
239	+ ( gg - cc ) * ql(k,j+1,i) + ( gg - dd ) * ql(k,j+1,i+1) &
240	) / ( 3.0 * gg )
241
242	ql_int_u = ( ( gg-aa ) * ql(k+1,j,i) + ( gg-bb ) * ql(k+1,j,i+1) &
243	+ ( gg-cc ) * ql(k+1,j+1,i) + ( gg-dd ) * ql(k+1,j+1,i+1) &
244	) / ( 3.0 * gg )
245
246	ql_int = ql_int_l + ( particles(n)%z - zu(k) ) / dz * &
247	( ql_int_u - ql_int_l )
248
249	!
250	!-- Calculate real temperature and saturation vapor pressure
251	p_int = hyp(k) + ( particles(n)%z - zu(k) ) / dz * ( hyp(k+1)-hyp(k) )
252	t_int = pt_int * ( p_int / 100000.0 )**0.286
253
254	e_s = 611.0 * EXP( l_d_rv * ( 3.6609E-3 - 1.0 / t_int ) )
255
256	!
257	!-- Current vapor pressure
258	e_a = q_int * p_int / ( 0.378 * q_int + 0.622 )
259
260	!
261	!-- Change in radius by condensation/evaporation
262	!-- ATTENTION: this is only an approximation for large radii
263	arg = particles(n)%radius*2 + 2.0 dt_3d * &
264	( e_a / e_s - 1.0 ) / &
265	( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &
266	thermal_conductivity_l + &
267	rho_l * r_v * t_int / diff_coeff_l / e_s )
268	IF ( arg < 1.0E-14 ) THEN
269	new_r = 1.0E-7
270	ELSE
271	new_r = SQRT( arg )
272	ENDIF
273
274	delta_r = new_r - particles(n)%radius
275
276	! NOTE: this is the correct formula (indipendent of radius).
277	! nevertheless, it give wrong results for large timesteps
278	! d_radius = 1.0 / particles(n)%radius
279	! delta_r = d_radius * ( e_a / e_s - 1.0 - 3.3E-7 / t_int * d_radius + &
280	! b_cond * d_radius**3 ) / &
281	! ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &
282	! thermal_conductivity_l + &
283	! rho_l * r_v * t_int / diff_coeff_l / e_s ) * dt_3d
284
285	! new_r = particles(n)%radius + delta_r
286	! IF ( new_r < 1.0E-7 ) new_r = 1.0E-7
287
288	!
289	!-- Sum up the change in volume of liquid water for the respective grid
290	!-- volume (this is needed later on for calculating the release of
291	!-- latent heat)
292	i = ( particles(n)%x + 0.5 * dx ) * ddx
293	j = ( particles(n)%y + 0.5 * dy ) * ddy
294	k = particles(n)%z / dz + 1 ! only exact if equidistant
295
296	ql_c(k,j,i) = ql_c(k,j,i) + particles(n)%weight_factor * &
297	rho_l * 1.33333333 * pi * &
298	( new_r3 - particles(n)%radius3 ) / &
299	( rho_surface * dx * dy * dz )
300	IF ( ql_c(k,j,i) > 100.0 ) THEN
301	print*,'+++ advec_particles k=',k,' j=',j,' i=',i, &
302	' ql_c=',ql_c(k,j,i), ' part(',n,')%wf=', &
303	particles(n)%weight_factor,' delta_r=',delta_r
304	#if defined( __parallel )
305	CALL MPI_ABORT( comm2d, 9999, ierr )
306	#else
307	STOP
308	#endif
309	ENDIF
310
311	!
312	!-- Change the droplet radius
313	IF ( ( new_r - particles(n)%radius ) < 0.0 .AND. new_r < 0.0 ) &
314	THEN
315	print*,'+++ advec_particles #1 k=',k,' j=',j,' i=',i, &
316	' e_s=',e_s, ' e_a=',e_a,' t_int=',t_int, &
317	' d_radius=',d_radius,' delta_r=',delta_r,&
318	' particle_radius=',particles(n)%radius
319	#if defined( __parallel )
320	CALL MPI_ABORT( comm2d, 9999, ierr )
321	#else
322	STOP
323	#endif
324	ENDIF
325	particles(n)%radius = new_r
326
327	!
328	!-- Sum up the total volume of liquid water (needed below for
329	!-- re-calculating the weighting factors)
330	ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * &
331	particles(n)%radius**3
332	ENDDO
333	CALL cpu_log( log_point_s(42), 'advec_part_cond', 'stop' )
334
335	!
336	!-- Particle growth by collision
337	CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' )
338
339	DO i = nxl, nxr
340	DO j = nys, nyn
341	DO k = nzb+1, nzt
342	!
343	!-- Collision requires at least two particles in the box
344	IF ( prt_count(k,j,i) > 1 ) THEN
345	!
346	!-- First, sort particles within the gridbox by their size,
347	!-- using Shell's method (see Numerical Recipes)
348	!-- NOTE: In case of using particle tails, the re-sorting of
349	!-- ---- tails would have to be included here!
350	psi = prt_start_index(k,j,i) - 1
351	inc = 1
352	DO WHILE ( inc <= prt_count(k,j,i) )
353	inc = 3 * inc + 1
354	ENDDO
355
356	DO WHILE ( inc > 1 )
357	inc = inc / 3
358	DO is = inc+1, prt_count(k,j,i)
359	tmp_particle = particles(psi+is)
360	js = is
361	DO WHILE ( particles(psi+js-inc)%radius > &
362	tmp_particle%radius )
363	particles(psi+js) = particles(psi+js-inc)
364	js = js - inc
365	IF ( js <= inc ) EXIT
366	ENDDO
367	particles(psi+js) = tmp_particle
368	ENDDO
369	ENDDO
370
371	!
372	!-- Calculate the mean radius of all those particles which
373	!-- are of smaller or equal size than the current particle
374	!-- and use this radius for calculating the collision efficiency
375	psi = prt_start_index(k,j,i)
376	s_r3 = 0.0
377	s_r4 = 0.0
378	DO n = psi, psi+prt_count(k,j,i)-1
379	!
380	!-- There may be some particles of size equal to the
381	!-- current particle but with larger index
382	sl_r3 = 0.0
383	sl_r4 = 0.0
384	DO is = n, psi+prt_count(k,j,i)-2
385	IF ( particles(is+1)%radius == &
386	particles(is)%radius ) THEN
387	sl_r3 = sl_r3 + particles(is+1)%radius**3
388	sl_r4 = sl_r4 + particles(is+1)%radius**4
389	ELSE
390	EXIT
391	ENDIF
392	ENDDO
393
394	IF ( ( s_r3 + sl_r3 ) > 0.0 ) THEN
395
396	mean_r = ( s_r4 + sl_r4 ) / ( s_r3 + sl_r3 )
397
398	CALL collision_efficiency( mean_r, &
399	particles(n)%radius, &
400	effective_coll_efficiency )
401
402	ELSE
403	effective_coll_efficiency = 0.0
404	ENDIF
405
406	!
407	!-- Contribution of the current particle to the next one
408	s_r3 = s_r3 + particles(n)%radius**3
409	s_r4 = s_r4 + particles(n)%radius**4
410
411	IF ( effective_coll_efficiency > 1.0 .OR. &
412	effective_coll_efficiency < 0.0 ) &
413	THEN
414	print*,'+++ advec_particles collision_efficiency ', &
415	'out of range:', effective_coll_efficiency
416	#if defined( __parallel )
417	CALL MPI_ABORT( comm2d, 9999, ierr )
418	#else
419	STOP
420	#endif
421	ENDIF
422
423	!
424	!-- Interpolation of ...
425	ii = particles(n)%x * ddx
426	jj = particles(n)%y * ddy
427	kk = ( particles(n)%z + 0.5 * dz ) / dz
428
429	x = particles(n)%x - ii * dx
430	y = particles(n)%y - jj * dy
431	aa = x2 + y2
432	bb = ( dx - x )2 + y2
433	cc = x2 + ( dy - y )2
434	dd = ( dx - x )2 + ( dy - y )2
435	gg = aa + bb + cc + dd
436
437	ql_int_l = ( ( gg-aa ) * ql(kk,jj,ii) + ( gg-bb ) * &
438	ql(kk,jj,ii+1) &
439	+ ( gg-cc ) * ql(kk,jj+1,ii) + ( gg-dd ) * &
440	ql(kk,jj+1,ii+1) &
441	) / ( 3.0 * gg )
442
443	ql_int_u = ( ( gg-aa ) * ql(kk+1,jj,ii) + ( gg-bb ) * &
444	ql(kk+1,jj,ii+1) &
445	+ ( gg-cc ) * ql(kk+1,jj+1,ii) + ( gg-dd ) * &
446	ql(kk+1,jj+1,ii+1) &
447	) / ( 3.0 * gg )
448
449	ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz * &
450	( ql_int_u - ql_int_l )
451
452	!
453	!-- Interpolate u velocity-component
454	ii = ( particles(n)%x + 0.5 * dx ) * ddx
455	jj = particles(n)%y * ddy
456	kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eq.dist
457
458	IF ( ( particles(n)%z - zu(kk) ) > ( 0.5*dz ) ) kk = kk+1
459
460	x = particles(n)%x + ( 0.5 - ii ) * dx
461	y = particles(n)%y - jj * dy
462	aa = x2 + y2
463	bb = ( dx - x )2 + y2
464	cc = x2 + ( dy - y )2
465	dd = ( dx - x )2 + ( dy - y )2
466	gg = aa + bb + cc + dd
467
468	u_int_l = ( ( gg-aa ) * u(kk,jj,ii) + ( gg-bb ) * &
469	u(kk,jj,ii+1) &
470	+ ( gg-cc ) * u(kk,jj+1,ii) + ( gg-dd ) * &
471	u(kk,jj+1,ii+1) &
472	) / ( 3.0 * gg ) - u_gtrans
473	IF ( kk+1 == nzt+1 ) THEN
474	u_int = u_int_l
475	ELSE
476	u_int_u = ( ( gg-aa ) * u(kk+1,jj,ii) + ( gg-bb ) * &
477	u(kk+1,jj,ii+1) &
478	+ ( gg-cc ) * u(kk+1,jj+1,ii) + ( gg-dd ) * &
479	u(kk+1,jj+1,ii+1) &
480	) / ( 3.0 * gg ) - u_gtrans
481	u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz * &
482	( u_int_u - u_int_l )
483	ENDIF
484
485	!
486	!-- Same procedure for interpolation of the v velocity-compo-
487	!-- nent (adopt index k from u velocity-component)
488	ii = particles(n)%x * ddx
489	jj = ( particles(n)%y + 0.5 * dy ) * ddy
490
491	x = particles(n)%x - ii * dx
492	y = particles(n)%y + ( 0.5 - jj ) * dy
493	aa = x2 + y2
494	bb = ( dx - x )2 + y2
495	cc = x2 + ( dy - y )2
496	dd = ( dx - x )2 + ( dy - y )2
497	gg = aa + bb + cc + dd
498
499	v_int_l = ( ( gg-aa ) * v(kk,jj,ii) + ( gg-bb ) * &
500	v(kk,jj,ii+1) &
501	+ ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) * &
502	v(kk,jj+1,ii+1) &
503	) / ( 3.0 * gg ) - v_gtrans
504	IF ( kk+1 == nzt+1 ) THEN
505	v_int = v_int_l
506	ELSE
507	v_int_u = ( ( gg-aa ) * v(kk+1,jj,ii) + ( gg-bb ) * &
508	v(kk+1,jj,ii+1) &
509	+ ( gg-cc ) * v(kk+1,jj+1,ii) + ( gg-dd ) * &
510	v(kk+1,jj+1,ii+1) &
511	) / ( 3.0 * gg ) - v_gtrans
512	v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz * &
513	( v_int_u - v_int_l )
514	ENDIF
515
516	!
517	!-- Same procedure for interpolation of the w velocity-compo-
518	!-- nent (adopt index i from v velocity-component)
519	jj = particles(n)%y * ddy
520	kk = particles(n)%z / dz
521
522	x = particles(n)%x - ii * dx
523	y = particles(n)%y - jj * dy
524	aa = x2 + y2
525	bb = ( dx - x )2 + y2
526	cc = x2 + ( dy - y )2
527	dd = ( dx - x )2 + ( dy - y )2
528	gg = aa + bb + cc + dd
529
530	w_int_l = ( ( gg-aa ) * w(kk,jj,ii) + ( gg-bb ) * &
531	w(kk,jj,ii+1) &
532	+ ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) * &
533	w(kk,jj+1,ii+1) &
534	) / ( 3.0 * gg )
535	IF ( kk+1 == nzt+1 ) THEN
536	w_int = w_int_l
537	ELSE
538	w_int_u = ( ( gg-aa ) * w(kk+1,jj,ii) + ( gg-bb ) * &
539	w(kk+1,jj,ii+1) &
540	+ ( gg-cc ) * w(kk+1,jj+1,ii) + ( gg-dd ) * &
541	w(kk+1,jj+1,ii+1) &
542	) / ( 3.0 * gg )
543	w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz * &
544	( w_int_u - w_int_l )
545	ENDIF
546
547	!
548	!-- Change in radius due to collision
549	delta_r = effective_coll_efficiency * &
550	ql_int * rho_surface / ( 1.0 - ql_int ) * &
551	0.25 / rho_l * &
552	SQRT( ( u_int - particles(n)%speed_x )**2 + &
553	( v_int - particles(n)%speed_y )**2 + &
554	( w_int - particles(n)%speed_z )**2 &
555	) * dt_3d
556
557	particles(n)%radius = particles(n)%radius + delta_r
558
559	ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%radius**3
560
561	ENDDO
562
563	ENDIF
564
565	!
566	!-- Re-calculate the weighting factor (total liquid water content
567	!-- must be conserved during collision)
568	IF ( ql_vp(k,j,i) /= 0.0 ) THEN
569
570	ql_vp(k,j,i) = ql_v(k,j,i) / ql_vp(k,j,i)
571	!
572	!-- Re-assign this weighting factor to the particles of the
573	!-- current gridbox
574	psi = prt_start_index(k,j,i)
575	DO n = psi, psi + prt_count(k,j,i)-1
576	particles(n)%weight_factor = ql_vp(k,j,i)
577	ENDDO
578
579	ENDIF
580
581	ENDDO
582	ENDDO
583	ENDDO
584
585	CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' )
586
587	ENDIF
588
589
590	!
591	!-- Particle advection.
592	!-- In case of including the SGS velocities, the LES timestep has probably
593	!-- to be split into several smaller timesteps because of the Lagrangian
594	!-- timescale condition. Because the number of timesteps to be carried out is
595	!-- not known at the beginning, these steps are carried out in an infinite loop
596	!-- with exit condition.
597	!
598	!-- If SGS velocities are used, gradients of the TKE have to be calculated and
599	!-- boundary conditions have to be set first. Also, horizontally averaged
600	!-- profiles of the SGS TKE and the resolved-scale velocity variances are
601	!-- needed.
602	IF ( use_sgs_for_particles ) THEN
603
604	!
605	!-- TKE gradient along x and y
606	DO i = nxl, nxr
607	DO j = nys, nyn
608	DO k = nzb, nzt+1
609
610	IF ( k <= nzb_s_inner(j,i-1) .AND. &
611	k > nzb_s_inner(j,i) .AND. &
612	k > nzb_s_inner(j,i+1) ) THEN
613	de_dx(k,j,i) = 2.0 * sgs_wfu_part * &
614	( e(k,j,i+1) - e(k,j,i) ) * ddx
615	ELSEIF ( k > nzb_s_inner(j,i-1) .AND. &
616	k > nzb_s_inner(j,i) .AND. &
617	k <= nzb_s_inner(j,i+1) ) THEN
618	de_dx(k,j,i) = 2.0 * sgs_wfu_part * &
619	( e(k,j,i) - e(k,j,i-1) ) * ddx
620	ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j,i+1) ) &
621	THEN
622	de_dx(k,j,i) = 0.0
623	ELSEIF ( k < nzb_s_inner(j,i-1) .AND. k < nzb_s_inner(j,i) ) &
624	THEN
625	de_dx(k,j,i) = 0.0
626	ELSE
627	de_dx(k,j,i) = sgs_wfu_part * &
628	( e(k,j,i+1) - e(k,j,i-1) ) * ddx
629	ENDIF
630
631	IF ( k <= nzb_s_inner(j-1,i) .AND. &
632	k > nzb_s_inner(j,i) .AND. &
633	k > nzb_s_inner(j+1,i) ) THEN
634	de_dy(k,j,i) = 2.0 * sgs_wfv_part * &
635	( e(k,j+1,i) - e(k,j,i) ) * ddy
636	ELSEIF ( k > nzb_s_inner(j-1,i) .AND. &
637	k > nzb_s_inner(j,i) .AND. &
638	k <= nzb_s_inner(j+1,i) ) THEN
639	de_dy(k,j,i) = 2.0 * sgs_wfv_part * &
640	( e(k,j,i) - e(k,j-1,i) ) * ddy
641	ELSEIF ( k < nzb_s_inner(j,i) .AND. k < nzb_s_inner(j+1,i) ) &
642	THEN
643	de_dy(k,j,i) = 0.0
644	ELSEIF ( k < nzb_s_inner(j-1,i) .AND. k < nzb_s_inner(j,i) ) &
645	THEN
646	de_dy(k,j,i) = 0.0
647	ELSE
648	de_dy(k,j,i) = sgs_wfv_part * &
649	( e(k,j+1,i) - e(k,j-1,i) ) * ddy
650	ENDIF
651
652	ENDDO
653	ENDDO
654	ENDDO
655
656	!
657	!-- TKE gradient along z, including bottom and top boundary conditions
658	DO i = nxl, nxr
659	DO j = nys, nyn
660
661	DO k = nzb_s_inner(j,i)+2, nzt-1
662	de_dz(k,j,i) = 2.0 * sgs_wfw_part * &
663	( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) )
664	ENDDO
665
666	k = nzb_s_inner(j,i)
667	de_dz(nzb:k,j,i) = 0.0
668	de_dz(k+1,j,i) = 2.0 * sgs_wfw_part * ( e(k+2,j,i) - e(k+1,j,i) ) &
669	/ ( zu(k+2) - zu(k+1) )
670	de_dz(nzt,j,i) = 0.0
671	de_dz(nzt+1,j,i) = 0.0
672	ENDDO
673	ENDDO
674
675	!
676	!-- Lateral boundary conditions
677	CALL exchange_horiz( de_dx )
678	CALL exchange_horiz( de_dy )
679	CALL exchange_horiz( de_dz )
680	CALL exchange_horiz( diss )
681
682	!
683	!-- Calculate the horizontally averaged profiles of SGS TKE and resolved
684	!-- velocity variances (they may have been already calculated in routine
685	!-- flow_statistics).
686	IF ( .NOT. flow_statistics_called ) THEN
687	!
688	!-- First calculate horizontally averaged profiles of the horizontal
689	!-- velocities.
690	sums_l(:,1,0) = 0.0
691	sums_l(:,2,0) = 0.0
692
693	DO i = nxl, nxr
694	DO j = nys, nyn
695	DO k = nzb_s_outer(j,i), nzt+1
696	sums_l(k,1,0) = sums_l(k,1,0) + u(k,j,i)
697	sums_l(k,2,0) = sums_l(k,2,0) + v(k,j,i)
698	ENDDO
699	ENDDO
700	ENDDO
701
702	#if defined( __parallel )
703	!
704	!-- Compute total sum from local sums
705	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, &
706	MPI_REAL, MPI_SUM, comm2d, ierr )
707	CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, &
708	MPI_REAL, MPI_SUM, comm2d, ierr )
709	#else
710	sums(:,1) = sums_l(:,1,0)
711	sums(:,2) = sums_l(:,2,0)
712	#endif
713
714	!
715	!-- Final values are obtained by division by the total number of grid
716	!-- points used for the summation.
717	hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0) ! u
718	hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0) ! v
719
720	!
721	!-- Now calculate the profiles of SGS TKE and the resolved-scale
722	!-- velocity variances
723	sums_l(:,8,0) = 0.0
724	sums_l(:,30,0) = 0.0
725	sums_l(:,31,0) = 0.0
726	sums_l(:,32,0) = 0.0
727	DO i = nxl, nxr
728	DO j = nys, nyn
729	DO k = nzb_s_outer(j,i), nzt+1
730	sums_l(k,8,0) = sums_l(k,8,0) + e(k,j,i)
731	sums_l(k,30,0) = sums_l(k,30,0) + &
732	( u(k,j,i) - hom(k,1,1,0) )**2
733	sums_l(k,31,0) = sums_l(k,31,0) + &
734	( v(k,j,i) - hom(k,1,2,0) )**2
735	sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)**2
736	ENDDO
737	ENDDO
738	ENDDO
739
740	#if defined( __parallel )
741	!
742	!-- Compute total sum from local sums
743	CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, &
744	MPI_REAL, MPI_SUM, comm2d, ierr )
745	CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, &
746	MPI_REAL, MPI_SUM, comm2d, ierr )
747	CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, &
748	MPI_REAL, MPI_SUM, comm2d, ierr )
749	CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
750	MPI_REAL, MPI_SUM, comm2d, ierr )
751
752	#else
753	sums(:,8) = sums_l(:,8,0)
754	sums(:,30) = sums_l(:,30,0)
755	sums(:,31) = sums_l(:,31,0)
756	sums(:,32) = sums_l(:,32,0)
757	#endif
758
759	!
760	!-- Final values are obtained by division by the total number of grid
761	!-- points used for the summation.
762	hom(:,1,8,0) = sums(:,8) / ngp_2dh_outer(:,0) ! e
763	hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0) ! u*2
764	hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0) ! v*2
765	hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0) ! w*2
766
767	ENDIF
768
769	ENDIF
770
771	!
772	!-- Initialize variables used for accumulating the number of particles
773	!-- exchanged between the subdomains during all sub-timesteps (if sgs
774	!-- velocities are included). These data are output further below on the
775	!-- particle statistics file.
776	trlp_count_sum = 0
777	trlp_count_recv_sum = 0
778	trrp_count_sum = 0
779	trrp_count_recv_sum = 0
780	trsp_count_sum = 0
781	trsp_count_recv_sum = 0
782	trnp_count_sum = 0
783	trnp_count_recv_sum = 0
784
785	!
786	!-- Initialize the variable storing the total time that a particle has advanced
787	!-- within the timestep procedure
788	particles(1:number_of_particles)%dt_sum = 0.0
789
790	!
791	!-- Timestep loop.
792	!-- This loop has to be repeated until the advection time of every particle
793	!-- (in the total domain!) has reached the LES timestep (dt_3d)
794	DO
795
796	CALL cpu_log( log_point_s(44), 'advec_part_advec', 'start' )
797
798	!
799	!-- Initialize the switch used for the loop exit condition checked at the
800	!-- end of this loop.
801	!-- If at least one particle has failed to reach the LES timestep, this
802	!-- switch will be set false.
803	dt_3d_reached_l = .TRUE.
804
805	!
806	!-- Initialize variables for the (sub-) timestep, i.e. for marking those
807	!-- particles to be deleted after the timestep
808	particle_mask = .TRUE.
809	deleted_particles = 0
810	trlp_count_recv = 0
811	trnp_count_recv = 0
812	trrp_count_recv = 0
813	trsp_count_recv = 0
814	IF ( use_particle_tails ) THEN
815	tail_mask = .TRUE.
816	deleted_tails = 0
817	ENDIF
818
819
820	DO n = 1, number_of_particles
821	!
822	!-- Move particles only if the LES timestep has not (approximately) been
823	!-- reached
824	IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 ) CYCLE
825
826	!
827	!-- Interpolate u velocity-component, determine left, front, bottom
828	!-- index of u-array
829	i = ( particles(n)%x + 0.5 * dx ) * ddx
830	j = particles(n)%y * ddy
831	k = ( particles(n)%z + 0.5 * dz ) / dz ! only exact if equidistant
832
833	!
834	!-- Interpolation of the velocity components in the xy-plane
835	x = particles(n)%x + ( 0.5 - i ) * dx
836	y = particles(n)%y - j * dy
837	aa = x2 + y2
838	bb = ( dx - x )2 + y2
839	cc = x2 + ( dy - y )2
840	dd = ( dx - x )2 + ( dy - y )2
841	gg = aa + bb + cc + dd
842
843	u_int_l = ( ( gg - aa ) * u(k,j,i) + ( gg - bb ) * u(k,j,i+1) &
844	+ ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
845	) / ( 3.0 * gg ) - u_gtrans
846	IF ( k+1 == nzt+1 ) THEN
847	u_int = u_int_l
848	ELSE
849	u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1) &
850	+ ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
851	) / ( 3.0 * gg ) - u_gtrans
852	u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
853	( u_int_u - u_int_l )
854	ENDIF
855
856	!
857	!-- Same procedure for interpolation of the v velocity-component (adopt
858	!-- index k from u velocity-component)
859	i = particles(n)%x * ddx
860	j = ( particles(n)%y + 0.5 * dy ) * ddy
861
862	x = particles(n)%x - i * dx
863	y = particles(n)%y + ( 0.5 - j ) * dy
864	aa = x2 + y2
865	bb = ( dx - x )2 + y2
866	cc = x2 + ( dy - y )2
867	dd = ( dx - x )2 + ( dy - y )2
868	gg = aa + bb + cc + dd
869
870	v_int_l = ( ( gg - aa ) * v(k,j,i) + ( gg - bb ) * v(k,j,i+1) &
871	+ ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
872	) / ( 3.0 * gg ) - v_gtrans
873	IF ( k+1 == nzt+1 ) THEN
874	v_int = v_int_l
875	ELSE
876	v_int_u = ( ( gg-aa ) * v(k+1,j,i) + ( gg-bb ) * v(k+1,j,i+1) &
877	+ ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
878	) / ( 3.0 * gg ) - v_gtrans
879	v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
880	( v_int_u - v_int_l )
881	ENDIF
882
883	!
884	!-- Same procedure for interpolation of the w velocity-component (adopt
885	!-- index i from v velocity-component)
886	IF ( vertical_particle_advection ) THEN
887	j = particles(n)%y * ddy
888	k = particles(n)%z / dz
889
890	x = particles(n)%x - i * dx
891	y = particles(n)%y - j * dy
892	aa = x2 + y2
893	bb = ( dx - x )2 + y2
894	cc = x2 + ( dy - y )2
895	dd = ( dx - x )2 + ( dy - y )2
896	gg = aa + bb + cc + dd
897
898	w_int_l = ( ( gg - aa ) * w(k,j,i) + ( gg - bb ) * w(k,j,i+1) &
899	+ ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
900	) / ( 3.0 * gg )
901	IF ( k+1 == nzt+1 ) THEN
902	w_int = w_int_l
903	ELSE
904	w_int_u = ( ( gg-aa ) * w(k+1,j,i) + &
905	( gg-bb ) * w(k+1,j,i+1) + &
906	( gg-cc ) * w(k+1,j+1,i) + &
907	( gg-dd ) * w(k+1,j+1,i+1) &
908	) / ( 3.0 * gg )
909	w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
910	( w_int_u - w_int_l )
911	ENDIF
912	ELSE
913	w_int = 0.0
914	ENDIF
915
916	!
917	!-- Interpolate and calculate quantities needed for calculating the SGS
918	!-- velocities
919	IF ( use_sgs_for_particles ) THEN
920	!
921	!-- Interpolate TKE
922	i = particles(n)%x * ddx
923	j = particles(n)%y * ddy
924	k = ( particles(n)%z + 0.5 * dz ) / dz ! only exact if eq.dist
925
926	IF ( topography == 'flat' ) THEN
927
928	x = particles(n)%x - i * dx
929	y = particles(n)%y - j * dy
930	aa = x2 + y2
931	bb = ( dx - x )2 + y2
932	cc = x2 + ( dy - y )2
933	dd = ( dx - x )2 + ( dy - y )2
934	gg = aa + bb + cc + dd
935
936	e_int_l = ( ( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
937	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
938	) / ( 3.0 * gg )
939
940	IF ( k+1 == nzt+1 ) THEN
941	e_int = e_int_l
942	ELSE
943	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
944	( gg - bb ) * e(k+1,j,i+1) + &
945	( gg - cc ) * e(k+1,j+1,i) + &
946	( gg - dd ) * e(k+1,j+1,i+1) &
947	) / ( 3.0 * gg )
948	e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
949	( e_int_u - e_int_l )
950	ENDIF
951
952	!
953	!-- Interpolate the TKE gradient along x (adopt incides i,j,k and
954	!-- all position variables from above (TKE))
955	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
956	( gg - bb ) * de_dx(k,j,i+1) + &
957	( gg - cc ) * de_dx(k,j+1,i) + &
958	( gg - dd ) * de_dx(k,j+1,i+1) &
959	) / ( 3.0 * gg )
960
961	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
962	de_dx_int = de_dx_int_l
963	ELSE
964	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
965	( gg - bb ) * de_dx(k+1,j,i+1) + &
966	( gg - cc ) * de_dx(k+1,j+1,i) + &
967	( gg - dd ) * de_dx(k+1,j+1,i+1) &
968	) / ( 3.0 * gg )
969	de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
970	( de_dx_int_u - de_dx_int_l )
971	ENDIF
972
973	!
974	!-- Interpolate the TKE gradient along y
975	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
976	( gg - bb ) * de_dy(k,j,i+1) + &
977	( gg - cc ) * de_dy(k,j+1,i) + &
978	( gg - dd ) * de_dy(k,j+1,i+1) &
979	) / ( 3.0 * gg )
980	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
981	de_dy_int = de_dy_int_l
982	ELSE
983	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
984	( gg - bb ) * de_dy(k+1,j,i+1) + &
985	( gg - cc ) * de_dy(k+1,j+1,i) + &
986	( gg - dd ) * de_dy(k+1,j+1,i+1) &
987	) / ( 3.0 * gg )
988	de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
989	( de_dy_int_u - de_dy_int_l )
990	ENDIF
991
992	!
993	!-- Interpolate the TKE gradient along z
994	IF ( particles(n)%z < 0.5 * dz ) THEN
995	de_dz_int = 0.0
996	ELSE
997	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
998	( gg - bb ) * de_dz(k,j,i+1) + &
999	( gg - cc ) * de_dz(k,j+1,i) + &
1000	( gg - dd ) * de_dz(k,j+1,i+1) &
1001	) / ( 3.0 * gg )
1002
1003	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
1004	de_dz_int = de_dz_int_l
1005	ELSE
1006	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
1007	( gg - bb ) * de_dz(k+1,j,i+1) + &
1008	( gg - cc ) * de_dz(k+1,j+1,i) + &
1009	( gg - dd ) * de_dz(k+1,j+1,i+1) &
1010	) / ( 3.0 * gg )
1011	de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
1012	( de_dz_int_u - de_dz_int_l )
1013	ENDIF
1014	ENDIF
1015
1016	!
1017	!-- Interpolate the dissipation of TKE
1018	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
1019	( gg - bb ) * diss(k,j,i+1) + &
1020	( gg - cc ) * diss(k,j+1,i) + &
1021	( gg - dd ) * diss(k,j+1,i+1) &
1022	) / ( 3.0 * gg )
1023
1024	IF ( k+1 == nzt+1 ) THEN
1025	diss_int = diss_int_l
1026	ELSE
1027	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
1028	( gg - bb ) * diss(k+1,j,i+1) + &
1029	( gg - cc ) * diss(k+1,j+1,i) + &
1030	( gg - dd ) * diss(k+1,j+1,i+1) &
1031	) / ( 3.0 * gg )
1032	diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
1033	( diss_int_u - diss_int_l )
1034	ENDIF
1035
1036	ELSE
1037
1038	!
1039	!-- In case that there are buildings it has to be determined
1040	!-- how many of the gridpoints defining the particle box are
1041	!-- situated within a building
1042	!-- gp_outside_of_building(1): i,j,k
1043	!-- gp_outside_of_building(2): i,j+1,k
1044	!-- gp_outside_of_building(3): i,j,k+1
1045	!-- gp_outside_of_building(4): i,j+1,k+1
1046	!-- gp_outside_of_building(5): i+1,j,k
1047	!-- gp_outside_of_building(6): i+1,j+1,k
1048	!-- gp_outside_of_building(7): i+1,j,k+1
1049	!-- gp_outside_of_building(8): i+1,j+1,k+1
1050
1051	gp_outside_of_building = 0
1052	location = 0.0
1053	num_gp = 0
1054
1055	IF ( k > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
1056	num_gp = num_gp + 1
1057	gp_outside_of_building(1) = 1
1058	location(num_gp,1) = i * dx
1059	location(num_gp,2) = j * dy
1060	location(num_gp,3) = k * dz - 0.5 * dz
1061	ei(num_gp) = e(k,j,i)
1062	dissi(num_gp) = diss(k,j,i)
1063	de_dxi(num_gp) = de_dx(k,j,i)
1064	de_dyi(num_gp) = de_dy(k,j,i)
1065	de_dzi(num_gp) = de_dz(k,j,i)
1066	ENDIF
1067
1068	IF ( k > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) &
1069	THEN
1070	num_gp = num_gp + 1
1071	gp_outside_of_building(2) = 1
1072	location(num_gp,1) = i * dx
1073	location(num_gp,2) = (j+1) * dy
1074	location(num_gp,3) = k * dz - 0.5 * dz
1075	ei(num_gp) = e(k,j+1,i)
1076	dissi(num_gp) = diss(k,j+1,i)
1077	de_dxi(num_gp) = de_dx(k,j+1,i)
1078	de_dyi(num_gp) = de_dy(k,j+1,i)
1079	de_dzi(num_gp) = de_dz(k,j+1,i)
1080	ENDIF
1081
1082	IF ( k+1 > nzb_s_inner(j,i) .OR. nzb_s_inner(j,i) == 0 ) THEN
1083	num_gp = num_gp + 1
1084	gp_outside_of_building(3) = 1
1085	location(num_gp,1) = i * dx
1086	location(num_gp,2) = j * dy
1087	location(num_gp,3) = (k+1) * dz - 0.5 * dz
1088	ei(num_gp) = e(k+1,j,i)
1089	dissi(num_gp) = diss(k+1,j,i)
1090	de_dxi(num_gp) = de_dx(k+1,j,i)
1091	de_dyi(num_gp) = de_dy(k+1,j,i)
1092	de_dzi(num_gp) = de_dz(k+1,j,i)
1093	ENDIF
1094
1095	IF ( k+1 > nzb_s_inner(j+1,i) .OR. nzb_s_inner(j+1,i) == 0 ) &
1096	THEN
1097	num_gp = num_gp + 1
1098	gp_outside_of_building(4) = 1
1099	location(num_gp,1) = i * dx
1100	location(num_gp,2) = (j+1) * dy
1101	location(num_gp,3) = (k+1) * dz - 0.5 * dz
1102	ei(num_gp) = e(k+1,j+1,i)
1103	dissi(num_gp) = diss(k+1,j+1,i)
1104	de_dxi(num_gp) = de_dx(k+1,j+1,i)
1105	de_dyi(num_gp) = de_dy(k+1,j+1,i)
1106	de_dzi(num_gp) = de_dz(k+1,j+1,i)
1107	ENDIF
1108
1109	IF ( k > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) &
1110	THEN
1111	num_gp = num_gp + 1
1112	gp_outside_of_building(5) = 1
1113	location(num_gp,1) = (i+1) * dx
1114	location(num_gp,2) = j * dy
1115	location(num_gp,3) = k * dz - 0.5 * dz
1116	ei(num_gp) = e(k,j,i+1)
1117	dissi(num_gp) = diss(k,j,i+1)
1118	de_dxi(num_gp) = de_dx(k,j,i+1)
1119	de_dyi(num_gp) = de_dy(k,j,i+1)
1120	de_dzi(num_gp) = de_dz(k,j,i+1)
1121	ENDIF
1122
1123	IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) &
1124	THEN
1125	num_gp = num_gp + 1
1126	gp_outside_of_building(6) = 1
1127	location(num_gp,1) = (i+1) * dx
1128	location(num_gp,2) = (j+1) * dy
1129	location(num_gp,3) = k * dz - 0.5 * dz
1130	ei(num_gp) = e(k,j+1,i+1)
1131	dissi(num_gp) = diss(k,j+1,i+1)
1132	de_dxi(num_gp) = de_dx(k,j+1,i+1)
1133	de_dyi(num_gp) = de_dy(k,j+1,i+1)
1134	de_dzi(num_gp) = de_dz(k,j+1,i+1)
1135	ENDIF
1136
1137	IF ( k+1 > nzb_s_inner(j,i+1) .OR. nzb_s_inner(j,i+1) == 0 ) &
1138	THEN
1139	num_gp = num_gp + 1
1140	gp_outside_of_building(7) = 1
1141	location(num_gp,1) = (i+1) * dx
1142	location(num_gp,2) = j * dy
1143	location(num_gp,3) = (k+1) * dz - 0.5 * dz
1144	ei(num_gp) = e(k+1,j,i+1)
1145	dissi(num_gp) = diss(k+1,j,i+1)
1146	de_dxi(num_gp) = de_dx(k+1,j,i+1)
1147	de_dyi(num_gp) = de_dy(k+1,j,i+1)
1148	de_dzi(num_gp) = de_dz(k+1,j,i+1)
1149	ENDIF
1150
1151	IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)&
1152	THEN
1153	num_gp = num_gp + 1
1154	gp_outside_of_building(8) = 1
1155	location(num_gp,1) = (i+1) * dx
1156	location(num_gp,2) = (j+1) * dy
1157	location(num_gp,3) = (k+1) * dz - 0.5 * dz
1158	ei(num_gp) = e(k+1,j+1,i+1)
1159	dissi(num_gp) = diss(k+1,j+1,i+1)
1160	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
1161	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
1162	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
1163	ENDIF
1164
1165	!
1166	!-- If all gridpoints are situated outside of a building, then the
1167	!-- ordinary interpolation scheme can be used.
1168	IF ( num_gp == 8 ) THEN
1169
1170	x = particles(n)%x - i * dx
1171	y = particles(n)%y - j * dy
1172	aa = x2 + y2
1173	bb = ( dx - x )2 + y2
1174	cc = x2 + ( dy - y )2
1175	dd = ( dx - x )2 + ( dy - y )2
1176	gg = aa + bb + cc + dd
1177
1178	e_int_l = (( gg-aa ) * e(k,j,i) + ( gg-bb ) * e(k,j,i+1) &
1179	+ ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)&
1180	) / ( 3.0 * gg )
1181
1182	IF ( k+1 == nzt+1 ) THEN
1183	e_int = e_int_l
1184	ELSE
1185	e_int_u = ( ( gg - aa ) * e(k+1,j,i) + &
1186	( gg - bb ) * e(k+1,j,i+1) + &
1187	( gg - cc ) * e(k+1,j+1,i) + &
1188	( gg - dd ) * e(k+1,j+1,i+1) &
1189	) / ( 3.0 * gg )
1190	e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
1191	( e_int_u - e_int_l )
1192	ENDIF
1193
1194	!
1195	!-- Interpolate the TKE gradient along x (adopt incides i,j,k
1196	!-- and all position variables from above (TKE))
1197	de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i) + &
1198	( gg - bb ) * de_dx(k,j,i+1) + &
1199	( gg - cc ) * de_dx(k,j+1,i) + &
1200	( gg - dd ) * de_dx(k,j+1,i+1) &
1201	) / ( 3.0 * gg )
1202
1203	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
1204	de_dx_int = de_dx_int_l
1205	ELSE
1206	de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i) + &
1207	( gg - bb ) * de_dx(k+1,j,i+1) + &
1208	( gg - cc ) * de_dx(k+1,j+1,i) + &
1209	( gg - dd ) * de_dx(k+1,j+1,i+1) &
1210	) / ( 3.0 * gg )
1211	de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / &
1212	dz * ( de_dx_int_u - de_dx_int_l )
1213	ENDIF
1214
1215	!
1216	!-- Interpolate the TKE gradient along y
1217	de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i) + &
1218	( gg - bb ) * de_dy(k,j,i+1) + &
1219	( gg - cc ) * de_dy(k,j+1,i) + &
1220	( gg - dd ) * de_dy(k,j+1,i+1) &
1221	) / ( 3.0 * gg )
1222	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
1223	de_dy_int = de_dy_int_l
1224	ELSE
1225	de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i) + &
1226	( gg - bb ) * de_dy(k+1,j,i+1) + &
1227	( gg - cc ) * de_dy(k+1,j+1,i) + &
1228	( gg - dd ) * de_dy(k+1,j+1,i+1) &
1229	) / ( 3.0 * gg )
1230	de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / &
1231	dz * ( de_dy_int_u - de_dy_int_l )
1232	ENDIF
1233
1234	!
1235	!-- Interpolate the TKE gradient along z
1236	IF ( particles(n)%z < 0.5 * dz ) THEN
1237	de_dz_int = 0.0
1238	ELSE
1239	de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i) + &
1240	( gg - bb ) * de_dz(k,j,i+1) + &
1241	( gg - cc ) * de_dz(k,j+1,i) + &
1242	( gg - dd ) * de_dz(k,j+1,i+1) &
1243	) / ( 3.0 * gg )
1244
1245	IF ( ( k+1 == nzt+1 ) .OR. ( k == nzb ) ) THEN
1246	de_dz_int = de_dz_int_l
1247	ELSE
1248	de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i) + &
1249	( gg - bb ) * de_dz(k+1,j,i+1) + &
1250	( gg - cc ) * de_dz(k+1,j+1,i) + &
1251	( gg - dd ) * de_dz(k+1,j+1,i+1) &
1252	) / ( 3.0 * gg )
1253	de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /&
1254	dz * ( de_dz_int_u - de_dz_int_l )
1255	ENDIF
1256	ENDIF
1257
1258	!
1259	!-- Interpolate the dissipation of TKE
1260	diss_int_l = ( ( gg - aa ) * diss(k,j,i) + &
1261	( gg - bb ) * diss(k,j,i+1) + &
1262	( gg - cc ) * diss(k,j+1,i) + &
1263	( gg - dd ) * diss(k,j+1,i+1) &
1264	) / ( 3.0 * gg )
1265
1266	IF ( k+1 == nzt+1 ) THEN
1267	diss_int = diss_int_l
1268	ELSE
1269	diss_int_u = ( ( gg - aa ) * diss(k+1,j,i) + &
1270	( gg - bb ) * diss(k+1,j,i+1) + &
1271	( gg - cc ) * diss(k+1,j+1,i) + &
1272	( gg - dd ) * diss(k+1,j+1,i+1) &
1273	) / ( 3.0 * gg )
1274	diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *&
1275	( diss_int_u - diss_int_l )
1276	ENDIF
1277
1278	ELSE
1279
1280	!
1281	!-- If wall between gridpoint 1 and gridpoint 5, then
1282	!-- Neumann boundary condition has to be applied
1283	IF ( gp_outside_of_building(1) == 1 .AND. &
1284	gp_outside_of_building(5) == 0 ) THEN
1285	num_gp = num_gp + 1
1286	location(num_gp,1) = i * dx + 0.5 * dx
1287	location(num_gp,2) = j * dy
1288	location(num_gp,3) = k * dz - 0.5 * dz
1289	ei(num_gp) = e(k,j,i)
1290	dissi(num_gp) = diss(k,j,i)
1291	de_dxi(num_gp) = 0.0
1292	de_dyi(num_gp) = de_dy(k,j,i)
1293	de_dzi(num_gp) = de_dz(k,j,i)
1294	ENDIF
1295
1296	IF ( gp_outside_of_building(5) == 1 .AND. &
1297	gp_outside_of_building(1) == 0 ) THEN
1298	num_gp = num_gp + 1
1299	location(num_gp,1) = i * dx + 0.5 * dx
1300	location(num_gp,2) = j * dy
1301	location(num_gp,3) = k * dz - 0.5 * dz
1302	ei(num_gp) = e(k,j,i+1)
1303	dissi(num_gp) = diss(k,j,i+1)
1304	de_dxi(num_gp) = 0.0
1305	de_dyi(num_gp) = de_dy(k,j,i+1)
1306	de_dzi(num_gp) = de_dz(k,j,i+1)
1307	ENDIF
1308
1309	!
1310	!-- If wall between gridpoint 5 and gridpoint 6, then
1311	!-- then Neumann boundary condition has to be applied
1312	IF ( gp_outside_of_building(5) == 1 .AND. &
1313	gp_outside_of_building(6) == 0 ) THEN
1314	num_gp = num_gp + 1
1315	location(num_gp,1) = (i+1) * dx
1316	location(num_gp,2) = j * dy + 0.5 * dy
1317	location(num_gp,3) = k * dz - 0.5 * dz
1318	ei(num_gp) = e(k,j,i+1)
1319	dissi(num_gp) = diss(k,j,i+1)
1320	de_dxi(num_gp) = de_dx(k,j,i+1)
1321	de_dyi(num_gp) = 0.0
1322	de_dzi(num_gp) = de_dz(k,j,i+1)
1323	ENDIF
1324
1325	IF ( gp_outside_of_building(6) == 1 .AND. &
1326	gp_outside_of_building(5) == 0 ) THEN
1327	num_gp = num_gp + 1
1328	location(num_gp,1) = (i+1) * dx
1329	location(num_gp,2) = j * dy + 0.5 * dy
1330	location(num_gp,3) = k * dz - 0.5 * dz
1331	ei(num_gp) = e(k,j+1,i+1)
1332	dissi(num_gp) = diss(k,j+1,i+1)
1333	de_dxi(num_gp) = de_dx(k,j+1,i+1)
1334	de_dyi(num_gp) = 0.0
1335	de_dzi(num_gp) = de_dz(k,j+1,i+1)
1336	ENDIF
1337
1338	!
1339	!-- If wall between gridpoint 2 and gridpoint 6, then
1340	!-- Neumann boundary condition has to be applied
1341	IF ( gp_outside_of_building(2) == 1 .AND. &
1342	gp_outside_of_building(6) == 0 ) THEN
1343	num_gp = num_gp + 1
1344	location(num_gp,1) = i * dx + 0.5 * dx
1345	location(num_gp,2) = (j+1) * dy
1346	location(num_gp,3) = k * dz - 0.5 * dz
1347	ei(num_gp) = e(k,j+1,i)
1348	dissi(num_gp) = diss(k,j+1,i)
1349	de_dxi(num_gp) = 0.0
1350	de_dyi(num_gp) = de_dy(k,j+1,i)
1351	de_dzi(num_gp) = de_dz(k,j+1,i)
1352	ENDIF
1353
1354	IF ( gp_outside_of_building(6) == 1 .AND. &
1355	gp_outside_of_building(2) == 0 ) THEN
1356	num_gp = num_gp + 1
1357	location(num_gp,1) = i * dx + 0.5 * dx
1358	location(num_gp,2) = (j+1) * dy
1359	location(num_gp,3) = k * dz - 0.5 * dz
1360	ei(num_gp) = e(k,j+1,i+1)
1361	dissi(num_gp) = diss(k,j+1,i+1)
1362	de_dxi(num_gp) = 0.0
1363	de_dyi(num_gp) = de_dy(k,j+1,i+1)
1364	de_dzi(num_gp) = de_dz(k,j+1,i+1)
1365	ENDIF
1366
1367	!
1368	!-- If wall between gridpoint 1 and gridpoint 2, then
1369	!-- Neumann boundary condition has to be applied
1370	IF ( gp_outside_of_building(1) == 1 .AND. &
1371	gp_outside_of_building(2) == 0 ) THEN
1372	num_gp = num_gp + 1
1373	location(num_gp,1) = i * dx
1374	location(num_gp,2) = j * dy + 0.5 * dy
1375	location(num_gp,3) = k * dz - 0.5 * dz
1376	ei(num_gp) = e(k,j,i)
1377	dissi(num_gp) = diss(k,j,i)
1378	de_dxi(num_gp) = de_dx(k,j,i)
1379	de_dyi(num_gp) = 0.0
1380	de_dzi(num_gp) = de_dz(k,j,i)
1381	ENDIF
1382
1383	IF ( gp_outside_of_building(2) == 1 .AND. &
1384	gp_outside_of_building(1) == 0 ) THEN
1385	num_gp = num_gp + 1
1386	location(num_gp,1) = i * dx
1387	location(num_gp,2) = j * dy + 0.5 * dy
1388	location(num_gp,3) = k * dz - 0.5 * dz
1389	ei(num_gp) = e(k,j+1,i)
1390	dissi(num_gp) = diss(k,j+1,i)
1391	de_dxi(num_gp) = de_dx(k,j+1,i)
1392	de_dyi(num_gp) = 0.0
1393	de_dzi(num_gp) = de_dz(k,j+1,i)
1394	ENDIF
1395
1396	!
1397	!-- If wall between gridpoint 3 and gridpoint 7, then
1398	!-- Neumann boundary condition has to be applied
1399	IF ( gp_outside_of_building(3) == 1 .AND. &
1400	gp_outside_of_building(7) == 0 ) THEN
1401	num_gp = num_gp + 1
1402	location(num_gp,1) = i * dx + 0.5 * dx
1403	location(num_gp,2) = j * dy
1404	location(num_gp,3) = k * dz + 0.5 * dz
1405	ei(num_gp) = e(k+1,j,i)
1406	dissi(num_gp) = diss(k+1,j,i)
1407	de_dxi(num_gp) = 0.0
1408	de_dyi(num_gp) = de_dy(k+1,j,i)
1409	de_dzi(num_gp) = de_dz(k+1,j,i)
1410	ENDIF
1411
1412	IF ( gp_outside_of_building(7) == 1 .AND. &
1413	gp_outside_of_building(3) == 0 ) THEN
1414	num_gp = num_gp + 1
1415	location(num_gp,1) = i * dx + 0.5 * dx
1416	location(num_gp,2) = j * dy
1417	location(num_gp,3) = k * dz + 0.5 * dz
1418	ei(num_gp) = e(k+1,j,i+1)
1419	dissi(num_gp) = diss(k+1,j,i+1)
1420	de_dxi(num_gp) = 0.0
1421	de_dyi(num_gp) = de_dy(k+1,j,i+1)
1422	de_dzi(num_gp) = de_dz(k+1,j,i+1)
1423	ENDIF
1424
1425	!
1426	!-- If wall between gridpoint 7 and gridpoint 8, then
1427	!-- Neumann boundary condition has to be applied
1428	IF ( gp_outside_of_building(7) == 1 .AND. &
1429	gp_outside_of_building(8) == 0 ) THEN
1430	num_gp = num_gp + 1
1431	location(num_gp,1) = (i+1) * dx
1432	location(num_gp,2) = j * dy + 0.5 * dy
1433	location(num_gp,3) = k * dz + 0.5 * dz
1434	ei(num_gp) = e(k+1,j,i+1)
1435	dissi(num_gp) = diss(k+1,j,i+1)
1436	de_dxi(num_gp) = de_dx(k+1,j,i+1)
1437	de_dyi(num_gp) = 0.0
1438	de_dzi(num_gp) = de_dz(k+1,j,i+1)
1439	ENDIF
1440
1441	IF ( gp_outside_of_building(8) == 1 .AND. &
1442	gp_outside_of_building(7) == 0 ) THEN
1443	num_gp = num_gp + 1
1444	location(num_gp,1) = (i+1) * dx
1445	location(num_gp,2) = j * dy + 0.5 * dy
1446	location(num_gp,3) = k * dz + 0.5 * dz
1447	ei(num_gp) = e(k+1,j+1,i+1)
1448	dissi(num_gp) = diss(k+1,j+1,i+1)
1449	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
1450	de_dyi(num_gp) = 0.0
1451	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
1452	ENDIF
1453
1454	!
1455	!-- If wall between gridpoint 4 and gridpoint 8, then
1456	!-- Neumann boundary condition has to be applied
1457	IF ( gp_outside_of_building(4) == 1 .AND. &
1458	gp_outside_of_building(8) == 0 ) THEN
1459	num_gp = num_gp + 1
1460	location(num_gp,1) = i * dx + 0.5 * dx
1461	location(num_gp,2) = (j+1) * dy
1462	location(num_gp,3) = k * dz + 0.5 * dz
1463	ei(num_gp) = e(k+1,j+1,i)
1464	dissi(num_gp) = diss(k+1,j+1,i)
1465	de_dxi(num_gp) = 0.0
1466	de_dyi(num_gp) = de_dy(k+1,j+1,i)
1467	de_dzi(num_gp) = de_dz(k+1,j+1,i)
1468	ENDIF
1469
1470	IF ( gp_outside_of_building(8) == 1 .AND. &
1471	gp_outside_of_building(4) == 0 ) THEN
1472	num_gp = num_gp + 1
1473	location(num_gp,1) = i * dx + 0.5 * dx
1474	location(num_gp,2) = (j+1) * dy
1475	location(num_gp,3) = k * dz + 0.5 * dz
1476	ei(num_gp) = e(k+1,j+1,i+1)
1477	dissi(num_gp) = diss(k+1,j+1,i+1)
1478	de_dxi(num_gp) = 0.0
1479	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
1480	de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
1481	ENDIF
1482
1483	!
1484	!-- If wall between gridpoint 3 and gridpoint 4, then
1485	!-- Neumann boundary condition has to be applied
1486	IF ( gp_outside_of_building(3) == 1 .AND. &
1487	gp_outside_of_building(4) == 0 ) THEN
1488	num_gp = num_gp + 1
1489	location(num_gp,1) = i * dx
1490	location(num_gp,2) = j * dy + 0.5 * dy
1491	location(num_gp,3) = k * dz + 0.5 * dz
1492	ei(num_gp) = e(k+1,j,i)
1493	dissi(num_gp) = diss(k+1,j,i)
1494	de_dxi(num_gp) = de_dx(k+1,j,i)
1495	de_dyi(num_gp) = 0.0
1496	de_dzi(num_gp) = de_dz(k+1,j,i)
1497	ENDIF
1498
1499	IF ( gp_outside_of_building(4) == 1 .AND. &
1500	gp_outside_of_building(3) == 0 ) THEN
1501	num_gp = num_gp + 1
1502	location(num_gp,1) = i * dx
1503	location(num_gp,2) = j * dy + 0.5 * dy
1504	location(num_gp,3) = k * dz + 0.5 * dz
1505	ei(num_gp) = e(k+1,j+1,i)
1506	dissi(num_gp) = diss(k+1,j+1,i)
1507	de_dxi(num_gp) = de_dx(k+1,j+1,i)
1508	de_dyi(num_gp) = 0.0
1509	de_dzi(num_gp) = de_dz(k+1,j+1,i)
1510	ENDIF
1511
1512	!
1513	!-- If wall between gridpoint 1 and gridpoint 3, then
1514	!-- Neumann boundary condition has to be applied
1515	!-- (only one case as only building beneath is possible)
1516	IF ( gp_outside_of_building(1) == 0 .AND. &
1517	gp_outside_of_building(3) == 1 ) THEN
1518	num_gp = num_gp + 1
1519	location(num_gp,1) = i * dx
1520	location(num_gp,2) = j * dy
1521	location(num_gp,3) = k * dz
1522	ei(num_gp) = e(k+1,j,i)
1523	dissi(num_gp) = diss(k+1,j,i)
1524	de_dxi(num_gp) = de_dx(k+1,j,i)
1525	de_dyi(num_gp) = de_dy(k+1,j,i)
1526	de_dzi(num_gp) = 0.0
1527	ENDIF
1528
1529	!
1530	!-- If wall between gridpoint 5 and gridpoint 7, then
1531	!-- Neumann boundary condition has to be applied
1532	!-- (only one case as only building beneath is possible)
1533	IF ( gp_outside_of_building(5) == 0 .AND. &
1534	gp_outside_of_building(7) == 1 ) THEN
1535	num_gp = num_gp + 1
1536	location(num_gp,1) = (i+1) * dx
1537	location(num_gp,2) = j * dy
1538	location(num_gp,3) = k * dz
1539	ei(num_gp) = e(k+1,j,i+1)
1540	dissi(num_gp) = diss(k+1,j,i+1)
1541	de_dxi(num_gp) = de_dx(k+1,j,i+1)
1542	de_dyi(num_gp) = de_dy(k+1,j,i+1)
1543	de_dzi(num_gp) = 0.0
1544	ENDIF
1545
1546	!
1547	!-- If wall between gridpoint 2 and gridpoint 4, then
1548	!-- Neumann boundary condition has to be applied
1549	!-- (only one case as only building beneath is possible)
1550	IF ( gp_outside_of_building(2) == 0 .AND. &
1551	gp_outside_of_building(4) == 1 ) THEN
1552	num_gp = num_gp + 1
1553	location(num_gp,1) = i * dx
1554	location(num_gp,2) = (j+1) * dy
1555	location(num_gp,3) = k * dz
1556	ei(num_gp) = e(k+1,j+1,i)
1557	dissi(num_gp) = diss(k+1,j+1,i)
1558	de_dxi(num_gp) = de_dx(k+1,j+1,i)
1559	de_dyi(num_gp) = de_dy(k+1,j+1,i)
1560	de_dzi(num_gp) = 0.0
1561	ENDIF
1562
1563	!
1564	!-- If wall between gridpoint 6 and gridpoint 8, then
1565	!-- Neumann boundary condition has to be applied
1566	!-- (only one case as only building beneath is possible)
1567	IF ( gp_outside_of_building(6) == 0 .AND. &
1568	gp_outside_of_building(8) == 1 ) THEN
1569	num_gp = num_gp + 1
1570	location(num_gp,1) = (i+1) * dx
1571	location(num_gp,2) = (j+1) * dy
1572	location(num_gp,3) = k * dz
1573	ei(num_gp) = e(k+1,j+1,i+1)
1574	dissi(num_gp) = diss(k+1,j+1,i+1)
1575	de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
1576	de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
1577	de_dzi(num_gp) = 0.0
1578	ENDIF
1579
1580	!
1581	!-- Carry out the interpolation
1582	IF ( num_gp == 1 ) THEN
1583	!
1584	!-- If only one of the gridpoints is situated outside of the
1585	!-- building, it follows that the values at the particle
1586	!-- location are the same as the gridpoint values
1587	e_int = ei(num_gp)
1588	diss_int = dissi(num_gp)
1589	de_dx_int = de_dxi(num_gp)
1590	de_dy_int = de_dyi(num_gp)
1591	de_dz_int = de_dzi(num_gp)
1592	ELSE IF ( num_gp > 1 ) THEN
1593
1594	d_sum = 0.0
1595	!
1596	!-- Evaluation of the distances between the gridpoints
1597	!-- contributing to the interpolated values, and the particle
1598	!-- location
1599	DO agp = 1, num_gp
1600	d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2 &
1601	+ ( particles(n)%y-location(agp,2) )**2 &
1602	+ ( particles(n)%z-location(agp,3) )**2
1603	d_sum = d_sum + d_gp_pl(agp)
1604	ENDDO
1605
1606	!
1607	!-- Finally the interpolation can be carried out
1608	e_int = 0.0
1609	diss_int = 0.0
1610	de_dx_int = 0.0
1611	de_dy_int = 0.0
1612	de_dz_int = 0.0
1613	DO agp = 1, num_gp
1614	e_int = e_int + ( d_sum - d_gp_pl(agp) ) * &
1615	ei(agp) / ( (num_gp-1) * d_sum )
1616	diss_int = diss_int + ( d_sum - d_gp_pl(agp) ) * &
1617	dissi(agp) / ( (num_gp-1) * d_sum )
1618	de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * &
1619	de_dxi(agp) / ( (num_gp-1) * d_sum )
1620	de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * &
1621	de_dyi(agp) / ( (num_gp-1) * d_sum )
1622	de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * &
1623	de_dzi(agp) / ( (num_gp-1) * d_sum )
1624	ENDDO
1625
1626	ENDIF
1627
1628	ENDIF
1629
1630	ENDIF
1631
1632	!
1633	!-- Vertically interpolate the horizontally averaged SGS TKE and
1634	!-- resolved-scale velocity variances and use the interpolated values
1635	!-- to calculate the coefficient fs, which is a measure of the ratio
1636	!-- of the subgrid-scale turbulent kinetic energy to the total amount
1637	!-- of turbulent kinetic energy.
1638	IF ( k == 0 ) THEN
1639	e_mean_int = hom(0,1,8,0)
1640	ELSE
1641	e_mean_int = hom(k,1,8,0) + &
1642	( hom(k+1,1,8,0) - hom(k,1,8,0) ) / &
1643	( zu(k+1) - zu(k) ) * &
1644	( particles(n)%z - zu(k) )
1645	ENDIF
1646
1647	kw = particles(n)%z / dz
1648
1649	IF ( k == 0 ) THEN
1650	aa = hom(k+1,1,30,0) * ( particles(n)%z / &
1651	( 0.5 * ( zu(k+1) - zu(k) ) ) )
1652	bb = hom(k+1,1,31,0) * ( particles(n)%z / &
1653	( 0.5 * ( zu(k+1) - zu(k) ) ) )
1654	cc = hom(kw+1,1,32,0) * ( particles(n)%z / &
1655	( 1.0 * ( zw(kw+1) - zw(kw) ) ) )
1656	ELSE
1657	aa = hom(k,1,30,0) + ( hom(k+1,1,30,0) - hom(k,1,30,0) ) * &
1658	( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
1659	bb = hom(k,1,31,0) + ( hom(k+1,1,31,0) - hom(k,1,31,0) ) * &
1660	( ( particles(n)%z - zu(k) ) / ( zu(k+1) - zu(k) ) )
1661	cc = hom(kw,1,32,0) + ( hom(kw+1,1,32,0)-hom(kw,1,32,0) ) *&
1662	( ( particles(n)%z - zw(kw) ) / ( zw(kw+1)-zw(kw) ) )
1663	ENDIF
1664
1665	vv_int = ( 1.0 / 3.0 ) * ( aa + bb + cc )
1666
1667	fs_int = ( 2.0 / 3.0 ) * e_mean_int / &
1668	( vv_int + ( 2.0 / 3.0 ) * e_mean_int )
1669
1670	!
1671	!-- Calculate the Lagrangian timescale according to the suggestion of
1672	!-- Weil et al. (2004).
1673	lagr_timescale = ( 4.0 * e_int ) / &
1674	( 3.0 * fs_int * c_0 * diss_int )
1675
1676	!
1677	!-- Calculate the next particle timestep. dt_gap is the time needed to
1678	!-- complete the current LES timestep.
1679	dt_gap = dt_3d - particles(n)%dt_sum
1680	dt_particle = MIN( dt_3d, 0.025 * lagr_timescale, dt_gap )
1681
1682	!
1683	!-- The particle timestep should not be too small in order to prevent
1684	!-- the number of particle timesteps of getting too large
1685	IF ( dt_particle < dt_min_part .AND. dt_min_part < dt_gap ) &
1686	THEN
1687	dt_particle = dt_min_part
1688	ENDIF
1689
1690	!
1691	!-- Calculate the SGS velocity components
1692	IF ( particles(n)%age == 0.0 ) THEN
1693	!
1694	!-- For new particles the SGS components are derived from the SGS
1695	!-- TKE. Limit the Gaussian random number to the interval
1696	!-- [-5.0sigma, 5.0sigma] in order to prevent the SGS velocities
1697	!-- from becoming unrealistically large.
1698	particles(n)%speed_x_sgs = SQRT( 2.0 * sgs_wfu_part * e_int ) *&
1699	( random_gauss( iran_part, 5.0 ) &
1700	- 1.0 )
1701	particles(n)%speed_y_sgs = SQRT( 2.0 * sgs_wfv_part * e_int ) *&
1702	( random_gauss( iran_part, 5.0 ) &
1703	- 1.0 )
1704	particles(n)%speed_z_sgs = SQRT( 2.0 * sgs_wfw_part * e_int ) *&
1705	( random_gauss( iran_part, 5.0 ) &
1706	- 1.0 )
1707
1708	ELSE
1709
1710	!
1711	!-- Restriction of the size of the new timestep: compared to the
1712	!-- previous timestep the increase must not exceed 200%
1713
1714	dt_particle_m = particles(n)%age - particles(n)%age_m
1715	IF ( dt_particle > 2.0 * dt_particle_m ) THEN
1716	dt_particle = 2.0 * dt_particle_m
1717	ENDIF
1718
1719	!
1720	!-- For old particles the SGS components are correlated with the
1721	!-- values from the previous timestep. Random numbers have also to
1722	!-- be limited (see above).
1723	!-- As negative values for the subgrid TKE are not allowed, the
1724	!-- change of the subgrid TKE with time cannot be smaller than
1725	!-- -e_int/dt_particle. This value is used as a lower boundary
1726	!-- value for the change of TKE
1727
1728	de_dt_min = - e_int / dt_particle
1729
1730	de_dt = ( e_int - particles(n)%e_m ) / dt_particle_m
1731
1732	IF ( de_dt < de_dt_min ) THEN
1733	de_dt = de_dt_min
1734	ENDIF
1735
1736	particles(n)%speed_x_sgs = particles(n)%speed_x_sgs - &
1737	fs_int * c_0 * diss_int * particles(n)%speed_x_sgs * &
1738	dt_particle / ( 4.0 * sgs_wfu_part * e_int ) + &
1739	( 2.0 * sgs_wfu_part * de_dt * &
1740	particles(n)%speed_x_sgs / &
1741	( 2.0 * sgs_wfu_part * e_int ) + de_dx_int &
1742	) * dt_particle / 2.0 + &
1743	SQRT( fs_int * c_0 * diss_int ) * &
1744	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
1745	SQRT( dt_particle )
1746
1747	particles(n)%speed_y_sgs = particles(n)%speed_y_sgs - &
1748	fs_int * c_0 * diss_int * particles(n)%speed_y_sgs * &
1749	dt_particle / ( 4.0 * sgs_wfv_part * e_int ) + &
1750	( 2.0 * sgs_wfv_part * de_dt * &
1751	particles(n)%speed_y_sgs / &
1752	( 2.0 * sgs_wfv_part * e_int ) + de_dy_int &
1753	) * dt_particle / 2.0 + &
1754	SQRT( fs_int * c_0 * diss_int ) * &
1755	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
1756	SQRT( dt_particle )
1757
1758	particles(n)%speed_z_sgs = particles(n)%speed_z_sgs - &
1759	fs_int * c_0 * diss_int * particles(n)%speed_z_sgs * &
1760	dt_particle / ( 4.0 * sgs_wfw_part * e_int ) + &
1761	( 2.0 * sgs_wfw_part * de_dt * &
1762	particles(n)%speed_z_sgs / &
1763	( 2.0 * sgs_wfw_part * e_int ) + de_dz_int &
1764	) * dt_particle / 2.0 + &
1765	SQRT( fs_int * c_0 * diss_int ) * &
1766	( random_gauss( iran_part, 5.0 ) - 1.0 ) * &
1767	SQRT( dt_particle )
1768
1769	ENDIF
1770
1771	u_int = u_int + particles(n)%speed_x_sgs
1772	v_int = v_int + particles(n)%speed_y_sgs
1773	w_int = w_int + particles(n)%speed_z_sgs
1774
1775	!
1776	!-- Store the SGS TKE of the current timelevel which is needed for
1777	!-- for calculating the SGS particle velocities at the next timestep
1778	particles(n)%e_m = e_int
1779
1780	ELSE
1781	!
1782	!-- If no SGS velocities are used, only the particle timestep has to
1783	!-- be set
1784	dt_particle = dt_3d
1785
1786	ENDIF
1787
1788	!
1789	!-- Remember the old age of the particle ( needed to prevent that a
1790	!-- particle crosses several PEs during one timestep and for the
1791	!-- evaluation of the subgrid particle velocity fluctuations )
1792	particles(n)%age_m = particles(n)%age
1793
1794
1795	!
1796	!-- Particle advection
1797	IF ( particle_groups(particles(n)%group)%density_ratio == 0.0 ) THEN
1798	!
1799	!-- Pure passive transport (without particle inertia)
1800	particles(n)%x = particles(n)%x + u_int * dt_particle
1801	particles(n)%y = particles(n)%y + v_int * dt_particle
1802	particles(n)%z = particles(n)%z + w_int * dt_particle
1803
1804	particles(n)%speed_x = u_int
1805	particles(n)%speed_y = v_int
1806	particles(n)%speed_z = w_int
1807
1808	ELSE
1809	!
1810	!-- Transport of particles with inertia
1811	particles(n)%x = particles(n)%x + particles(n)%speed_x * &
1812	dt_particle
1813	particles(n)%y = particles(n)%y + particles(n)%speed_y * &
1814	dt_particle
1815	particles(n)%z = particles(n)%z + particles(n)%speed_z * &
1816	dt_particle
1817
1818	!
1819	!-- Update of the particle velocity
1820	dens_ratio = particle_groups(particles(n)%group)%density_ratio
1821	IF ( cloud_droplets ) THEN
1822	exp_arg = 4.5 * dens_ratio * molecular_viscosity / &
1823	( particles(n)%radius )**2 / &
1824	( 1.0 + 0.15 * ( 2.0 * particles(n)%radius * &
1825	SQRT( ( u_int - particles(n)%speed_x )**2 + &
1826	( v_int - particles(n)%speed_y )**2 + &
1827	( w_int - particles(n)%speed_z )**2 ) / &
1828	molecular_viscosity )**0.687 &
1829	)
1830	exp_term = EXP( -exp_arg * dt_particle )
1831	ELSEIF ( use_sgs_for_particles ) THEN
1832	exp_arg = particle_groups(particles(n)%group)%exp_arg
1833	exp_term = EXP( -exp_arg * dt_particle )
1834	ELSE
1835	exp_arg = particle_groups(particles(n)%group)%exp_arg
1836	exp_term = particle_groups(particles(n)%group)%exp_term
1837	ENDIF
1838	particles(n)%speed_x = particles(n)%speed_x * exp_term + &
1839	u_int * ( 1.0 - exp_term )
1840	particles(n)%speed_y = particles(n)%speed_y * exp_term + &
1841	v_int * ( 1.0 - exp_term )
1842	particles(n)%speed_z = particles(n)%speed_z * exp_term + &
1843	( w_int - ( 1.0 - dens_ratio ) * g / exp_arg ) &
1844	* ( 1.0 - exp_term )
1845	ENDIF
1846
1847	!
1848	!-- Increment the particle age and the total time that the particle
1849	!-- has advanced within the particle timestep procedure
1850	particles(n)%age = particles(n)%age + dt_particle
1851	particles(n)%dt_sum = particles(n)%dt_sum + dt_particle
1852
1853	!
1854	!-- Check whether there is still a particle that has not yet completed
1855	!-- the total LES timestep
1856	IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8 ) THEN
1857	dt_3d_reached_l = .FALSE.
1858	ENDIF
1859
1860	ENDDO ! advection loop
1861
1862	!
1863	!-- Particle reflection from walls
1864	CALL particle_boundary_conds
1865
1866	!
1867	!-- User-defined actions after the evaluation of the new particle position
1868	CALL user_advec_particles
1869
1870	!
1871	!-- Find out, if all particles on every PE have completed the LES timestep
1872	!-- and set the switch corespondingly
1873	#if defined( __parallel )
1874	CALL MPI_ALLREDUCE( dt_3d_reached_l, dt_3d_reached, 1, MPI_LOGICAL, &
1875	MPI_LAND, comm2d, ierr )
1876	#else
1877	dt_3d_reached = dt_3d_reached_l
1878	#endif
1879
1880	CALL cpu_log( log_point_s(44), 'advec_part_advec', 'stop' )
1881
1882	!
1883	!-- Increment time since last release
1884	IF ( dt_3d_reached ) time_prel = time_prel + dt_3d
1885
1886	! WRITE ( 9, * ) '*** advec_particles: ##0.4'
1887	! CALL local_flush( 9 )
1888	! nd = 0
1889	! DO n = 1, number_of_particles
1890	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
1891	! ENDDO
1892	! IF ( nd /= deleted_particles ) THEN
1893	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
1894	! CALL local_flush( 9 )
1895	! CALL MPI_ABORT( comm2d, 9999, ierr )
1896	! ENDIF
1897
1898	!
1899	!-- If necessary, release new set of particles
1900	IF ( time_prel >= dt_prel .AND. end_time_prel > simulated_time .AND. &
1901	dt_3d_reached ) THEN
1902
1903	!
1904	!-- Check, if particle storage must be extended
1905	IF ( number_of_particles + number_of_initial_particles > &
1906	maximum_number_of_particles ) THEN
1907	IF ( netcdf_output ) THEN
1908	PRINT*, '+++ advec_particles: maximum_number_of_particles ', &
1909	'needs to be increased'
1910	PRINT*, ' but this is not allowed with ', &
1911	'NetCDF output switched on'
1912	#if defined( __parallel )
1913	CALL MPI_ABORT( comm2d, 9999, ierr )
1914	#else
1915	CALL local_stop
1916	#endif
1917	ELSE
1918	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory dt_prel'
1919	! CALL local_flush( 9 )
1920	CALL allocate_prt_memory( number_of_initial_particles )
1921	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory dt_prel'
1922	! CALL local_flush( 9 )
1923	ENDIF
1924	ENDIF
1925
1926	!
1927	!-- Check, if tail storage must be extended
1928	IF ( use_particle_tails ) THEN
1929	IF ( number_of_tails + number_of_initial_tails > &
1930	maximum_number_of_tails ) THEN
1931	IF ( netcdf_output ) THEN
1932	PRINT*, '+++ advec_particles: maximum_number_of_tails ', &
1933	'needs to be increased'
1934	PRINT*, ' but this is not allowed wi', &
1935	'th NetCDF output switched on'
1936	#if defined( __parallel )
1937	CALL MPI_ABORT( comm2d, 9999, ierr )
1938	#else
1939	CALL local_stop
1940	#endif
1941	ELSE
1942	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory dt_prel'
1943	! CALL local_flush( 9 )
1944	CALL allocate_tail_memory( number_of_initial_tails )
1945	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory dt_prel'
1946	! CALL local_flush( 9 )
1947	ENDIF
1948	ENDIF
1949	ENDIF
1950
1951	!
1952	!-- The MOD function allows for changes in the output interval with
1953	!-- restart runs.
1954	time_prel = MOD( time_prel, MAX( dt_prel, dt_3d ) )
1955	IF ( number_of_initial_particles /= 0 ) THEN
1956	is = number_of_particles+1
1957	ie = number_of_particles+number_of_initial_particles
1958	particles(is:ie) = initial_particles(1:number_of_initial_particles)
1959	!
1960	!-- Add random fluctuation to particle positions. Particles should
1961	!-- remain in the subdomain.
1962	IF ( random_start_position ) THEN
1963	DO n = is, ie
1964	IF ( psl(particles(n)%group) /= psr(particles(n)%group) ) &
1965	THEN
1966	particles(n)%x = particles(n)%x + &
1967	( random_function( iran_part ) - 0.5 ) *&
1968	pdx(particles(n)%group)
1969	IF ( particles(n)%x <= ( nxl - 0.5 ) * dx ) THEN
1970	particles(n)%x = ( nxl - 0.4999999999 ) * dx
1971	ELSEIF ( particles(n)%x >= ( nxr + 0.5 ) * dx ) THEN
1972	particles(n)%x = ( nxr + 0.4999999999 ) * dx
1973	ENDIF
1974	ENDIF
1975	IF ( pss(particles(n)%group) /= psn(particles(n)%group) ) &
1976	THEN
1977	particles(n)%y = particles(n)%y + &
1978	( random_function( iran_part ) - 0.5 ) *&
1979	pdy(particles(n)%group)
1980	IF ( particles(n)%y <= ( nys - 0.5 ) * dy ) THEN
1981	particles(n)%y = ( nys - 0.4999999999 ) * dy
1982	ELSEIF ( particles(n)%y >= ( nyn + 0.5 ) * dy ) THEN
1983	particles(n)%y = ( nyn + 0.4999999999 ) * dy
1984	ENDIF
1985	ENDIF
1986	IF ( psb(particles(n)%group) /= pst(particles(n)%group) ) &
1987	THEN
1988	particles(n)%z = particles(n)%z + &
1989	( random_function( iran_part ) - 0.5 ) *&
1990	pdz(particles(n)%group)
1991	ENDIF
1992	ENDDO
1993	ENDIF
1994
1995	!
1996	!-- Set the beginning of the new particle tails and their age
1997	IF ( use_particle_tails ) THEN
1998	DO n = is, ie
1999	!
2000	!-- New particles which should have a tail, already have got a
2001	!-- provisional tail id unequal zero (see init_particles)
2002	IF ( particles(n)%tail_id /= 0 ) THEN
2003	number_of_tails = number_of_tails + 1
2004	nn = number_of_tails
2005	particles(n)%tail_id = nn ! set the final tail id
2006	particle_tail_coordinates(1,1,nn) = particles(n)%x
2007	particle_tail_coordinates(1,2,nn) = particles(n)%y
2008	particle_tail_coordinates(1,3,nn) = particles(n)%z
2009	particle_tail_coordinates(1,4,nn) = particles(n)%color
2010	particles(n)%tailpoints = 1
2011	IF ( minimum_tailpoint_distance /= 0.0 ) THEN
2012	particle_tail_coordinates(2,1,nn) = particles(n)%x
2013	particle_tail_coordinates(2,2,nn) = particles(n)%y
2014	particle_tail_coordinates(2,3,nn) = particles(n)%z
2015	particle_tail_coordinates(2,4,nn) = particles(n)%color
2016	particle_tail_coordinates(1:2,5,nn) = 0.0
2017	particles(n)%tailpoints = 2
2018	ENDIF
2019	ENDIF
2020	ENDDO
2021	ENDIF
2022	! WRITE ( 9, * ) '*** advec_particles: after setting the beginning of new tails'
2023	! CALL local_flush( 9 )
2024
2025	number_of_particles = number_of_particles + &
2026	number_of_initial_particles
2027	ENDIF
2028
2029	ENDIF
2030
2031	! WRITE ( 9, * ) '*** advec_particles: ##0.5'
2032	! CALL local_flush( 9 )
2033	! nd = 0
2034	! DO n = 1, number_of_particles
2035	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
2036	! ENDDO
2037	! IF ( nd /= deleted_particles ) THEN
2038	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
2039	! CALL local_flush( 9 )
2040	! CALL MPI_ABORT( comm2d, 9999, ierr )
2041	! ENDIF
2042
2043	! IF ( number_of_particles /= number_of_tails ) THEN
2044	! WRITE (9,*) '--- advec_particles: #2'
2045	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
2046	! CALL local_flush( 9 )
2047	! ENDIF
2048	! DO n = 1, number_of_particles
2049	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
2050	! THEN
2051	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
2052	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
2053	! CALL local_flush( 9 )
2054	! CALL MPI_ABORT( comm2d, 9999, ierr )
2055	! ENDIF
2056	! ENDDO
2057
2058	#if defined( __parallel )
2059	!
2060	!-- As soon as one particle has moved beyond the boundary of the domain, it
2061	!-- is included in the relevant transfer arrays and marked for subsequent
2062	!-- deletion on this PE.
2063	!-- First run for crossings in x direction. Find out first the number of
2064	!-- particles to be transferred and allocate temporary arrays needed to store
2065	!-- them.
2066	!-- For a one-dimensional decomposition along y, no transfer is necessary,
2067	!-- because the particle remains on the PE.
2068	trlp_count = 0
2069	trlpt_count = 0
2070	trrp_count = 0
2071	trrpt_count = 0
2072	IF ( pdims(1) /= 1 ) THEN
2073	!
2074	!-- First calculate the storage necessary for sending and receiving the
2075	!-- data
2076	DO n = 1, number_of_particles
2077	i = ( particles(n)%x + 0.5 * dx ) * ddx
2078	!
2079	!-- Above calculation does not work for indices less than zero
2080	IF ( particles(n)%x < -0.5 * dx ) i = -1
2081
2082	IF ( i < nxl ) THEN
2083	trlp_count = trlp_count + 1
2084	IF ( particles(n)%tail_id /= 0 ) trlpt_count = trlpt_count + 1
2085	ELSEIF ( i > nxr ) THEN
2086	trrp_count = trrp_count + 1
2087	IF ( particles(n)%tail_id /= 0 ) trrpt_count = trrpt_count + 1
2088	ENDIF
2089	ENDDO
2090	IF ( trlp_count == 0 ) trlp_count = 1
2091	IF ( trlpt_count == 0 ) trlpt_count = 1
2092	IF ( trrp_count == 0 ) trrp_count = 1
2093	IF ( trrpt_count == 0 ) trrpt_count = 1
2094
2095	ALLOCATE( trlp(trlp_count), trrp(trrp_count) )
2096
2097	trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
2098	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
2099	0.0, 0, 0, 0, 0 )
2100	trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
2101	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
2102	0.0, 0, 0, 0, 0 )
2103
2104	IF ( use_particle_tails ) THEN
2105	ALLOCATE( trlpt(maximum_number_of_tailpoints,5,trlpt_count), &
2106	trrpt(maximum_number_of_tailpoints,5,trrpt_count) )
2107	tlength = maximum_number_of_tailpoints * 5
2108	ENDIF
2109
2110	trlp_count = 0
2111	trlpt_count = 0
2112	trrp_count = 0
2113	trrpt_count = 0
2114
2115	ENDIF
2116
2117	! WRITE ( 9, * ) '*** advec_particles: ##1'
2118	! CALL local_flush( 9 )
2119	! nd = 0
2120	! DO n = 1, number_of_particles
2121	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
2122	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
2123	! THEN
2124	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
2125	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
2126	! CALL local_flush( 9 )
2127	! CALL MPI_ABORT( comm2d, 9999, ierr )
2128	! ENDIF
2129	! ENDDO
2130	! IF ( nd /= deleted_particles ) THEN
2131	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
2132	! CALL local_flush( 9 )
2133	! CALL MPI_ABORT( comm2d, 9999, ierr )
2134	! ENDIF
2135
2136	DO n = 1, number_of_particles
2137
2138	nn = particles(n)%tail_id
2139
2140	i = ( particles(n)%x + 0.5 * dx ) * ddx
2141	!
2142	!-- Above calculation does not work for indices less than zero
2143	IF ( particles(n)%x < - 0.5 * dx ) i = -1
2144
2145	IF ( i < nxl ) THEN
2146	IF ( i < 0 ) THEN
2147	!
2148	!-- Apply boundary condition along x
2149	IF ( ibc_par_lr == 0 ) THEN
2150	!
2151	!-- Cyclic condition
2152	IF ( pdims(1) == 1 ) THEN
2153	particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
2154	particles(n)%origin_x = ( nx + 1 ) * dx + &
2155	particles(n)%origin_x
2156	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2157	i = particles(n)%tailpoints
2158	particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx &
2159	+ particle_tail_coordinates(1:i,1,nn)
2160	ENDIF
2161	ELSE
2162	trlp_count = trlp_count + 1
2163	trlp(trlp_count) = particles(n)
2164	trlp(trlp_count)%x = ( nx + 1 ) * dx + trlp(trlp_count)%x
2165	trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x + &
2166	( nx + 1 ) * dx
2167	particle_mask(n) = .FALSE.
2168	deleted_particles = deleted_particles + 1
2169
2170	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2171	trlpt_count = trlpt_count + 1
2172	trlpt(:,:,trlpt_count) = &
2173	particle_tail_coordinates(:,:,nn)
2174	trlpt(:,1,trlpt_count) = ( nx + 1 ) * dx + &
2175	trlpt(:,1,trlpt_count)
2176	tail_mask(nn) = .FALSE.
2177	deleted_tails = deleted_tails + 1
2178	ENDIF
2179	ENDIF
2180
2181	ELSEIF ( ibc_par_lr == 1 ) THEN
2182	!
2183	!-- Particle absorption
2184	particle_mask(n) = .FALSE.
2185	deleted_particles = deleted_particles + 1
2186	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2187	tail_mask(nn) = .FALSE.
2188	deleted_tails = deleted_tails + 1
2189	ENDIF
2190
2191	ELSEIF ( ibc_par_lr == 2 ) THEN
2192	!
2193	!-- Particle reflection
2194	particles(n)%x = -particles(n)%x
2195	particles(n)%speed_x = -particles(n)%speed_x
2196
2197	ENDIF
2198	ELSE
2199	!
2200	!-- Store particle data in the transfer array, which will be send
2201	!-- to the neighbouring PE
2202	trlp_count = trlp_count + 1
2203	trlp(trlp_count) = particles(n)
2204	particle_mask(n) = .FALSE.
2205	deleted_particles = deleted_particles + 1
2206
2207	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2208	trlpt_count = trlpt_count + 1
2209	trlpt(:,:,trlpt_count) = particle_tail_coordinates(:,:,nn)
2210	tail_mask(nn) = .FALSE.
2211	deleted_tails = deleted_tails + 1
2212	ENDIF
2213	ENDIF
2214
2215	ELSEIF ( i > nxr ) THEN
2216	IF ( i > nx ) THEN
2217	!
2218	!-- Apply boundary condition along x
2219	IF ( ibc_par_lr == 0 ) THEN
2220	!
2221	!-- Cyclic condition
2222	IF ( pdims(1) == 1 ) THEN
2223	particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
2224	particles(n)%origin_x = particles(n)%origin_x - &
2225	( nx + 1 ) * dx
2226	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2227	i = particles(n)%tailpoints
2228	particle_tail_coordinates(1:i,1,nn) = - ( nx+1 ) * dx &
2229	+ particle_tail_coordinates(1:i,1,nn)
2230	ENDIF
2231	ELSE
2232	trrp_count = trrp_count + 1
2233	trrp(trrp_count) = particles(n)
2234	trrp(trrp_count)%x = trrp(trrp_count)%x - ( nx + 1 ) * dx
2235	trrp(trrp_count)%origin_x = trrp(trrp_count)%origin_x - &
2236	( nx + 1 ) * dx
2237	particle_mask(n) = .FALSE.
2238	deleted_particles = deleted_particles + 1
2239
2240	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2241	trrpt_count = trrpt_count + 1
2242	trrpt(:,:,trrpt_count) = &
2243	particle_tail_coordinates(:,:,nn)
2244	trrpt(:,1,trrpt_count) = trrpt(:,1,trrpt_count) - &
2245	( nx + 1 ) * dx
2246	tail_mask(nn) = .FALSE.
2247	deleted_tails = deleted_tails + 1
2248	ENDIF
2249	ENDIF
2250
2251	ELSEIF ( ibc_par_lr == 1 ) THEN
2252	!
2253	!-- Particle absorption
2254	particle_mask(n) = .FALSE.
2255	deleted_particles = deleted_particles + 1
2256	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2257	tail_mask(nn) = .FALSE.
2258	deleted_tails = deleted_tails + 1
2259	ENDIF
2260
2261	ELSEIF ( ibc_par_lr == 2 ) THEN
2262	!
2263	!-- Particle reflection
2264	particles(n)%x = 2 * ( nx * dx ) - particles(n)%x
2265	particles(n)%speed_x = -particles(n)%speed_x
2266
2267	ENDIF
2268	ELSE
2269	!
2270	!-- Store particle data in the transfer array, which will be send
2271	!-- to the neighbouring PE
2272	trrp_count = trrp_count + 1
2273	trrp(trrp_count) = particles(n)
2274	particle_mask(n) = .FALSE.
2275	deleted_particles = deleted_particles + 1
2276
2277	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2278	trrpt_count = trrpt_count + 1
2279	trrpt(:,:,trrpt_count) = particle_tail_coordinates(:,:,nn)
2280	tail_mask(nn) = .FALSE.
2281	deleted_tails = deleted_tails + 1
2282	ENDIF
2283	ENDIF
2284
2285	ENDIF
2286	ENDDO
2287
2288	! WRITE ( 9, * ) '*** advec_particles: ##2'
2289	! CALL local_flush( 9 )
2290	! nd = 0
2291	! DO n = 1, number_of_particles
2292	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
2293	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
2294	! THEN
2295	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
2296	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
2297	! CALL local_flush( 9 )
2298	! CALL MPI_ABORT( comm2d, 9999, ierr )
2299	! ENDIF
2300	! ENDDO
2301	! IF ( nd /= deleted_particles ) THEN
2302	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
2303	! CALL local_flush( 9 )
2304	! CALL MPI_ABORT( comm2d, 9999, ierr )
2305	! ENDIF
2306
2307	!
2308	!-- Send left boundary, receive right boundary (but first exchange how many
2309	!-- and check, if particle storage must be extended)
2310	IF ( pdims(1) /= 1 ) THEN
2311
2312	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'start' )
2313	CALL MPI_SENDRECV( trlp_count, 1, MPI_INTEGER, pleft, 0, &
2314	trrp_count_recv, 1, MPI_INTEGER, pright, 0, &
2315	comm2d, status, ierr )
2316
2317	IF ( number_of_particles + trrp_count_recv > &
2318	maximum_number_of_particles ) &
2319	THEN
2320	IF ( netcdf_output ) THEN
2321	PRINT*, '+++ advec_particles: maximum_number_of_particles ', &
2322	'needs to be increased'
2323	PRINT*, ' but this is not allowed with ', &
2324	'NetCDF output switched on'
2325	CALL MPI_ABORT( comm2d, 9999, ierr )
2326	ELSE
2327	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trrp'
2328	! CALL local_flush( 9 )
2329	CALL allocate_prt_memory( trrp_count_recv )
2330	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trrp'
2331	! CALL local_flush( 9 )
2332	ENDIF
2333	ENDIF
2334
2335	CALL MPI_SENDRECV( trlp(1)%age, trlp_count, mpi_particle_type, &
2336	pleft, 1, particles(number_of_particles+1)%age, &
2337	trrp_count_recv, mpi_particle_type, pright, 1, &
2338	comm2d, status, ierr )
2339
2340	IF ( use_particle_tails ) THEN
2341
2342	CALL MPI_SENDRECV( trlpt_count, 1, MPI_INTEGER, pleft, 0, &
2343	trrpt_count_recv, 1, MPI_INTEGER, pright, 0, &
2344	comm2d, status, ierr )
2345
2346	IF ( number_of_tails+trrpt_count_recv > maximum_number_of_tails ) &
2347	THEN
2348	IF ( netcdf_output ) THEN
2349	PRINT*, '+++ advec_particles: maximum_number_of_tails ', &
2350	'needs to be increased'
2351	PRINT*, ' but this is not allowed wi', &
2352	'th NetCDF output switched on'
2353	CALL MPI_ABORT( comm2d, 9999, ierr )
2354	ELSE
2355	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trrpt'
2356	! CALL local_flush( 9 )
2357	CALL allocate_tail_memory( trrpt_count_recv )
2358	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trrpt'
2359	! CALL local_flush( 9 )
2360	ENDIF
2361	ENDIF
2362
2363	CALL MPI_SENDRECV( trlpt(1,1,1), trlpt_count*tlength, MPI_REAL, &
2364	pleft, 1, &
2365	particle_tail_coordinates(1,1,number_of_tails+1), &
2366	trrpt_count_recv*tlength, MPI_REAL, pright, 1, &
2367	comm2d, status, ierr )
2368	!
2369	!-- Update the tail ids for the transferred particles
2370	nn = number_of_tails
2371	DO n = number_of_particles+1, number_of_particles+trrp_count_recv
2372	IF ( particles(n)%tail_id /= 0 ) THEN
2373	nn = nn + 1
2374	particles(n)%tail_id = nn
2375	ENDIF
2376	ENDDO
2377
2378	ENDIF
2379
2380	number_of_particles = number_of_particles + trrp_count_recv
2381	number_of_tails = number_of_tails + trrpt_count_recv
2382	! IF ( number_of_particles /= number_of_tails ) THEN
2383	! WRITE (9,*) '--- advec_particles: #3'
2384	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
2385	! CALL local_flush( 9 )
2386	! ENDIF
2387
2388	!
2389	!-- Send right boundary, receive left boundary
2390	CALL MPI_SENDRECV( trrp_count, 1, MPI_INTEGER, pright, 0, &
2391	trlp_count_recv, 1, MPI_INTEGER, pleft, 0, &
2392	comm2d, status, ierr )
2393
2394	IF ( number_of_particles + trlp_count_recv > &
2395	maximum_number_of_particles ) &
2396	THEN
2397	IF ( netcdf_output ) THEN
2398	PRINT*, '+++ advec_particles: maximum_number_of_particles ', &
2399	'needs to be increased'
2400	PRINT*, ' but this is not allowed with ', &
2401	'NetCDF output switched on'
2402	CALL MPI_ABORT( comm2d, 9999, ierr )
2403	ELSE
2404	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trlp'
2405	! CALL local_flush( 9 )
2406	CALL allocate_prt_memory( trlp_count_recv )
2407	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trlp'
2408	! CALL local_flush( 9 )
2409	ENDIF
2410	ENDIF
2411
2412	CALL MPI_SENDRECV( trrp(1)%age, trrp_count, mpi_particle_type, &
2413	pright, 1, particles(number_of_particles+1)%age, &
2414	trlp_count_recv, mpi_particle_type, pleft, 1, &
2415	comm2d, status, ierr )
2416
2417	IF ( use_particle_tails ) THEN
2418
2419	CALL MPI_SENDRECV( trrpt_count, 1, MPI_INTEGER, pright, 0, &
2420	trlpt_count_recv, 1, MPI_INTEGER, pleft, 0, &
2421	comm2d, status, ierr )
2422
2423	IF ( number_of_tails+trlpt_count_recv > maximum_number_of_tails ) &
2424	THEN
2425	IF ( netcdf_output ) THEN
2426	PRINT*, '+++ advec_particles: maximum_number_of_tails ', &
2427	'needs to be increased'
2428	PRINT*, ' but this is not allowed wi', &
2429	'th NetCDF output switched on'
2430	CALL MPI_ABORT( comm2d, 9999, ierr )
2431	ELSE
2432	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trlpt'
2433	! CALL local_flush( 9 )
2434	CALL allocate_tail_memory( trlpt_count_recv )
2435	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trlpt'
2436	! CALL local_flush( 9 )
2437	ENDIF
2438	ENDIF
2439
2440	CALL MPI_SENDRECV( trrpt(1,1,1), trrpt_count*tlength, MPI_REAL, &
2441	pright, 1, &
2442	particle_tail_coordinates(1,1,number_of_tails+1), &
2443	trlpt_count_recv*tlength, MPI_REAL, pleft, 1, &
2444	comm2d, status, ierr )
2445	!
2446	!-- Update the tail ids for the transferred particles
2447	nn = number_of_tails
2448	DO n = number_of_particles+1, number_of_particles+trlp_count_recv
2449	IF ( particles(n)%tail_id /= 0 ) THEN
2450	nn = nn + 1
2451	particles(n)%tail_id = nn
2452	ENDIF
2453	ENDDO
2454
2455	ENDIF
2456
2457	number_of_particles = number_of_particles + trlp_count_recv
2458	number_of_tails = number_of_tails + trlpt_count_recv
2459	! IF ( number_of_particles /= number_of_tails ) THEN
2460	! WRITE (9,*) '--- advec_particles: #4'
2461	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
2462	! CALL local_flush( 9 )
2463	! ENDIF
2464
2465	IF ( use_particle_tails ) THEN
2466	DEALLOCATE( trlpt, trrpt )
2467	ENDIF
2468	DEALLOCATE( trlp, trrp )
2469
2470	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'pause' )
2471
2472	ENDIF
2473
2474	! WRITE ( 9, * ) '*** advec_particles: ##3'
2475	! CALL local_flush( 9 )
2476	! nd = 0
2477	! DO n = 1, number_of_particles
2478	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
2479	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
2480	! THEN
2481	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
2482	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
2483	! CALL local_flush( 9 )
2484	! CALL MPI_ABORT( comm2d, 9999, ierr )
2485	! ENDIF
2486	! ENDDO
2487	! IF ( nd /= deleted_particles ) THEN
2488	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
2489	! CALL local_flush( 9 )
2490	! CALL MPI_ABORT( comm2d, 9999, ierr )
2491	! ENDIF
2492
2493	!
2494	!-- Check whether particles have crossed the boundaries in y direction. Note
2495	!-- that this case can also apply to particles that have just been received
2496	!-- from the adjacent right or left PE.
2497	!-- Find out first the number of particles to be transferred and allocate
2498	!-- temporary arrays needed to store them.
2499	!-- For a one-dimensional decomposition along x, no transfer is necessary,
2500	!-- because the particle remains on the PE.
2501	trsp_count = 0
2502	trspt_count = 0
2503	trnp_count = 0
2504	trnpt_count = 0
2505	IF ( pdims(2) /= 1 ) THEN
2506	!
2507	!-- First calculate the storage necessary for sending and receiving the
2508	!-- data
2509	DO n = 1, number_of_particles
2510	IF ( particle_mask(n) ) THEN
2511	j = ( particles(n)%y + 0.5 * dy ) * ddy
2512	!
2513	!-- Above calculation does not work for indices less than zero
2514	IF ( particles(n)%y < -0.5 * dy ) j = -1
2515
2516	IF ( j < nys ) THEN
2517	trsp_count = trsp_count + 1
2518	IF ( particles(n)%tail_id /= 0 ) trspt_count = trspt_count+1
2519	ELSEIF ( j > nyn ) THEN
2520	trnp_count = trnp_count + 1
2521	IF ( particles(n)%tail_id /= 0 ) trnpt_count = trnpt_count+1
2522	ENDIF
2523	ENDIF
2524	ENDDO
2525	IF ( trsp_count == 0 ) trsp_count = 1
2526	IF ( trspt_count == 0 ) trspt_count = 1
2527	IF ( trnp_count == 0 ) trnp_count = 1
2528	IF ( trnpt_count == 0 ) trnpt_count = 1
2529
2530	ALLOCATE( trsp(trsp_count), trnp(trnp_count) )
2531
2532	trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
2533	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
2534	0.0, 0, 0, 0, 0 )
2535	trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
2536	0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
2537	0.0, 0, 0, 0, 0 )
2538
2539	IF ( use_particle_tails ) THEN
2540	ALLOCATE( trspt(maximum_number_of_tailpoints,5,trspt_count), &
2541	trnpt(maximum_number_of_tailpoints,5,trnpt_count) )
2542	tlength = maximum_number_of_tailpoints * 5
2543	ENDIF
2544
2545	trsp_count = 0
2546	trspt_count = 0
2547	trnp_count = 0
2548	trnpt_count = 0
2549
2550	ENDIF
2551
2552	! WRITE ( 9, * ) '*** advec_particles: ##4'
2553	! CALL local_flush( 9 )
2554	! nd = 0
2555	! DO n = 1, number_of_particles
2556	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
2557	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
2558	! THEN
2559	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
2560	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
2561	! CALL local_flush( 9 )
2562	! CALL MPI_ABORT( comm2d, 9999, ierr )
2563	! ENDIF
2564	! ENDDO
2565	! IF ( nd /= deleted_particles ) THEN
2566	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
2567	! CALL local_flush( 9 )
2568	! CALL MPI_ABORT( comm2d, 9999, ierr )
2569	! ENDIF
2570
2571	DO n = 1, number_of_particles
2572
2573	nn = particles(n)%tail_id
2574	!
2575	!-- Only those particles that have not been marked as 'deleted' may be
2576	!-- moved.
2577	IF ( particle_mask(n) ) THEN
2578	j = ( particles(n)%y + 0.5 * dy ) * ddy
2579	!
2580	!-- Above calculation does not work for indices less than zero
2581	IF ( particles(n)%y < -0.5 * dy ) j = -1
2582
2583	IF ( j < nys ) THEN
2584	IF ( j < 0 ) THEN
2585	!
2586	!-- Apply boundary condition along y
2587	IF ( ibc_par_ns == 0 ) THEN
2588	!
2589	!-- Cyclic condition
2590	IF ( pdims(2) == 1 ) THEN
2591	particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
2592	particles(n)%origin_y = ( ny + 1 ) * dy + &
2593	particles(n)%origin_y
2594	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2595	i = particles(n)%tailpoints
2596	particle_tail_coordinates(1:i,2,nn) = ( ny+1 ) * dy&
2597	+ particle_tail_coordinates(1:i,2,nn)
2598	ENDIF
2599	ELSE
2600	trsp_count = trsp_count + 1
2601	trsp(trsp_count) = particles(n)
2602	trsp(trsp_count)%y = ( ny + 1 ) * dy + &
2603	trsp(trsp_count)%y
2604	trsp(trsp_count)%origin_y = trsp(trsp_count)%origin_y &
2605	+ ( ny + 1 ) * dy
2606	particle_mask(n) = .FALSE.
2607	deleted_particles = deleted_particles + 1
2608
2609	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2610	trspt_count = trspt_count + 1
2611	trspt(:,:,trspt_count) = &
2612	particle_tail_coordinates(:,:,nn)
2613	trspt(:,2,trspt_count) = ( ny + 1 ) * dy + &
2614	trspt(:,2,trspt_count)
2615	tail_mask(nn) = .FALSE.
2616	deleted_tails = deleted_tails + 1
2617	ENDIF
2618	ENDIF
2619
2620	ELSEIF ( ibc_par_ns == 1 ) THEN
2621	!
2622	!-- Particle absorption
2623	particle_mask(n) = .FALSE.
2624	deleted_particles = deleted_particles + 1
2625	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2626	tail_mask(nn) = .FALSE.
2627	deleted_tails = deleted_tails + 1
2628	ENDIF
2629
2630	ELSEIF ( ibc_par_ns == 2 ) THEN
2631	!
2632	!-- Particle reflection
2633	particles(n)%y = -particles(n)%y
2634	particles(n)%speed_y = -particles(n)%speed_y
2635
2636	ENDIF
2637	ELSE
2638	!
2639	!-- Store particle data in the transfer array, which will be send
2640	!-- to the neighbouring PE
2641	trsp_count = trsp_count + 1
2642	trsp(trsp_count) = particles(n)
2643	particle_mask(n) = .FALSE.
2644	deleted_particles = deleted_particles + 1
2645
2646	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2647	trspt_count = trspt_count + 1
2648	trspt(:,:,trspt_count) = particle_tail_coordinates(:,:,nn)
2649	tail_mask(nn) = .FALSE.
2650	deleted_tails = deleted_tails + 1
2651	ENDIF
2652	ENDIF
2653
2654	ELSEIF ( j > nyn ) THEN
2655	IF ( j > ny ) THEN
2656	!
2657	!-- Apply boundary condition along x
2658	IF ( ibc_par_ns == 0 ) THEN
2659	!
2660	!-- Cyclic condition
2661	IF ( pdims(2) == 1 ) THEN
2662	particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
2663	particles(n)%origin_y = particles(n)%origin_y - &
2664	( ny + 1 ) * dy
2665	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2666	i = particles(n)%tailpoints
2667	particle_tail_coordinates(1:i,2,nn) = - (ny+1) * dy&
2668	+ particle_tail_coordinates(1:i,2,nn)
2669	ENDIF
2670	ELSE
2671	trnp_count = trnp_count + 1
2672	trnp(trnp_count) = particles(n)
2673	trnp(trnp_count)%y = trnp(trnp_count)%y - &
2674	( ny + 1 ) * dy
2675	trnp(trnp_count)%origin_y = trnp(trnp_count)%origin_y &
2676	- ( ny + 1 ) * dy
2677	particle_mask(n) = .FALSE.
2678	deleted_particles = deleted_particles + 1
2679
2680	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2681	trnpt_count = trnpt_count + 1
2682	trnpt(:,:,trnpt_count) = &
2683	particle_tail_coordinates(:,:,nn)
2684	trnpt(:,2,trnpt_count) = trnpt(:,2,trnpt_count) - &
2685	( ny + 1 ) * dy
2686	tail_mask(nn) = .FALSE.
2687	deleted_tails = deleted_tails + 1
2688	ENDIF
2689	ENDIF
2690
2691	ELSEIF ( ibc_par_ns == 1 ) THEN
2692	!
2693	!-- Particle absorption
2694	particle_mask(n) = .FALSE.
2695	deleted_particles = deleted_particles + 1
2696	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2697	tail_mask(nn) = .FALSE.
2698	deleted_tails = deleted_tails + 1
2699	ENDIF
2700
2701	ELSEIF ( ibc_par_ns == 2 ) THEN
2702	!
2703	!-- Particle reflection
2704	particles(n)%y = 2 * ( ny * dy ) - particles(n)%y
2705	particles(n)%speed_y = -particles(n)%speed_y
2706
2707	ENDIF
2708	ELSE
2709	!
2710	!-- Store particle data in the transfer array, which will be send
2711	!-- to the neighbouring PE
2712	trnp_count = trnp_count + 1
2713	trnp(trnp_count) = particles(n)
2714	particle_mask(n) = .FALSE.
2715	deleted_particles = deleted_particles + 1
2716
2717	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2718	trnpt_count = trnpt_count + 1
2719	trnpt(:,:,trnpt_count) = particle_tail_coordinates(:,:,nn)
2720	tail_mask(nn) = .FALSE.
2721	deleted_tails = deleted_tails + 1
2722	ENDIF
2723	ENDIF
2724
2725	ENDIF
2726	ENDIF
2727	ENDDO
2728
2729	! WRITE ( 9, * ) '*** advec_particles: ##5'
2730	! CALL local_flush( 9 )
2731	! nd = 0
2732	! DO n = 1, number_of_particles
2733	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
2734	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
2735	! THEN
2736	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
2737	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
2738	! CALL local_flush( 9 )
2739	! CALL MPI_ABORT( comm2d, 9999, ierr )
2740	! ENDIF
2741	! ENDDO
2742	! IF ( nd /= deleted_particles ) THEN
2743	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
2744	! CALL local_flush( 9 )
2745	! CALL MPI_ABORT( comm2d, 9999, ierr )
2746	! ENDIF
2747
2748	!
2749	!-- Send front boundary, receive back boundary (but first exchange how many
2750	!-- and check, if particle storage must be extended)
2751	IF ( pdims(2) /= 1 ) THEN
2752
2753	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'continue' )
2754	CALL MPI_SENDRECV( trsp_count, 1, MPI_INTEGER, psouth, 0, &
2755	trnp_count_recv, 1, MPI_INTEGER, pnorth, 0, &
2756	comm2d, status, ierr )
2757
2758	IF ( number_of_particles + trnp_count_recv > &
2759	maximum_number_of_particles ) &
2760	THEN
2761	IF ( netcdf_output ) THEN
2762	PRINT*, '+++ advec_particles: maximum_number_of_particles ', &
2763	'needs to be increased'
2764	PRINT*, ' but this is not allowed with ', &
2765	'NetCDF output switched on'
2766	CALL MPI_ABORT( comm2d, 9999, ierr )
2767	ELSE
2768	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trnp'
2769	! CALL local_flush( 9 )
2770	CALL allocate_prt_memory( trnp_count_recv )
2771	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trnp'
2772	! CALL local_flush( 9 )
2773	ENDIF
2774	ENDIF
2775
2776	CALL MPI_SENDRECV( trsp(1)%age, trsp_count, mpi_particle_type, &
2777	psouth, 1, particles(number_of_particles+1)%age, &
2778	trnp_count_recv, mpi_particle_type, pnorth, 1, &
2779	comm2d, status, ierr )
2780
2781	IF ( use_particle_tails ) THEN
2782
2783	CALL MPI_SENDRECV( trspt_count, 1, MPI_INTEGER, psouth, 0, &
2784	trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, &
2785	comm2d, status, ierr )
2786
2787	IF ( number_of_tails+trnpt_count_recv > maximum_number_of_tails ) &
2788	THEN
2789	IF ( netcdf_output ) THEN
2790	PRINT*, '+++ advec_particles: maximum_number_of_tails ', &
2791	'needs to be increased'
2792	PRINT*, ' but this is not allowed wi', &
2793	'th NetCDF output switched on'
2794	CALL MPI_ABORT( comm2d, 9999, ierr )
2795	ELSE
2796	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trnpt'
2797	! CALL local_flush( 9 )
2798	CALL allocate_tail_memory( trnpt_count_recv )
2799	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trnpt'
2800	! CALL local_flush( 9 )
2801	ENDIF
2802	ENDIF
2803
2804	CALL MPI_SENDRECV( trspt(1,1,1), trspt_count*tlength, MPI_REAL, &
2805	psouth, 1, &
2806	particle_tail_coordinates(1,1,number_of_tails+1), &
2807	trnpt_count_recv*tlength, MPI_REAL, pnorth, 1, &
2808	comm2d, status, ierr )
2809	!
2810	!-- Update the tail ids for the transferred particles
2811	nn = number_of_tails
2812	DO n = number_of_particles+1, number_of_particles+trnp_count_recv
2813	IF ( particles(n)%tail_id /= 0 ) THEN
2814	nn = nn + 1
2815	particles(n)%tail_id = nn
2816	ENDIF
2817	ENDDO
2818
2819	ENDIF
2820
2821	number_of_particles = number_of_particles + trnp_count_recv
2822	number_of_tails = number_of_tails + trnpt_count_recv
2823	! IF ( number_of_particles /= number_of_tails ) THEN
2824	! WRITE (9,*) '--- advec_particles: #5'
2825	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
2826	! CALL local_flush( 9 )
2827	! ENDIF
2828
2829	!
2830	!-- Send back boundary, receive front boundary
2831	CALL MPI_SENDRECV( trnp_count, 1, MPI_INTEGER, pnorth, 0, &
2832	trsp_count_recv, 1, MPI_INTEGER, psouth, 0, &
2833	comm2d, status, ierr )
2834
2835	IF ( number_of_particles + trsp_count_recv > &
2836	maximum_number_of_particles ) &
2837	THEN
2838	IF ( netcdf_output ) THEN
2839	PRINT*, '+++ advec_particles: maximum_number_of_particles ', &
2840	'needs to be increased'
2841	PRINT*, ' but this is not allowed with ', &
2842	'NetCDF output switched on'
2843	CALL MPI_ABORT( comm2d, 9999, ierr )
2844	ELSE
2845	! WRITE ( 9, * ) '*** advec_particles: before allocate_prt_memory trsp'
2846	! CALL local_flush( 9 )
2847	CALL allocate_prt_memory( trsp_count_recv )
2848	! WRITE ( 9, * ) '*** advec_particles: after allocate_prt_memory trsp'
2849	! CALL local_flush( 9 )
2850	ENDIF
2851	ENDIF
2852
2853	CALL MPI_SENDRECV( trnp(1)%age, trnp_count, mpi_particle_type, &
2854	pnorth, 1, particles(number_of_particles+1)%age, &
2855	trsp_count_recv, mpi_particle_type, psouth, 1, &
2856	comm2d, status, ierr )
2857
2858	IF ( use_particle_tails ) THEN
2859
2860	CALL MPI_SENDRECV( trnpt_count, 1, MPI_INTEGER, pnorth, 0, &
2861	trspt_count_recv, 1, MPI_INTEGER, psouth, 0, &
2862	comm2d, status, ierr )
2863
2864	IF ( number_of_tails+trspt_count_recv > maximum_number_of_tails ) &
2865	THEN
2866	IF ( netcdf_output ) THEN
2867	PRINT*, '+++ advec_particles: maximum_number_of_tails ', &
2868	'needs to be increased'
2869	PRINT*, ' but this is not allowed wi', &
2870	'th NetCDF output switched on'
2871	CALL MPI_ABORT( comm2d, 9999, ierr )
2872	ELSE
2873	! WRITE ( 9, * ) '*** advec_particles: before allocate_tail_memory trspt'
2874	! CALL local_flush( 9 )
2875	CALL allocate_tail_memory( trspt_count_recv )
2876	! WRITE ( 9, * ) '*** advec_particles: after allocate_tail_memory trspt'
2877	! CALL local_flush( 9 )
2878	ENDIF
2879	ENDIF
2880
2881	CALL MPI_SENDRECV( trnpt(1,1,1), trnpt_count*tlength, MPI_REAL, &
2882	pnorth, 1, &
2883	particle_tail_coordinates(1,1,number_of_tails+1), &
2884	trspt_count_recv*tlength, MPI_REAL, psouth, 1, &
2885	comm2d, status, ierr )
2886	!
2887	!-- Update the tail ids for the transferred particles
2888	nn = number_of_tails
2889	DO n = number_of_particles+1, number_of_particles+trsp_count_recv
2890	IF ( particles(n)%tail_id /= 0 ) THEN
2891	nn = nn + 1
2892	particles(n)%tail_id = nn
2893	ENDIF
2894	ENDDO
2895
2896	ENDIF
2897
2898	number_of_particles = number_of_particles + trsp_count_recv
2899	number_of_tails = number_of_tails + trspt_count_recv
2900	! IF ( number_of_particles /= number_of_tails ) THEN
2901	! WRITE (9,*) '--- advec_particles: #6'
2902	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
2903	! CALL local_flush( 9 )
2904	! ENDIF
2905
2906	IF ( use_particle_tails ) THEN
2907	DEALLOCATE( trspt, trnpt )
2908	ENDIF
2909	DEALLOCATE( trsp, trnp )
2910
2911	CALL cpu_log( log_point_s(23), 'sendrcv_particles', 'stop' )
2912
2913	ENDIF
2914
2915	! WRITE ( 9, * ) '*** advec_particles: ##6'
2916	! CALL local_flush( 9 )
2917	! nd = 0
2918	! DO n = 1, number_of_particles
2919	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
2920	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
2921	! THEN
2922	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
2923	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
2924	! CALL local_flush( 9 )
2925	! CALL MPI_ABORT( comm2d, 9999, ierr )
2926	! ENDIF
2927	! ENDDO
2928	! IF ( nd /= deleted_particles ) THEN
2929	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
2930	! CALL local_flush( 9 )
2931	! CALL MPI_ABORT( comm2d, 9999, ierr )
2932	! ENDIF
2933
2934	#else
2935
2936	!
2937	!-- Apply boundary conditions
2938	DO n = 1, number_of_particles
2939
2940	nn = particles(n)%tail_id
2941
2942	IF ( particles(n)%x < -0.5 * dx ) THEN
2943
2944	IF ( ibc_par_lr == 0 ) THEN
2945	!
2946	!-- Cyclic boundary. Relevant coordinate has to be changed.
2947	particles(n)%x = ( nx + 1 ) * dx + particles(n)%x
2948	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2949	i = particles(n)%tailpoints
2950	particle_tail_coordinates(1:i,1,nn) = ( nx + 1 ) * dx + &
2951	particle_tail_coordinates(1:i,1,nn)
2952	ENDIF
2953	ELSEIF ( ibc_par_lr == 1 ) THEN
2954	!
2955	!-- Particle absorption
2956	particle_mask(n) = .FALSE.
2957	deleted_particles = deleted_particles + 1
2958	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2959	tail_mask(nn) = .FALSE.
2960	deleted_tails = deleted_tails + 1
2961	ENDIF
2962	ELSEIF ( ibc_par_lr == 2 ) THEN
2963	!
2964	!-- Particle reflection
2965	particles(n)%x = -dx - particles(n)%x
2966	particles(n)%speed_x = -particles(n)%speed_x
2967	ENDIF
2968
2969	ELSEIF ( particles(n)%x >= ( nx + 0.5 ) * dx ) THEN
2970
2971	IF ( ibc_par_lr == 0 ) THEN
2972	!
2973	!-- Cyclic boundary. Relevant coordinate has to be changed.
2974	particles(n)%x = particles(n)%x - ( nx + 1 ) * dx
2975	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2976	i = particles(n)%tailpoints
2977	particle_tail_coordinates(1:i,1,nn) = - ( nx + 1 ) * dx + &
2978	particle_tail_coordinates(1:i,1,nn)
2979	ENDIF
2980	ELSEIF ( ibc_par_lr == 1 ) THEN
2981	!
2982	!-- Particle absorption
2983	particle_mask(n) = .FALSE.
2984	deleted_particles = deleted_particles + 1
2985	IF ( use_particle_tails .AND. nn /= 0 ) THEN
2986	tail_mask(nn) = .FALSE.
2987	deleted_tails = deleted_tails + 1
2988	ENDIF
2989	ELSEIF ( ibc_par_lr == 2 ) THEN
2990	!
2991	!-- Particle reflection
2992	particles(n)%x = ( nx + 1 ) * dx - particles(n)%x
2993	particles(n)%speed_x = -particles(n)%speed_x
2994	ENDIF
2995
2996	ENDIF
2997
2998	IF ( particles(n)%y < -0.5 * dy ) THEN
2999
3000	IF ( ibc_par_ns == 0 ) THEN
3001	!
3002	!-- Cyclic boundary. Relevant coordinate has to be changed.
3003	particles(n)%y = ( ny + 1 ) * dy + particles(n)%y
3004	IF ( use_particle_tails .AND. nn /= 0 ) THEN
3005	i = particles(n)%tailpoints
3006	particle_tail_coordinates(1:i,2,nn) = ( ny + 1 ) * dy + &
3007	particle_tail_coordinates(1:i,2,nn)
3008	ENDIF
3009	ELSEIF ( ibc_par_ns == 1 ) THEN
3010	!
3011	!-- Particle absorption
3012	particle_mask(n) = .FALSE.
3013	deleted_particles = deleted_particles + 1
3014	IF ( use_particle_tails .AND. nn /= 0 ) THEN
3015	tail_mask(nn) = .FALSE.
3016	deleted_tails = deleted_tails + 1
3017	ENDIF
3018	ELSEIF ( ibc_par_ns == 2 ) THEN
3019	!
3020	!-- Particle reflection
3021	particles(n)%y = -dy - particles(n)%y
3022	particles(n)%speed_y = -particles(n)%speed_y
3023	ENDIF
3024
3025	ELSEIF ( particles(n)%y >= ( ny + 0.5 ) * dy ) THEN
3026
3027	IF ( ibc_par_ns == 0 ) THEN
3028	!
3029	!-- Cyclic boundary. Relevant coordinate has to be changed.
3030	particles(n)%y = particles(n)%y - ( ny + 1 ) * dy
3031	IF ( use_particle_tails .AND. nn /= 0 ) THEN
3032	i = particles(n)%tailpoints
3033	particle_tail_coordinates(1:i,2,nn) = - ( ny + 1 ) * dy + &
3034	particle_tail_coordinates(1:i,2,nn)
3035	ENDIF
3036	ELSEIF ( ibc_par_ns == 1 ) THEN
3037	!
3038	!-- Particle absorption
3039	particle_mask(n) = .FALSE.
3040	deleted_particles = deleted_particles + 1
3041	IF ( use_particle_tails .AND. nn /= 0 ) THEN
3042	tail_mask(nn) = .FALSE.
3043	deleted_tails = deleted_tails + 1
3044	ENDIF
3045	ELSEIF ( ibc_par_ns == 2 ) THEN
3046	!
3047	!-- Particle reflection
3048	particles(n)%y = ( ny + 1 ) * dy - particles(n)%y
3049	particles(n)%speed_y = -particles(n)%speed_y
3050	ENDIF
3051
3052	ENDIF
3053	ENDDO
3054
3055	#endif
3056
3057	!
3058	!-- Apply boundary conditions to those particles that have crossed the top or
3059	!-- bottom boundary and delete those particles, which are older than allowed
3060	DO n = 1, number_of_particles
3061
3062	nn = particles(n)%tail_id
3063
3064	!
3065	!-- Stop if particles have moved further than the length of one
3066	!-- PE subdomain
3067	IF ( ABS(particles(n)%speed_x) > &
3068	((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m) .OR. &
3069	ABS(particles(n)%speed_y) > &
3070	((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) ) THEN
3071
3072	PRINT*, '+++ advec_particles: particle too fast. n = ', n
3073	#if defined( __parallel )
3074	CALL MPI_ABORT( comm2d, 9999, ierr )
3075	#else
3076	CALL local_stop
3077	#endif
3078	ENDIF
3079
3080	IF ( particles(n)%age > particle_maximum_age .AND. &
3081	particle_mask(n) ) &
3082	THEN
3083	particle_mask(n) = .FALSE.
3084	deleted_particles = deleted_particles + 1
3085	IF ( use_particle_tails .AND. nn /= 0 ) THEN
3086	tail_mask(nn) = .FALSE.
3087	deleted_tails = deleted_tails + 1
3088	ENDIF
3089	ENDIF
3090
3091	IF ( particles(n)%z >= zu(nz) .AND. particle_mask(n) ) THEN
3092	IF ( ibc_par_t == 1 ) THEN
3093	!
3094	!-- Particle absorption
3095	particle_mask(n) = .FALSE.
3096	deleted_particles = deleted_particles + 1
3097	IF ( use_particle_tails .AND. nn /= 0 ) THEN
3098	tail_mask(nn) = .FALSE.
3099	deleted_tails = deleted_tails + 1
3100	ENDIF
3101	ELSEIF ( ibc_par_t == 2 ) THEN
3102	!
3103	!-- Particle reflection
3104	particles(n)%z = 2.0 * zu(nz) - particles(n)%z
3105	particles(n)%speed_z = -particles(n)%speed_z
3106	IF ( use_sgs_for_particles .AND. &
3107	particles(n)%speed_z_sgs > 0.0 ) THEN
3108	particles(n)%speed_z_sgs = -particles(n)%speed_z_sgs
3109	ENDIF
3110	IF ( use_particle_tails .AND. nn /= 0 ) THEN
3111	particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
3112	particle_tail_coordinates(1,3,nn)
3113	ENDIF
3114	ENDIF
3115	ENDIF
3116	IF ( particles(n)%z < 0.0 .AND. particle_mask(n) ) THEN
3117	IF ( ibc_par_b == 1 ) THEN
3118	!
3119	!-- Particle absorption
3120	particle_mask(n) = .FALSE.
3121	deleted_particles = deleted_particles + 1
3122	IF ( use_particle_tails .AND. nn /= 0 ) THEN
3123	tail_mask(nn) = .FALSE.
3124	deleted_tails = deleted_tails + 1
3125	ENDIF
3126	ELSEIF ( ibc_par_b == 2 ) THEN
3127	!
3128	!-- Particle reflection
3129	particles(n)%z = -particles(n)%z
3130	particles(n)%speed_z = -particles(n)%speed_z
3131	IF ( use_sgs_for_particles .AND. &
3132	particles(n)%speed_z_sgs < 0.0 ) THEN
3133	particles(n)%speed_z_sgs = -particles(n)%speed_z_sgs
3134	ENDIF
3135	IF ( use_particle_tails .AND. nn /= 0 ) THEN
3136	particle_tail_coordinates(1,3,nn) = 2.0 * zu(nz) - &
3137	particle_tail_coordinates(1,3,nn)
3138	ENDIF
3139	IF ( use_particle_tails .AND. nn /= 0 ) THEN
3140	particle_tail_coordinates(1,3,nn) = &
3141	-particle_tail_coordinates(1,3,nn)
3142	ENDIF
3143	ENDIF
3144	ENDIF
3145	ENDDO
3146
3147	! WRITE ( 9, * ) '*** advec_particles: ##7'
3148	! CALL local_flush( 9 )
3149	! nd = 0
3150	! DO n = 1, number_of_particles
3151	! IF ( .NOT. particle_mask(n) ) nd = nd + 1
3152	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
3153	! THEN
3154	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
3155	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
3156	! CALL local_flush( 9 )
3157	! CALL MPI_ABORT( comm2d, 9999, ierr )
3158	! ENDIF
3159	! ENDDO
3160	! IF ( nd /= deleted_particles ) THEN
3161	! WRITE (9,) '** nd=',nd,' deleted_particles=',deleted_particles
3162	! CALL local_flush( 9 )
3163	! CALL MPI_ABORT( comm2d, 9999, ierr )
3164	! ENDIF
3165
3166	!
3167	!-- Pack particles (eliminate those marked for deletion),
3168	!-- determine new number of particles
3169	IF ( number_of_particles > 0 .AND. deleted_particles > 0 ) THEN
3170	nn = 0
3171	nd = 0
3172	DO n = 1, number_of_particles
3173	IF ( particle_mask(n) ) THEN
3174	nn = nn + 1
3175	particles(nn) = particles(n)
3176	ELSE
3177	nd = nd + 1
3178	ENDIF
3179	ENDDO
3180	! IF ( nd /= deleted_particles ) THEN
3181	! WRITE (9,) '** advec_part nd=',nd,' deleted_particles=',deleted_particles
3182	! CALL local_flush( 9 )
3183	! CALL MPI_ABORT( comm2d, 9999, ierr )
3184	! ENDIF
3185
3186	number_of_particles = number_of_particles - deleted_particles
3187	!
3188	!-- Pack the tails, store the new tail ids and re-assign it to the
3189	!-- respective
3190	!-- particles
3191	IF ( use_particle_tails ) THEN
3192	nn = 0
3193	nd = 0
3194	DO n = 1, number_of_tails
3195	IF ( tail_mask(n) ) THEN
3196	nn = nn + 1
3197	particle_tail_coordinates(:,:,nn) = &
3198	particle_tail_coordinates(:,:,n)
3199	new_tail_id(n) = nn
3200	ELSE
3201	nd = nd + 1
3202	! WRITE (9,*) '+++ n=',n,' (of ',number_of_tails,' #oftails)'
3203	! WRITE (9,*) ' id=',new_tail_id(n)
3204	! CALL local_flush( 9 )
3205	ENDIF
3206	ENDDO
3207	ENDIF
3208
3209	! IF ( nd /= deleted_tails .AND. use_particle_tails ) THEN
3210	! WRITE (9,) '** advec_part nd=',nd,' deleted_tails=',deleted_tails
3211	! CALL local_flush( 9 )
3212	! CALL MPI_ABORT( comm2d, 9999, ierr )
3213	! ENDIF
3214
3215	number_of_tails = number_of_tails - deleted_tails
3216
3217	! nn = 0
3218	DO n = 1, number_of_particles
3219	IF ( particles(n)%tail_id /= 0 ) THEN
3220	! nn = nn + 1
3221	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id > number_of_tails ) THEN
3222	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
3223	! WRITE (9,*) ' tail_id=',particles(n)%tail_id
3224	! WRITE (9,*) ' new_tail_id=', new_tail_id(particles(n)%tail_id), &
3225	! ' of (',number_of_tails,')'
3226	! CALL local_flush( 9 )
3227	! ENDIF
3228	particles(n)%tail_id = new_tail_id(particles(n)%tail_id)
3229	ENDIF
3230	ENDDO
3231
3232	! IF ( nn /= number_of_tails .AND. use_particle_tails ) THEN
3233	! WRITE (9,) '** advec_part #of_tails=',number_of_tails,' nn=',nn
3234	! CALL local_flush( 9 )
3235	! DO n = 1, number_of_particles
3236	! WRITE (9,*) 'prt# ',n,' tail_id=',particles(n)%tail_id, &
3237	! ' x=',particles(n)%x, ' y=',particles(n)%y, &
3238	! ' z=',particles(n)%z
3239	! ENDDO
3240	! CALL MPI_ABORT( comm2d, 9999, ierr )
3241	! ENDIF
3242
3243	ENDIF
3244
3245	! IF ( number_of_particles /= number_of_tails ) THEN
3246	! WRITE (9,*) '--- advec_particles: #7'
3247	! WRITE (9,*) ' #of p=',number_of_particles,' #of t=',number_of_tails
3248	! CALL local_flush( 9 )
3249	! ENDIF
3250	! WRITE ( 9, * ) '*** advec_particles: ##8'
3251	! CALL local_flush( 9 )
3252	! DO n = 1, number_of_particles
3253	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
3254	! THEN
3255	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
3256	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
3257	! CALL local_flush( 9 )
3258	! CALL MPI_ABORT( comm2d, 9999, ierr )
3259	! ENDIF
3260	! ENDDO
3261
3262	! WRITE ( 9, * ) '*** advec_particles: ##9'
3263	! CALL local_flush( 9 )
3264	! DO n = 1, number_of_particles
3265	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
3266	! THEN
3267	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
3268	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
3269	! CALL local_flush( 9 )
3270	! CALL MPI_ABORT( comm2d, 9999, ierr )
3271	! ENDIF
3272	! ENDDO
3273
3274	!
3275	!-- Accumulate the number of particles transferred between the subdomains
3276	#if defined( __parallel )
3277	trlp_count_sum = trlp_count_sum + trlp_count
3278	trlp_count_recv_sum = trlp_count_recv_sum + trlp_count_recv
3279	trrp_count_sum = trrp_count_sum + trrp_count
3280	trrp_count_recv_sum = trrp_count_recv_sum + trrp_count_recv
3281	trsp_count_sum = trsp_count_sum + trsp_count
3282	trsp_count_recv_sum = trsp_count_recv_sum + trsp_count_recv
3283	trnp_count_sum = trnp_count_sum + trnp_count
3284	trnp_count_recv_sum = trnp_count_recv_sum + trnp_count_recv
3285	#endif
3286
3287	IF ( dt_3d_reached ) EXIT
3288
3289	ENDDO ! timestep loop
3290
3291
3292	!
3293	!-- Sort particles in the sequence the gridboxes are stored in the memory
3294	time_sort_particles = time_sort_particles + dt_3d
3295	IF ( time_sort_particles >= dt_sort_particles ) THEN
3296	CALL sort_particles
3297	time_sort_particles = MOD( time_sort_particles, &
3298	MAX( dt_sort_particles, dt_3d ) )
3299	ENDIF
3300
3301
3302	!
3303	!-- Re-evaluate the weighting factors. After advection, particles within a
3304	!-- grid box may have different weighting factors if some have been advected
3305	!-- from a neighbouring box. The factors are re-evaluated so that they are
3306	!-- the same for all particles of one box. This procedure must conserve the
3307	!-- liquid water content within one box.
3308	IF ( cloud_droplets ) THEN
3309
3310	CALL cpu_log( log_point_s(45), 'advec_part_reeval_we', 'start' )
3311
3312	ql = 0.0; ql_v = 0.0; ql_vp = 0.0
3313
3314	!
3315	!-- Re-calculate the weighting factors and calculate the liquid water content
3316	DO i = nxl, nxr
3317	DO j = nys, nyn
3318	DO k = nzb, nzt+1
3319
3320	!
3321	!-- Calculate the total volume of particles in the boxes (ql_vp) as
3322	!-- well as the real volume (ql_v, weighting factor has to be
3323	!-- included)
3324	psi = prt_start_index(k,j,i)
3325	DO n = psi, psi+prt_count(k,j,i)-1
3326	ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%radius**3
3327
3328	ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * &
3329	particles(n)%radius**3
3330	ENDDO
3331
3332	!
3333	!-- Re-calculate the weighting factors and calculate the liquid
3334	!-- water content
3335	IF ( ql_vp(k,j,i) /= 0.0 ) THEN
3336	ql_vp(k,j,i) = ql_v(k,j,i) / ql_vp(k,j,i)
3337	ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi * &
3338	ql_v(k,j,i) / &
3339	( rho_surface * dx * dy * dz )
3340	ELSE
3341	ql(k,j,i) = 0.0
3342	ENDIF
3343
3344	!
3345	!-- Re-assign the weighting factor to the particles
3346	DO n = psi, psi+prt_count(k,j,i)-1
3347	particles(n)%weight_factor = ql_vp(k,j,i)
3348	ENDDO
3349
3350	ENDDO
3351	ENDDO
3352	ENDDO
3353
3354	CALL cpu_log( log_point_s(45), 'advec_part_reeval_we', 'stop' )
3355
3356	ENDIF
3357
3358	!
3359	!-- Set particle attributes defined by the user
3360	CALL user_particle_attributes
3361	! WRITE ( 9, * ) '*** advec_particles: ##10'
3362	! CALL local_flush( 9 )
3363	! DO n = 1, number_of_particles
3364	! IF ( particles(n)%tail_id<0 .OR. particles(n)%tail_id>number_of_tails ) &
3365	! THEN
3366	! WRITE (9,*) '+++ n=',n,' (of ',number_of_particles,')'
3367	! WRITE (9,*) ' id=',particles(n)%tail_id,' of (',number_of_tails,')'
3368	! CALL local_flush( 9 )
3369	! CALL MPI_ABORT( comm2d, 9999, ierr )
3370	! ENDIF
3371	! ENDDO
3372
3373	!
3374	!-- If necessary, add the actual particle positions to the particle tails
3375	IF ( use_particle_tails ) THEN
3376
3377	distance = 0.0
3378	DO n = 1, number_of_particles
3379
3380	nn = particles(n)%tail_id
3381
3382	IF ( nn /= 0 ) THEN
3383	!
3384	!-- Calculate the distance between the actual particle position and the
3385	!-- next tailpoint
3386	! WRITE ( 9, * ) '*** advec_particles: ##10.1 nn=',nn
3387	! CALL local_flush( 9 )
3388	IF ( minimum_tailpoint_distance /= 0.0 ) THEN
3389	distance = ( particle_tail_coordinates(1,1,nn) - &
3390	particle_tail_coordinates(2,1,nn) )**2 + &
3391	( particle_tail_coordinates(1,2,nn) - &
3392	particle_tail_coordinates(2,2,nn) )**2 + &
3393	( particle_tail_coordinates(1,3,nn) - &
3394	particle_tail_coordinates(2,3,nn) )**2
3395	ENDIF
3396	! WRITE ( 9, * ) '*** advec_particles: ##10.2'
3397	! CALL local_flush( 9 )
3398	!
3399	!-- First, increase the index of all existings tailpoints by one
3400	IF ( distance >= minimum_tailpoint_distance ) THEN
3401	DO i = particles(n)%tailpoints, 1, -1
3402	particle_tail_coordinates(i+1,:,nn) = &
3403	particle_tail_coordinates(i,:,nn)
3404	ENDDO
3405	!
3406	!-- Increase the counter which contains the number of tailpoints.
3407	!-- This must always be smaller than the given maximum number of
3408	!-- tailpoints because otherwise the index bounds of
3409	!-- particle_tail_coordinates would be exceeded
3410	IF ( particles(n)%tailpoints < maximum_number_of_tailpoints-1 )&
3411	THEN
3412	particles(n)%tailpoints = particles(n)%tailpoints + 1
3413	ENDIF
3414	ENDIF
3415	! WRITE ( 9, * ) '*** advec_particles: ##10.3'
3416	! CALL local_flush( 9 )
3417	!
3418	!-- In any case, store the new point at the beginning of the tail
3419	particle_tail_coordinates(1,1,nn) = particles(n)%x
3420	particle_tail_coordinates(1,2,nn) = particles(n)%y
3421	particle_tail_coordinates(1,3,nn) = particles(n)%z
3422	particle_tail_coordinates(1,4,nn) = particles(n)%color
3423	! WRITE ( 9, * ) '*** advec_particles: ##10.4'
3424	! CALL local_flush( 9 )
3425	!
3426	!-- Increase the age of the tailpoints
3427	IF ( minimum_tailpoint_distance /= 0.0 ) THEN
3428	particle_tail_coordinates(2:particles(n)%tailpoints,5,nn) = &
3429	particle_tail_coordinates(2:particles(n)%tailpoints,5,nn) + &
3430	dt_3d
3431	!
3432	!-- Delete the last tailpoint, if it has exceeded its maximum age
3433	IF ( particle_tail_coordinates(particles(n)%tailpoints,5,nn) > &
3434	maximum_tailpoint_age ) THEN
3435	particles(n)%tailpoints = particles(n)%tailpoints - 1
3436	ENDIF
3437	ENDIF
3438	! WRITE ( 9, * ) '*** advec_particles: ##10.5'
3439	! CALL local_flush( 9 )
3440
3441	ENDIF
3442
3443	ENDDO
3444
3445	ENDIF
3446	! WRITE ( 9, * ) '*** advec_particles: ##11'
3447	! CALL local_flush( 9 )
3448
3449	!
3450	!-- Write particle statistics on file
3451	IF ( write_particle_statistics ) THEN
3452	CALL check_open( 80 )
3453	#if defined( __parallel )
3454	WRITE ( 80, 8000 ) current_timestep_number+1, simulated_time+dt_3d, &
3455	number_of_particles, pleft, trlp_count_sum, &
3456	trlp_count_recv_sum, pright, trrp_count_sum, &
3457	trrp_count_recv_sum, psouth, trsp_count_sum, &
3458	trsp_count_recv_sum, pnorth, trnp_count_sum, &
3459	trnp_count_recv_sum, maximum_number_of_particles
3460	CALL close_file( 80 )
3461	#else
3462	WRITE ( 80, 8000 ) current_timestep_number+1, simulated_time+dt_3d, &
3463	number_of_particles, maximum_number_of_particles
3464	#endif
3465	ENDIF
3466
3467	CALL cpu_log( log_point(25), 'advec_particles', 'stop' )
3468
3469	!
3470	!-- Formats
3471	8000 FORMAT (I6,1X,F7.2,4X,I6,5X,4(I3,1X,I4,'/',I4,2X),6X,I6)
3472
3473	END SUBROUTINE advec_particles
3474
3475
3476	SUBROUTINE allocate_prt_memory( number_of_new_particles )
3477
3478	!------------------------------------------------------------------------------!
3479	! Description:
3480	! ------------
3481	! Extend particle memory
3482	!------------------------------------------------------------------------------!
3483
3484	USE particle_attributes
3485
3486	IMPLICIT NONE
3487
3488	INTEGER :: new_maximum_number, number_of_new_particles
3489
3490	LOGICAL, DIMENSION(:), ALLOCATABLE :: tmp_particle_mask
3491
3492	TYPE(particle_type), DIMENSION(:), ALLOCATABLE :: tmp_particles
3493
3494
3495	new_maximum_number = maximum_number_of_particles + &
3496	MAX( 5*number_of_new_particles, number_of_initial_particles )
3497
3498	IF ( write_particle_statistics ) THEN
3499	CALL check_open( 80 )
3500	WRITE ( 80, '(''*** Request: '', I7, '' new_maximum_number(prt)'')' ) &
3501	new_maximum_number
3502	CALL close_file( 80 )
3503	ENDIF
3504
3505	ALLOCATE( tmp_particles(maximum_number_of_particles), &
3506	tmp_particle_mask(maximum_number_of_particles) )
3507	tmp_particles = particles
3508	tmp_particle_mask = particle_mask
3509
3510	DEALLOCATE( particles, particle_mask )
3511	ALLOCATE( particles(new_maximum_number), &
3512	particle_mask(new_maximum_number) )
3513	maximum_number_of_particles = new_maximum_number
3514
3515	particles(1:number_of_particles) = tmp_particles(1:number_of_particles)
3516	particle_mask(1:number_of_particles) = &
3517	tmp_particle_mask(1:number_of_particles)
3518	particle_mask(number_of_particles+1:maximum_number_of_particles) = .TRUE.
3519	DEALLOCATE( tmp_particles, tmp_particle_mask )
3520
3521	END SUBROUTINE allocate_prt_memory
3522
3523
3524	SUBROUTINE allocate_tail_memory( number_of_new_tails )
3525
3526	!------------------------------------------------------------------------------!
3527	! Description:
3528	! ------------
3529	! Extend tail memory
3530	!------------------------------------------------------------------------------!
3531
3532	USE particle_attributes
3533
3534	IMPLICIT NONE
3535
3536	INTEGER :: new_maximum_number, number_of_new_tails
3537
3538	LOGICAL, DIMENSION(maximum_number_of_tails) :: tmp_tail_mask
3539
3540	REAL, DIMENSION(maximum_number_of_tailpoints,5,maximum_number_of_tails) :: &
3541	tmp_tail
3542
3543
3544	new_maximum_number = maximum_number_of_tails + &
3545	MAX( 5*number_of_new_tails, number_of_initial_tails )
3546
3547	IF ( write_particle_statistics ) THEN
3548	CALL check_open( 80 )
3549	WRITE ( 80, '(''*** Request: '', I5, '' new_maximum_number(tails)'')' ) &
3550	new_maximum_number
3551	CALL close_file( 80 )
3552	ENDIF
3553	WRITE (9,) '** Request: ',new_maximum_number,' new_maximum_number(tails)'
3554	! CALL local_flush( 9 )
3555
3556	tmp_tail(:,:,1:number_of_tails) = &
3557	particle_tail_coordinates(:,:,1:number_of_tails)
3558	tmp_tail_mask(1:number_of_tails) = tail_mask(1:number_of_tails)
3559
3560	DEALLOCATE( new_tail_id, particle_tail_coordinates, tail_mask )
3561	ALLOCATE( new_tail_id(new_maximum_number), &
3562	particle_tail_coordinates(maximum_number_of_tailpoints,5, &
3563	new_maximum_number), &
3564	tail_mask(new_maximum_number) )
3565	maximum_number_of_tails = new_maximum_number
3566
3567	particle_tail_coordinates = 0.0
3568	particle_tail_coordinates(:,:,1:number_of_tails) = &
3569	tmp_tail(:,:,1:number_of_tails)
3570	tail_mask(1:number_of_tails) = tmp_tail_mask(1:number_of_tails)
3571	tail_mask(number_of_tails+1:maximum_number_of_tails) = .TRUE.
3572
3573	END SUBROUTINE allocate_tail_memory
3574
3575
3576	SUBROUTINE output_particles_netcdf
3577	#if defined( __netcdf )
3578
3579	USE control_parameters
3580	USE netcdf_control
3581	USE particle_attributes
3582
3583	IMPLICIT NONE
3584
3585
3586	CALL check_open( 108 )
3587
3588	!
3589	!-- Update the NetCDF time axis
3590	prt_time_count = prt_time_count + 1
3591
3592	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_time_prt, (/ simulated_time /), &
3593	start = (/ prt_time_count /), count = (/ 1 /) )
3594	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 1 )
3595
3596	!
3597	!-- Output the real number of particles used
3598	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_rnop_prt, &
3599	(/ number_of_particles /), &
3600	start = (/ prt_time_count /), count = (/ 1 /) )
3601	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 2 )
3602
3603	!
3604	!-- Output all particle attributes
3605	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(1), particles%age, &
3606	start = (/ 1, prt_time_count /), &
3607	count = (/ maximum_number_of_particles /) )
3608	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 3 )
3609
3610	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(2), particles%dvrp_psize, &
3611	start = (/ 1, prt_time_count /), &
3612	count = (/ maximum_number_of_particles /) )
3613	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 4 )
3614
3615	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(3), particles%origin_x, &
3616	start = (/ 1, prt_time_count /), &
3617	count = (/ maximum_number_of_particles /) )
3618	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 5 )
3619
3620	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(4), particles%origin_y, &
3621	start = (/ 1, prt_time_count /), &
3622	count = (/ maximum_number_of_particles /) )
3623	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 6 )
3624
3625	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(5), particles%origin_z, &
3626	start = (/ 1, prt_time_count /), &
3627	count = (/ maximum_number_of_particles /) )
3628	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 7 )
3629
3630	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(6), particles%radius, &
3631	start = (/ 1, prt_time_count /), &
3632	count = (/ maximum_number_of_particles /) )
3633	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 8 )
3634
3635	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(7), particles%speed_x, &
3636	start = (/ 1, prt_time_count /), &
3637	count = (/ maximum_number_of_particles /) )
3638	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 9 )
3639
3640	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(8), particles%speed_y, &
3641	start = (/ 1, prt_time_count /), &
3642	count = (/ maximum_number_of_particles /) )
3643	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 10 )
3644
3645	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(9), particles%speed_z, &
3646	start = (/ 1, prt_time_count /), &
3647	count = (/ maximum_number_of_particles /) )
3648	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 11 )
3649
3650	nc_stat = NF90_PUT_VAR( id_set_prt,id_var_prt(10),particles%weight_factor,&
3651	start = (/ 1, prt_time_count /), &
3652	count = (/ maximum_number_of_particles /) )
3653	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 12 )
3654
3655	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(11), particles%x, &
3656	start = (/ 1, prt_time_count /), &
3657	count = (/ maximum_number_of_particles /) )
3658	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 13 )
3659
3660	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(12), particles%y, &
3661	start = (/ 1, prt_time_count /), &
3662	count = (/ maximum_number_of_particles /) )
3663	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 14 )
3664
3665	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(13), particles%z, &
3666	start = (/ 1, prt_time_count /), &
3667	count = (/ maximum_number_of_particles /) )
3668	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 15 )
3669
3670	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(14), particles%color, &
3671	start = (/ 1, prt_time_count /), &
3672	count = (/ maximum_number_of_particles /) )
3673	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 16 )
3674
3675	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(15), particles%group, &
3676	start = (/ 1, prt_time_count /), &
3677	count = (/ maximum_number_of_particles /) )
3678	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 17 )
3679
3680	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(16), particles%tailpoints, &
3681	start = (/ 1, prt_time_count /), &
3682	count = (/ maximum_number_of_particles /) )
3683	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 18 )
3684
3685	nc_stat = NF90_PUT_VAR( id_set_prt, id_var_prt(17), particles%tail_id, &
3686	start = (/ 1, prt_time_count /), &
3687	count = (/ maximum_number_of_particles /) )
3688	IF (nc_stat /= NF90_NOERR) CALL handle_netcdf_error( 19 )
3689
3690	#endif
3691	END SUBROUTINE output_particles_netcdf
3692
3693
3694	SUBROUTINE write_particles
3695
3696	!------------------------------------------------------------------------------!
3697	! Description:
3698	! ------------
3699	! Write particle data on restart file
3700	!------------------------------------------------------------------------------!
3701
3702	USE control_parameters
3703	USE particle_attributes
3704	USE pegrid
3705
3706	IMPLICIT NONE
3707
3708	CHARACTER (LEN=10) :: particle_binary_version
3709
3710	!
3711	!-- First open the output unit.
3712	IF ( myid_char == '' ) THEN
3713	OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT'//myid_char, &
3714	FORM='UNFORMATTED')
3715	ELSE
3716	IF ( myid == 0 ) CALL local_system( 'mkdir PARTICLE_RESTART_DATA_OUT' )
3717	#if defined( __parallel )
3718	!
3719	!-- Set a barrier in order to allow that thereafter all other processors
3720	!-- in the directory created by PE0 can open their file
3721	CALL MPI_BARRIER( comm2d, ierr )
3722	#endif
3723	OPEN ( 90, FILE='PARTICLE_RESTART_DATA_OUT/'//myid_char, &
3724	FORM='UNFORMATTED' )
3725	ENDIF
3726
3727	!
3728	!-- Write the version number of the binary format.
3729	!-- Attention: After changes to the following output commands the version
3730	!-- --------- number of the variable particle_binary_version must be changed!
3731	!-- Also, the version number and the list of arrays to be read in
3732	!-- init_particles must be adjusted accordingly.
3733	particle_binary_version = '3.0'
3734	WRITE ( 90 ) particle_binary_version
3735
3736	!
3737	!-- Write some particle parameters, the size of the particle arrays as well as
3738	!-- other dvrp-plot variables.
3739	WRITE ( 90 ) bc_par_b, bc_par_lr, bc_par_ns, bc_par_t, &
3740	maximum_number_of_particles, maximum_number_of_tailpoints, &
3741	maximum_number_of_tails, number_of_initial_particles, &
3742	number_of_particles, number_of_particle_groups, &
3743	number_of_tails, particle_groups, time_prel, &
3744	time_write_particle_data, uniform_particles
3745
3746	IF ( number_of_initial_particles /= 0 ) WRITE ( 90 ) initial_particles
3747
3748	WRITE ( 90 ) prt_count, prt_start_index
3749	WRITE ( 90 ) particles
3750
3751	IF ( use_particle_tails ) THEN
3752	WRITE ( 90 ) particle_tail_coordinates
3753	ENDIF
3754
3755	CLOSE ( 90 )
3756
3757	END SUBROUTINE write_particles
3758
3759
3760	SUBROUTINE collision_efficiency( mean_r, r, e)
3761	!------------------------------------------------------------------------------!
3762	! Description:
3763	! ------------
3764	! Interpolate collision efficiency from table
3765	!------------------------------------------------------------------------------!
3766
3767	IMPLICIT NONE
3768
3769	INTEGER :: i, j, k
3770
3771	LOGICAL, SAVE :: first = .TRUE.
3772
3773	REAL :: aa, bb, cc, dd, dx, dy, e, gg, mean_r, mean_rm, r, rm, &
3774	x, y
3775
3776	REAL, DIMENSION(1:9), SAVE :: collected_r = 0.0
3777	REAL, DIMENSION(1:19), SAVE :: collector_r = 0.0
3778	REAL, DIMENSION(1:9,1:19), SAVE :: ef = 0.0
3779
3780	mean_rm = mean_r * 1.0E06
3781	rm = r * 1.0E06
3782
3783	IF ( first ) THEN
3784	collected_r = (/ 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0, 25.0 /)
3785	collector_r = (/ 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 80.0, 100.0, 150.0,&
3786	200.0, 300.0, 400.0, 500.0, 600.0, 1000.0, 1400.0, &
3787	1800.0, 2400.0, 3000.0 /)
3788	ef(:,1) = (/0.017, 0.027, 0.037, 0.052, 0.052, 0.052, 0.052, 0.0, 0.0 /)
3789	ef(:,2) = (/0.001, 0.016, 0.027, 0.060, 0.12, 0.17, 0.17, 0.17, 0.0 /)
3790	ef(:,3) = (/0.001, 0.001, 0.02, 0.13, 0.28, 0.37, 0.54, 0.55, 0.47/)
3791	ef(:,4) = (/0.001, 0.001, 0.02, 0.23, 0.4, 0.55, 0.7, 0.75, 0.75/)
3792	ef(:,5) = (/0.01, 0.01, 0.03, 0.3, 0.4, 0.58, 0.73, 0.75, 0.79/)
3793	ef(:,6) = (/0.01, 0.01, 0.13, 0.38, 0.57, 0.68, 0.80, 0.86, 0.91/)
3794	ef(:,7) = (/0.01, 0.085, 0.23, 0.52, 0.68, 0.76, 0.86, 0.92, 0.95/)
3795	ef(:,8) = (/0.01, 0.14, 0.32, 0.60, 0.73, 0.81, 0.90, 0.94, 0.96/)
3796	ef(:,9) = (/0.025, 0.25, 0.43, 0.66, 0.78, 0.83, 0.92, 0.95, 0.96/)
3797	ef(:,10)= (/0.039, 0.3, 0.46, 0.69, 0.81, 0.87, 0.93, 0.95, 0.96/)
3798	ef(:,11)= (/0.095, 0.33, 0.51, 0.72, 0.82, 0.87, 0.93, 0.96, 0.97/)
3799	ef(:,12)= (/0.098, 0.36, 0.51, 0.73, 0.83, 0.88, 0.93, 0.96, 0.97/)
3800	ef(:,13)= (/0.1, 0.36, 0.52, 0.74, 0.83, 0.88, 0.93, 0.96, 0.97/)
3801	ef(:,14)= (/0.17, 0.4, 0.54, 0.72, 0.83, 0.88, 0.94, 0.98, 1.0 /)
3802	ef(:,15)= (/0.15, 0.37, 0.52, 0.74, 0.82, 0.88, 0.94, 0.98, 1.0 /)
3803	ef(:,16)= (/0.11, 0.34, 0.49, 0.71, 0.83, 0.88, 0.94, 0.95, 1.0 /)
3804	ef(:,17)= (/0.08, 0.29, 0.45, 0.68, 0.8, 0.86, 0.96, 0.94, 1.0 /)
3805	ef(:,18)= (/0.04, 0.22, 0.39, 0.62, 0.75, 0.83, 0.92, 0.96, 1.0 /)
3806	ef(:,19)= (/0.02, 0.16, 0.33, 0.55, 0.71, 0.81, 0.90, 0.94, 1.0 /)
3807	ENDIF
3808
3809	DO k = 1, 8
3810	IF ( collected_r(k) <= mean_rm ) i = k
3811	ENDDO
3812
3813	DO k = 1, 18
3814	IF ( collector_r(k) <= rm ) j = k
3815	ENDDO
3816
3817	IF ( rm < 10.0 ) THEN
3818	e = 0.0
3819	ELSEIF ( mean_rm < 2.0 ) THEN
3820	e = 0.001
3821	ELSEIF ( mean_rm >= 25.0 ) THEN
3822	IF( j <= 3 ) e = 0.55
3823	IF( j == 4 ) e = 0.8
3824	IF( j == 5 ) e = 0.9
3825	IF( j >=6 ) e = 1.0
3826	ELSEIF ( rm >= 3000.0 ) THEN
3827	e = 1.0
3828	ELSE
3829	x = mean_rm - collected_r(i)
3830	y = rm - collected_r(j)
3831	dx = collected_r(i+1) - collected_r(i)
3832	dy = collector_r(j+1) - collector_r(j)
3833	aa = x2 + y2
3834	bb = ( dx - x )2 + y2
3835	cc = x2 + ( dy - y )2
3836	dd = ( dx - x )2 + ( dy - y )2
3837	gg = aa + bb + cc + dd
3838
3839	e = ( (gg-aa)ef(i,j) + (gg-bb)ef(i+1,j) + (gg-cc)*ef(i,j+1) + &
3840	(gg-dd)ef(i+1,j+1) ) / (3.0gg)
3841	ENDIF
3842
3843	END SUBROUTINE collision_efficiency
3844
3845
3846
3847	SUBROUTINE sort_particles
3848
3849	!------------------------------------------------------------------------------!
3850	! Description:
3851	! ------------
3852	! Sort particles in the sequence the grid boxes are stored in memory
3853	!------------------------------------------------------------------------------!
3854
3855	USE arrays_3d
3856	USE control_parameters
3857	USE cpulog
3858	USE grid_variables
3859	USE indices
3860	USE interfaces
3861	USE particle_attributes
3862
3863	IMPLICIT NONE
3864
3865	INTEGER :: i, ilow, j, k, n
3866
3867	TYPE(particle_type), DIMENSION(1:number_of_particles) :: particles_temp
3868
3869
3870	CALL cpu_log( log_point_s(47), 'sort_particles', 'start' )
3871
3872	!
3873	!-- Initialize the array used for counting and indexing the particles
3874	prt_count = 0
3875
3876	!
3877	!-- Count the particles per gridbox
3878	DO n = 1, number_of_particles
3879
3880	i = ( particles(n)%x + 0.5 * dx ) * ddx
3881	j = ( particles(n)%y + 0.5 * dy ) * ddy
3882	k = particles(n)%z / dz + 1 ! only exact if equidistant
3883
3884	prt_count(k,j,i) = prt_count(k,j,i) + 1
3885
3886	IF ( i < nxl .OR. i > nxr .OR. j < nys .OR. j > nyn .OR. k < nzb+1 .OR. &
3887	k > nzt ) THEN
3888	PRINT*, '+++ sort_particles: particle out of range: i=', i, ' j=', &
3889	j, ' k=', k
3890	PRINT*, ' nxl=', nxl, ' nxr=', nxr, &
3891	' nys=', nys, ' nyn=', nyn, &
3892	' nzb=', nzb, ' nzt=', nzt
3893	ENDIF
3894
3895	ENDDO
3896
3897	!
3898	!-- Calculate the lower indices of those ranges of the particles-array
3899	!-- containing particles which belong to the same gridpox i,j,k
3900	ilow = 1
3901	DO i = nxl, nxr
3902	DO j = nys, nyn
3903	DO k = nzb+1, nzt
3904	prt_start_index(k,j,i) = ilow
3905	ilow = ilow + prt_count(k,j,i)
3906	ENDDO
3907	ENDDO
3908	ENDDO
3909
3910	!
3911	!-- Sorting the particles
3912	DO n = 1, number_of_particles
3913
3914	i = ( particles(n)%x + 0.5 * dx ) * ddx
3915	j = ( particles(n)%y + 0.5 * dy ) * ddy
3916	k = particles(n)%z / dz + 1 ! only exact if equidistant
3917
3918	particles_temp(prt_start_index(k,j,i)) = particles(n)
3919
3920	prt_start_index(k,j,i) = prt_start_index(k,j,i) + 1
3921
3922	ENDDO
3923
3924	particles(1:number_of_particles) = particles_temp
3925
3926	!
3927	!-- Reset the index array to the actual start position
3928	DO i = nxl, nxr
3929	DO j = nys, nyn
3930	DO k = nzb+1, nzt
3931	prt_start_index(k,j,i) = prt_start_index(k,j,i) - prt_count(k,j,i)
3932	ENDDO
3933	ENDDO
3934	ENDDO
3935
3936	CALL cpu_log( log_point_s(47), 'sort_particles', 'stop' )
3937
3938	END SUBROUTINE sort_particles

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |