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

Last change on this file since 799 was 799, checked in by franke, 12 years ago

Implementation of Wang collision kernel and bugfix for calculation of vpt, pt_p, and ec in case of cloud droplets

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