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

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

small bugfixes in SGS part of adved_particles

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