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

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

New:
---
compare_palm_logs is additionally compiled with mbuild -u (Makefile in trunk/UTIL)

make options (mopts) to be set by configuration file implemented (mrun, mbuild)

humidity=.T. is now usable for runs with topography. wall_humidityflux and
wall_scalarflux are the corresponding new parin arrays.
(check_parameters, init_3d_model, parin)

Large scale vertical motion (subsidence/ascent) can be applied to the
prognostic equation for the potential temperature. (check_parameters, header,
Makefile, modules, parin, prognostic_equations, read_var_list, subsidence,
write_var_list)

A simple method for installing and running palm (with limited features)
has been added. (Makefile, palm_simple_install, palm_simple_run)

Changed:


2d-decomposition is default for Cray-XT machines (init_pegrid)

var_ts is replaced by dots_max (modules,init_3d_model)

Every cloud droplet has now an own weighting factor and can be deleted due to
collisions. Condensation and collision of cloud droplets are adjusted
accordingly. (advec_particles)

Collision efficiency for large cloud droplets has changed according to table of
Rogers and Yau. (collision_efficiency)

Errors:


Bugfix for generating serial jobs (subjob)

Bugfix: index problem concerning gradient_level indices removed (header)

Dimension of array stat in cascade change to prevent type problems with
mpi2 libraries (poisfft_hybrid)

Loop was split to make runs reproducible when using ifort compiler.
(disturb_field)

Bugfix: exchange of ghost points for prho included (time_integration)

Bugfix: calculation of time-averaged surface heatfluxes (sum_up_3d_data)

Bugfix: calculation of precipitation_rate (calc_precipitation)

Bugfix: initial data assignments to some dvrp arrays changed due to error
messages from gfortran compiler (modules)

Bugfix: calculation of cloud droplet velocity (advec_particles)

Bugfix: transfer of particles at south/left edge (advec_particles)

Bugfix: calculation of collision_efficiency (collision_efficiency)

Bugfix: initialisation of var_mod (subsidence)

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