source: palm/tags/release-3.1b/SOURCE/advec_particles.f90

Last change on this file was 4, checked in by raasch, 17 years ago

Id keyword set as property for all *.f90 files

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