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

Last change on this file since 150 was 150, checked in by raasch, 16 years ago

particle advection allowed for ocean runs

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