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

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

updating comments and rc-file

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