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

Last change on this file since 757 was 668, checked in by suehring, 14 years ago

last commit documented

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