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

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

last commit documented

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