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

Last change on this file since 208 was 198, checked in by raasch, 16 years ago

file headers updated for the next release 3.5

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