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

Last change on this file since 75 was 75, checked in by raasch, 14 years ago

preliminary update for changes concerning non-cyclic boundary conditions

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