Home

Context Navigation

advec_particles.f90 @ 622

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

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:

Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:

Property svn:keywords set to Id

File size: 171.2 KB

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

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

Context Navigation

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

Download in other formats: