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

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

preliminary updates for implementing buildings in poismg

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