Home

Context Navigation

advec_particles.f90 @ 63

Last change on this file since 63 was 63, checked in by raasch, 18 years ago
preliminary changes concerning update of BC-scheme
Property svn:keywords set to `Id`
File size: 163.2 KB

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

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

Download in other formats:

�