Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |