Home

Context Navigation

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

Last change on this file since 76 was 75, checked in by raasch, 17 years ago
preliminary update for changes concerning non-cyclic boundary conditions
Property svn:keywords set to `Id`
File size: 163.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |