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

Last change on this file since 792 was 792, checked in by raasch, 11 years ago

further adjustments for speedup of particle code

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