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

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

thermal_conductivity_l and diff_coeff_l now depend on temperature and pressure

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