Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |