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

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

preliminary update of further changes, not running

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