Home

Context Navigation

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

Last change on this file since 114 was 114, checked in by raasch, 17 years ago
preliminary updates for implementing buildings in poismg
Property svn:keywords set to `Id`
File size: 163.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |