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

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

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:


Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:


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