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

Last change on this file since 116 was 116, checked in by raasch, 17 years ago

further preliminary updates concerning particle sorting and documentation

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