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

Last change on this file since 836 was 836, checked in by raasch, 10 years ago

last commit documented

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