Home

Context Navigation

advec_particles.f90 @ 825

Last change on this file since 825 was 825, checked in by raasch, 13 years ago

New:
---

droplet growth by condensation may include curvature and solution effects.
Steered by new inipar-parameter curvature_solution_effects.
(advec_particles, check_parameters, header, init_cloud_physics, init_particles, modules, parin, read_var_list, write_var_list)

mean/min/max particle radius added as output quantity. (data_output_ptseries, modules)

Changed:

Initialisation of temporary particle array for resorting removed.
(advec_particles)

particle attributes speed_x|y|z_sgs renamed rvar1|2|3.
(advec_particles, data_output_ptseries, modules, init_particles, particle_boundary_conds)

routine wang_kernel and respective module renamed lpm_collision_kernels.
Package (particle) parameters wang_collision_kernel and turbulence_effects_on_collision
replaced by parameter collision_kernel.
(Makefile, advec_particles, check_parameters, diffusion_e, init_3d_model, modules, package_parin, time_integration, new: lpm_collision_kernels, deleted: wang_kernel)

Errors:

Property svn:keywords set to Id

File size: 184.6 KB

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

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

Context Navigation

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

Download in other formats: