Home

Context Navigation

advec_particles.f90 @ 413

Last change on this file since 413 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

Property svn:keywords set to Id

File size: 165.8 KB

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

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

Context Navigation

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

Download in other formats: