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

Last change on this file since 230 was 230, checked in by heinze, 14 years ago

Output of messages replaced by message handling routine

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