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

Last change on this file since 292 was 276, checked in by raasch, 15 years ago

bugfix and change concerning particle advection

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