Home

Context Navigation

advec_particles.f90 @ 229

Last change on this file since 229 was 229, checked in by raasch, 16 years ago
bugfixes concerning particle tails
Property svn:keywords set to `Id`
File size: 165.4 KB

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