Home

Context Navigation

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

Last change on this file since 60 was 60, checked in by raasch, 17 years ago
preliminary update of further changes, running
Property svn:keywords set to `Id`
File size: 162.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |