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

Last change on this file since 550 was 520, checked in by raasch, 14 years ago

last commit documented

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