Home

Context Navigation

advec_particles.f90 @ 828

Last change on this file since 828 was 828, checked in by raasch, 12 years ago

New:
---

Changed:

Optimization of collision kernels. Collision tables can be calculated once at
simulation start for defined radius (and dissipation) classes instead of
re-calculating them at every timestep and for the particle ensemble in
every gridbox.
For this purpose the particle feature color is renamed class.
New parpar parameters radius_classes and dissipation_classes.
(Makefile, advec_particles, check_parameters, data_output_dvrp, header, init_particles, lpm_collision_kernels, modules, package_parin, set_particle_attributes)

Lower limit for droplet radius changed from 1E-7 to 1E-8.
(advec_particles)

Complete re-formatting of collision code (including changes in names of
variables, modules and subroutines).
(advec_particles, lpm_collision_kernels)

Errors:

Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
(advec_particles, lpm_collision_kernels)

Property svn:keywords set to Id

File size: 190.5 KB

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

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

Context Navigation

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

Download in other formats: