Home

Context Navigation

advec_particles.f90 @ 845

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

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

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

Context Navigation

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

Download in other formats: