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

Last change on this file since 392 was 392, checked in by raasch, 12 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 raasch $
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 = x**2          + y**2
233          bb = ( dx - x )**2 + y**2
234          cc = x**2          + ( 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_r**3 - particles(n)%radius**3 ) / &
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 = x**2          + y**2
444                      bb = ( dx - x )**2 + y**2
445                      cc = x**2          + ( 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 = x**2          + y**2
475                      bb = ( dx - x )**2 + y**2
476                      cc = x**2          + ( 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 = x**2          + y**2
506                      bb = ( dx - x )**2 + y**2
507                      cc = x**2          + ( 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 = x**2          + y**2
537                      bb = ( dx - x )**2 + y**2
538                      cc = x**2          + ( 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 = x**2          + y**2
855          bb = ( dx - x )**2 + y**2
856          cc = x**2          + ( 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 = x**2          + y**2
882          bb = ( dx - x )**2 + y**2
883          cc = x**2          + ( 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 = x**2          + y**2
910             bb = ( dx - x )**2 + y**2
911             cc = x**2          + ( 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                = particles(n)%x - i * dx
947                y  = particles(n)%y - j * dy
948                aa = x**2          + y**2
949                bb = ( dx - x )**2 + y**2
950                cc = x**2          + ( 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                   = particles(n)%x - i * dx
1189                   y  = particles(n)%y - j * dy
1190                   aa = x**2          + y**2
1191                   bb = ( dx - x )**2 + y**2
1192                   cc = x**2          + ( 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.0*sigma, 5.0*sigma] 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)%<=  ( nxl - 0.5 ) * dx )  THEN
1980                         particles(n)%x = ( nxl - 0.4999999999 ) * dx
1981                      ELSEIF ( particles(n)%>=  ( 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)%<=  ( nys - 0.5 ) * dy )  THEN
1991                         particles(n)%y = ( nys - 0.4999999999 ) * dy
1992                      ELSEIF ( particles(n)%>=  ( 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                         = 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
34848000 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       = 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 = x**2 + y**2
3847       bb = ( dx - x )**2 + y**2
3848       cc = x**2 + ( 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.0*gg)
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.