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

Last change on this file since 226 was 226, checked in by raasch, 15 years ago

preparations for the next release

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