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

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

bugfixes in advec_particles and init_particles

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