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

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

last commit documented

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