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

Last change on this file since 73 was 64, checked in by raasch, 17 years ago

sendrecv adjustments in advec_particles

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