Home

Context Navigation

advec_particles.f90 @ 759

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

New:
---

The number of parallel I/O operations can be limited with new mrun-option -w.
(advec_particles, data_output_2d, data_output_3d, header, init_grid, init_pegrid, init_3d_model, modules, palm, parin, write_3d_binary)

Changed:

mrun option -T is obligatory

Errors:

Bugfix: No zero assignments to volume_flow_initial and volume_flow_area in
case of normal restart runs. (init_3d_model)

initialization of u_0, v_0. This is just to avoid access of uninitialized
memory in exchange_horiz_2d, which causes respective error messages
when the Intel thread checker (inspector) is used. (production_e)

Bugfix for ts limitation (prandtl_fluxes)

Property svn:keywords set to Id

File size: 171.7 KB

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

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

Context Navigation

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

Download in other formats: