source: palm/tags/release-3.2a/SOURCE/advec_particles.f90 @ 1887

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

vorlaeufige Standalone-Version fuer Linux-Cluster

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