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

Last change on this file since 519 was 519, checked in by raasch, 14 years ago

NetCDF4 support for particle data; special character allowed for NetCDF variable names

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