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

Last change on this file since 336 was 336, checked in by raasch, 12 years ago

several small bugfixes; some more dvrp changes

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