Home

Context Navigation

advec_particles.f90 @ 145

Last change on this file since 145 was 139, checked in by raasch, 17 years ago

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:

Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:

Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

Property svn:keywords set to Id

File size: 164.5 KB

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

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

Context Navigation

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

Download in other formats: