Ignore:
Timestamp:
Mar 15, 2012 10:35:09 AM (10 years ago)
Author:
raasch
Message:

Changed:


Original routine advec_particles split into several new subroutines and renamed
lpm.
init_particles renamed lpm_init
user_advec_particles renamed user_lpm_advec,
particle_boundary_conds renamed lpm_boundary_conds,
set_particle_attributes renamed lpm_set_attributes,
user_init_particles renamed user_lpm_init,
user_particle_attributes renamed user_lpm_set_attributes
(Makefile, lpm_droplet_collision, lpm_droplet_condensation, init_3d_model, modules, palm, read_var_list, time_integration, write_var_list, deleted: advec_particles, init_particles, particle_boundary_conds, set_particle_attributes, user_advec_particles, user_init_particles, user_particle_attributes, new: lpm, lpm_advec, lpm_boundary_conds, lpm_calc_liquid_water_content, lpm_data_output_particles, lpm_droplet_collision, lpm_drollet_condensation, lpm_exchange_horiz, lpm_extend_particle_array, lpm_extend_tails, lpm_extend_tail_array, lpm_init, lpm_init_sgs_tke, lpm_pack_arrays, lpm_read_restart_file, lpm_release_set, lpm_set_attributes, lpm_sort_arrays, lpm_write_exchange_statistics, lpm_write_restart_file, user_lpm_advec, user_lpm_init, user_lpm_set_attributes

File:
1 moved

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/lpm.f90

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