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

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

last commit documented

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