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

Last change on this file since 420 was 420, checked in by franke, 13 years ago

collision of cloud droplets has changed in advec_particles
bugfixes in advec_particles and collision_efficiency
change in disturb_field to make runs reproducible on HLRN

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