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

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

New:
---

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

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

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

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

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

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

Changed:


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

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

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

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

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

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

non-zero wall_heatflux,

  • it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

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

Errors:


Bugfix: summation of sums_l_l in diffusivities.

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

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

Bugfix: model_string needed a default value (combine_plot_fields)

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

Bugfix: assignment of fluxes at walls

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

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

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