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

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

preliminary update of further changes, advec_particles is not running!

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