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

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

last commit documented

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