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

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

comments prepared for 3.1c

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