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

Last change on this file since 225 was 212, checked in by raasch, 15 years ago

manual updated for changes in the user interface

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