Home

Context Navigation

advec_particles.f90 @ 826

Last change on this file since 826 was 826, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 184.7 KB

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