source: palm/tags/release-3.2/SOURCE/advec_particles.f90 @ 3968

Last change on this file since 3968 was 77, checked in by raasch, 17 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

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