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

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

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:


Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:


Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

Bugfix: Rayleigh damping for ocean fixed.

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