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

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

further updates for dvr output, bugfix in advec_particles concerning particle boundary condition

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