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

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

New:
---

The number of parallel I/O operations can be limited with new mrun-option -w.
(advec_particles, data_output_2d, data_output_3d, header, init_grid, init_pegrid, init_3d_model, modules, palm, parin, write_3d_binary)

Changed:


mrun option -T is obligatory

Errors:


Bugfix: No zero assignments to volume_flow_initial and volume_flow_area in
case of normal restart runs. (init_3d_model)

initialization of u_0, v_0. This is just to avoid access of uninitialized
memory in exchange_horiz_2d, which causes respective error messages
when the Intel thread checker (inspector) is used. (production_e)

Bugfix for ts limitation (prandtl_fluxes)

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