Index: /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4120)
+++ /palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4121)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Implementation of an simple method for interpolating the velocities to
+! particle position
+! 
+! 4114 2019-07-23 14:09:27Z schwenkel
 ! Bugfix: Added working precision for if statement
 ! 
@@ -72,5 +76,5 @@
 ! with scalar grid point of same index.
 ! Renamed module and subroutines: lpm_pack_arrays_mod -> lpm_pack_and_sort_mod
-! lpm_pack_all_arrays -> lpm_sort_in_subboxes, lpm_pack_arrays -> lpm_pack
+! lpm_pack_all_arrays -> lpm_sort_and_delete, lpm_pack_arrays -> lpm_pack
 ! lpm_sort -> lpm_sort_timeloop_done
 ! 
@@ -169,5 +173,5 @@
 ! Description:
 ! ------------
-!> 
+!>
 !------------------------------------------------------------------------------!
  MODULE lagrangian_particle_model_mod
@@ -178,5 +182,5 @@
         ONLY:  de_dx, de_dy, de_dz, dzw, zu, zw,  ql_c, ql_v, ql_vp, hyp,      &
                pt, q, exner, ql, diss, e, u, v, w, km, ql_1, ql_2, pt_p, q_p,  &
-               d_exner
+               d_exner, u_p, v_p, w_p
  
     USE averaging,                                                             &
@@ -263,14 +267,16 @@
     IMPLICIT NONE
 
-    CHARACTER(LEN=15) ::  aero_species = 'nacl'                    !< aerosol species
-    CHARACTER(LEN=15) ::  aero_type    = 'maritime'                !< aerosol type
-    CHARACTER(LEN=15) ::  bc_par_lr    = 'cyclic'                  !< left/right boundary condition
-    CHARACTER(LEN=15) ::  bc_par_ns    = 'cyclic'                  !< north/south boundary condition
-    CHARACTER(LEN=15) ::  bc_par_b     = 'reflect'                 !< bottom boundary condition
-    CHARACTER(LEN=15) ::  bc_par_t     = 'absorb'                  !< top boundary condition
-    CHARACTER(LEN=15) ::  collision_kernel   = 'none'              !< collision kernel    
-
-    CHARACTER(LEN=5)  ::  splitting_function = 'gamma'             !< function for calculation critical weighting factor
-    CHARACTER(LEN=5)  ::  splitting_mode     = 'const'             !< splitting mode
+    CHARACTER(LEN=15) ::  aero_species = 'nacl'                   !< aerosol species
+    CHARACTER(LEN=15) ::  aero_type    = 'maritime'               !< aerosol type
+    CHARACTER(LEN=15) ::  bc_par_lr    = 'cyclic'                 !< left/right boundary condition
+    CHARACTER(LEN=15) ::  bc_par_ns    = 'cyclic'                 !< north/south boundary condition
+    CHARACTER(LEN=15) ::  bc_par_b     = 'reflect'                !< bottom boundary condition
+    CHARACTER(LEN=15) ::  bc_par_t     = 'absorb'                 !< top boundary condition
+    CHARACTER(LEN=15) ::  collision_kernel   = 'none'             !< collision kernel
+
+    CHARACTER(LEN=5)  ::  splitting_function = 'gamma'            !< function for calculation critical weighting factor
+    CHARACTER(LEN=5)  ::  splitting_mode     = 'const'            !< splitting mode
+
+    CHARACTER(LEN=25) ::  particle_interpolation = 'trilinear'    !< interpolation method for calculatin the particle
 
     INTEGER(iwp) ::  deleted_particles = 0                        !< number of deleted particles per time step    
@@ -496,7 +502,7 @@
     END INTERFACE dealloc_particles_array
 
-    INTERFACE lpm_sort_in_subboxes
-       MODULE PROCEDURE lpm_sort_in_subboxes
-    END INTERFACE lpm_sort_in_subboxes
+    INTERFACE lpm_sort_and_delete
+       MODULE PROCEDURE lpm_sort_and_delete
+    END INTERFACE lpm_sort_and_delete
 
     INTERFACE lpm_sort_timeloop_done
@@ -549,4 +555,5 @@
        particles_per_point, &
        particle_advection_start, &
+       particle_interpolation, &
        particle_maximum_age, &
        pdx, &
@@ -608,4 +615,5 @@
        particles_per_point, &
        particle_advection_start, &
+       particle_interpolation, &
        particle_maximum_age, &
        pdx, &
@@ -849,4 +857,17 @@
     IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
 
+!
+!-- Subgrid scale velocites with the simple interpolation method for resolved
+!-- velocites is not implemented for passive particles. However, for cloud
+!-- it can be combined as the sgs-velocites for active particles are
+!-- calculated differently, i.e. no subboxes are needed.
+    IF ( .NOT. TRIM(particle_interpolation) == 'trilinear'  .AND.              &
+       use_sgs_for_particles .AND.  .NOT. cloud_droplets )  THEN
+          message_string = 'subrgrid scale velocities in combination with ' // &
+                           'simple interpolation method is not '            // &
+                           'implemented'
+          CALL message( 'check_parameters', 'PA0659', 1, 2, 0, 6, 0 )
+    ENDIF
+
  END SUBROUTINE
  
@@ -1067,5 +1088,5 @@
        CASE ( 'reflect' )
           ibc_par_t = 2
-          
+
        CASE ( 'nested' )
           ibc_par_t = 3
@@ -1087,5 +1108,5 @@
        CASE ( 'reflect' )
           ibc_par_lr = 2
-          
+
        CASE ( 'nested' )
           ibc_par_lr = 3
@@ -1107,5 +1128,5 @@
        CASE ( 'reflect' )
           ibc_par_ns = 2
-          
+
        CASE ( 'nested' )
           ibc_par_ns = 3
@@ -1546,5 +1567,5 @@
 !-- which is needed for a fast interpolation of the LES fields on the particle
 !-- position.
-    CALL lpm_sort_in_subboxes
+    CALL lpm_sort_and_delete
 !
 !-- Determine the current number of particles
@@ -1623,30 +1644,30 @@
 !-- Tropospheric aerosols. Published in Aerosol-Cloud-Climate Interactions.
     IF ( TRIM(aero_type) .EQ. 'polar' )  THEN
-       na        = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6
-       rm        = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6
+       na        = (/ 2.17e1, 1.86e-1, 3.04e-4 /) * 1.0E6_wp
+       rm        = (/ 0.0689, 0.375, 4.29 /) * 1.0E-6_wp
        log_sigma = (/ 0.245, 0.300, 0.291 /)
     ELSEIF ( TRIM(aero_type) .EQ. 'background' )  THEN
-       na        = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6
-       rm        = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6
+       na        = (/ 1.29e2, 5.97e1, 6.35e1 /) * 1.0E6_wp
+       rm        = (/ 0.0036, 0.127, 0.259 /) * 1.0E-6_wp
        log_sigma = (/ 0.645, 0.253, 0.425 /)
     ELSEIF ( TRIM(aero_type) .EQ. 'maritime' )  THEN
-       na        = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6
-       rm        = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6
+       na        = (/ 1.33e2, 6.66e1, 3.06e0 /) * 1.0E6_wp
+       rm        = (/ 0.0039, 0.133, 0.29 /) * 1.0E-6_wp
        log_sigma = (/ 0.657, 0.210, 0.396 /)
     ELSEIF ( TRIM(aero_type) .EQ. 'continental' )  THEN
-       na        = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6
-       rm        = (/ 0.01, 0.058, 0.9 /) * 1.0E-6
+       na        = (/ 3.20e3, 2.90e3, 3.00e-1 /) * 1.0E6_wp
+       rm        = (/ 0.01, 0.058, 0.9 /) * 1.0E-6_wp
        log_sigma = (/ 0.161, 0.217, 0.380 /)
     ELSEIF ( TRIM(aero_type) .EQ. 'desert' )  THEN
-       na        = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6
-       rm        = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6
+       na        = (/ 7.26e2, 1.14e3, 1.78e-1 /) * 1.0E6_wp
+       rm        = (/ 0.001, 0.0188, 10.8 /) * 1.0E-6_wp
        log_sigma = (/ 0.247, 0.770, 0.438 /)
     ELSEIF ( TRIM(aero_type) .EQ. 'rural' )  THEN
-       na        = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6
-       rm        = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6
+       na        = (/ 6.65e3, 1.47e2, 1.99e3 /) * 1.0E6_wp
+       rm        = (/ 0.00739, 0.0269, 0.0419 /) * 1.0E-6_wp
        log_sigma = (/ 0.225, 0.557, 0.266 /)
     ELSEIF ( TRIM(aero_type) .EQ. 'urban' )  THEN
-       na        = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6
-       rm        = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6
+       na        = (/ 9.93e4, 1.11e3, 3.64e4 /) * 1.0E6_wp
+       rm        = (/ 0.00651, 0.00714, 0.0248 /) * 1.0E-6_wp
        log_sigma = (/ 0.245, 0.666, 0.337 /)
     ELSEIF ( TRIM(aero_type) .EQ. 'user' )  THEN
@@ -1679,10 +1700,10 @@
                 particles(n)%aux1          = r_mid
                 particles(n)%weight_factor =                                           &
-                   ( na(1) / ( SQRT( 2.0 * pi ) * log_sigma(1) ) *                     &
-                     EXP( - LOG10( r_mid / rm(1) )**2 / ( 2.0 * log_sigma(1)**2 ) ) +  &
-                     na(2) / ( SQRT( 2.0 * pi ) * log_sigma(2) ) *                     &
-                     EXP( - LOG10( r_mid / rm(2) )**2 / ( 2.0 * log_sigma(2)**2 ) ) +  &
-                     na(3) / ( SQRT( 2.0 * pi ) * log_sigma(3) ) *                     &
-                     EXP( - LOG10( r_mid / rm(3) )**2 / ( 2.0 * log_sigma(3)**2 ) )    &
+                   ( na(1) / ( SQRT( 2.0_wp * pi ) * log_sigma(1) ) *                     &
+                     EXP( - LOG10( r_mid / rm(1) )**2 / ( 2.0_wp * log_sigma(1)**2 ) ) +  &
+                     na(2) / ( SQRT( 2.0_wp * pi ) * log_sigma(2) ) *                     &
+                     EXP( - LOG10( r_mid / rm(2) )**2 / ( 2.0_wp * log_sigma(2)**2 ) ) +  &
+                     na(3) / ( SQRT( 2.0_wp * pi ) * log_sigma(3) ) *                     &
+                     EXP( - LOG10( r_mid / rm(3) )**2 / ( 2.0_wp * log_sigma(3)**2 ) )    &
                    ) * ( LOG10(r_r) - LOG10(r_l) ) * ( dx * dy * dzw(kp) )
 
@@ -1700,5 +1721,5 @@
 !
 !--             Unnecessary particles will be deleted
-                IF ( particles(n)%weight_factor .LE. 0.0 )  particles(n)%particle_mask = .FALSE.
+                IF ( particles(n)%weight_factor .LE. 0.0_wp )  particles(n)%particle_mask = .FALSE.
 
              ENDDO
@@ -1741,6 +1762,6 @@
 
  END SUBROUTINE lpm_init_aerosols
- 
- 
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -1752,5 +1773,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lpm_init_sgs_tke
-    
+
     USE statistics,                                                            &
         ONLY:  flow_statistics_called, hom, sums, sums_l
@@ -2013,5 +2034,5 @@
     LOGICAL            ::  first_loop_stride  !<
 
-    
+
     SELECT CASE ( location )
 
@@ -2098,5 +2119,5 @@
              CALL cpu_log( log_point_s(44), 'lpm_advec', 'start' )
              CALL cpu_log( log_point_s(44), 'lpm_advec', 'pause' )
-             
+
 !
 !--          If particle advection includes SGS velocity components, calculate the
@@ -2116,5 +2137,4 @@
              IF ( .NOT. first_loop_stride  .AND.  use_sgs_for_particles )            &
                 CALL lpm_sort_timeloop_done
-
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -2164,5 +2184,5 @@
 !
 !--                   Particle advection
-                      CALL lpm_advec(i,j,k)
+                      CALL lpm_advec(TRIM(particle_interpolation),i,j,k)
 !
 !--                   Particle reflection from walls. Only applied if the particles 
@@ -2203,5 +2223,4 @@
                 ENDDO
              ENDDO
-
              steps = steps + 1
              dt_3d_reached_l = ALL(grid_particles(:,:,:)%time_loop_done)
@@ -2216,5 +2235,4 @@
              dt_3d_reached = dt_3d_reached_l
 #endif
-
              CALL cpu_log( log_point_s(44), 'lpm_advec', 'stop' )
 
@@ -2237,11 +2255,11 @@
 
 !
-!--          IF .FALSE., lpm_sort_in_subboxes is done inside pcmp
+!--          IF .FALSE., lpm_sort_and_delete is done inside pcmp
              IF ( .NOT. dt_3d_reached .OR. .NOT. nested_run )   THEN   
 !
 !--             Pack particles (eliminate those marked for deletion),
 !--             determine new number of particles
-                CALL lpm_sort_in_subboxes
-!
+                CALL lpm_sort_and_delete
+
 !--             Initialize variables for the next (sub-) timestep, i.e., for marking
 !--             those particles to be deleted after the timestep
@@ -2253,5 +2271,5 @@
              first_loop_stride = .FALSE.
           ENDDO   ! timestep loop
-!    
+!
 !--       in case of nested runs do the transfer of particles after every full model time step
           IF ( nested_run )   THEN
@@ -2260,5 +2278,5 @@
              CALL pmcp_p_delete_particles_in_fine_grid_area
 
-             CALL lpm_sort_in_subboxes
+             CALL lpm_sort_and_delete
 
              deleted_particles = 0
@@ -2284,5 +2302,5 @@
 
           CALL cpu_log( log_point(25), 'lpm', 'stop' )
-          
+
 ! !
 ! !--       Output of particle time series
@@ -2295,5 +2313,5 @@
 !              ENDIF
 !           ENDIF
-    
+
        CASE DEFAULT
           CONTINUE
@@ -2940,6 +2958,6 @@
 !-- Must be called to sort particles into blocks, which is needed for a fast
 !-- interpolation of the LES fields on the particle position.
-    CALL lpm_sort_in_subboxes
-       
+    CALL lpm_sort_and_delete
+
 
  END SUBROUTINE lpm_rrd_local_particles
@@ -2949,8 +2967,8 @@
                               nxr_on_file, nynf, nync, nyn_on_file, nysf,  &
                               nysc, nys_on_file, tmp_3d, found )
-                              
-       
+
+
    USE control_parameters,                                                 &
-       ONLY: length, restart_string           
+       ONLY: length, restart_string
 
     INTEGER(iwp) ::  k               !<
@@ -2972,12 +2990,12 @@
     REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
 
-    
+
     found = .TRUE.
 
     SELECT CASE ( restart_string(1:length) )
-    
+
        CASE ( 'iran' ) ! matching random numbers is still unresolved issue
           IF ( k == 1 )  READ ( 13 )  iran, iran_part
-          
+
         CASE ( 'pc_av' )
            IF ( .NOT. ALLOCATED( pc_av ) )  THEN
@@ -3154,27 +3172,41 @@
 !> grid scale) as well as the boundary conditions of particles. As a next step
 !> this submodule should be excluded as an own file. 
-!------------------------------------------------------------------------------!   
- SUBROUTINE lpm_advec (ip,jp,kp)
-
+!------------------------------------------------------------------------------!
+ SUBROUTINE lpm_advec (interpolation_method,ip,jp,kp)
+
+    CHARACTER (LEN=*), INTENT(IN) ::  interpolation_method !<
     LOGICAL ::  subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall
 
     INTEGER(iwp) ::  i                           !< index variable along x
-    INTEGER(iwp) ::  ip                          !< index variable along x 
-    INTEGER(iwp) ::  j                           !< index variable along y 
-    INTEGER(iwp) ::  jp                          !< index variable along y 
-    INTEGER(iwp) ::  k                           !< index variable along z 
+    INTEGER(iwp) ::  i_next                      !< index variable along x
+    INTEGER(iwp) ::  ip                          !< index variable along x
+    INTEGER(iwp) ::  iteration_steps = 1         !< amount of iterations steps for corrector step
+    INTEGER(iwp) ::  j                           !< index variable along y
+    INTEGER(iwp) ::  j_next                      !< index variable along y
+    INTEGER(iwp) ::  jp                          !< index variable along y
+    INTEGER(iwp) ::  k                           !< index variable along z
     INTEGER(iwp) ::  k_wall                      !< vertical index of topography top 
-    INTEGER(iwp) ::  kp                          !< index variable along z 
+    INTEGER(iwp) ::  kp                          !< index variable along z
+    INTEGER(iwp) ::  k_next                      !< index variable along z
     INTEGER(iwp) ::  kw                          !< index variable along z
+    INTEGER(iwp) ::  kkw                         !< index variable along z
     INTEGER(iwp) ::  n                           !< loop variable over all particles in a grid box
     INTEGER(iwp) ::  nb                          !< block number particles are sorted in
-    INTEGER(iwp) ::  surf_start                  !< Index on surface data-type for current grid box 
+    INTEGER(iwp) ::  particle_end                !< end index for partilce loop
+    INTEGER(iwp) ::  particle_start              !< start index for particle loop
+    INTEGER(iwp) ::  surf_start                  !< Index on surface data-type for current grid box
+    INTEGER(iwp) ::  subbox_end                  !< end index for loop over subboxes in particle advection
+    INTEGER(iwp) ::  subbox_start                !< start index for loop over subboxes in particle advection
+    INTEGER(iwp) ::  nn                          !< loop variable over iterations steps
 
     INTEGER(iwp), DIMENSION(0:7) ::  start_index !< start particle index for current block
     INTEGER(iwp), DIMENSION(0:7) ::  end_index   !< start particle index for current block
 
-    REAL(wp) ::  aa                 !< dummy argument for horizontal particle interpolation 
-    REAL(wp) ::  bb                 !< dummy argument for horizontal particle interpolation 
-    REAL(wp) ::  cc                 !< dummy argument for horizontal particle interpolation 
+    REAL(wp) ::  aa                 !< dummy argument for horizontal particle interpolation
+    REAL(wp) ::  alpha              !< interpolation facor for x-direction
+
+    REAL(wp) ::  bb                 !< dummy argument for horizontal particle interpolation
+    REAL(wp) ::  beta               !< interpolation facor for y-direction
+    REAL(wp) ::  cc                 !< dummy argument for horizontal particle interpolation
     REAL(wp) ::  d_z_p_z0           !< inverse of interpolation length for logarithmic interpolation 
     REAL(wp) ::  dd                 !< dummy argument for horizontal particle interpolation 
@@ -3197,4 +3229,5 @@
     REAL(wp) ::  exp_arg            !< argument in the exponent - particle radius
     REAL(wp) ::  exp_term           !< exponent term
+    REAL(wp) ::  gamma              !< interpolation facor for z-direction
     REAL(wp) ::  gg                 !< dummy argument for horizontal particle interpolation 
     REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
@@ -3208,4 +3241,5 @@
     REAL(wp) ::  u_int_l            !< x/y-interpolated u-component at particle position at lower vertical level
     REAL(wp) ::  u_int_u            !< x/y-interpolated u-component at particle position at upper vertical level
+    REAL(wp) ::  unext              !< calculated particle u-velocity of corrector step
     REAL(wp) ::  us_int             !< friction velocity at particle grid box
     REAL(wp) ::  usws_int           !< surface momentum flux (u component) at particle grid box
@@ -3213,11 +3247,16 @@
     REAL(wp) ::  v_int_u            !< x/y-interpolated v-component at particle position at upper vertical level
     REAL(wp) ::  vsws_int           !< surface momentum flux (u component) at particle grid box
+    REAL(wp) ::  vnext              !< calculated particle v-velocity of corrector step
     REAL(wp) ::  vv_int             !< dummy to compute interpolated mean SGS TKE, used to scale SGS advection 
     REAL(wp) ::  w_int_l            !< x/y-interpolated w-component at particle position at lower vertical level
     REAL(wp) ::  w_int_u            !< x/y-interpolated w-component at particle position at upper vertical level
+    REAL(wp) ::  wnext              !< calculated particle w-velocity of corrector step
     REAL(wp) ::  w_s                !< terminal velocity of droplets
     REAL(wp) ::  x                  !< dummy argument for horizontal particle interpolation 
-    REAL(wp) ::  y                  !< dummy argument for horizontal particle interpolation 
+    REAL(wp) ::  xp                 !< calculated particle position in x of predictor step
+    REAL(wp) ::  y                  !< dummy argument for horizontal particle interpolation
+    REAL(wp) ::  yp                 !< calculated particle position in y of predictor step
     REAL(wp) ::  z_p                !< surface layer height (0.5 dz)
+    REAL(wp) ::  zp                 !< calculated particle position in z of predictor step
 
     REAL(wp), PARAMETER ::  a_rog = 9.65_wp      !< parameter for fall velocity
@@ -3249,263 +3288,356 @@
     REAL(wp), DIMENSION(number_of_particles) ::  zv             !< z-position
 
-    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers 
+    REAL(wp), DIMENSION(number_of_particles, 3) ::  rg !< vector of Gaussian distributed random numbers
 
     CALL cpu_log( log_point_s(44), 'lpm_advec', 'continue' )
-
 !
 !-- Determine height of Prandtl layer and distance between Prandtl-layer
-!-- height and horizontal mean roughness height, which are required for 
-!-- vertical logarithmic interpolation of horizontal particle speeds 
+!-- height and horizontal mean roughness height, which are required for
+!-- vertical logarithmic interpolation of horizontal particle speeds
 !-- (for particles below first vertical grid level).
     z_p      = zu(nzb+1) - zw(nzb)
     d_z_p_z0 = 1.0_wp / ( z_p - z0_av_global )
 
-    start_index = grid_particles(kp,jp,ip)%start_index
-    end_index   = grid_particles(kp,jp,ip)%end_index
-
     xv = particles(1:number_of_particles)%x
     yv = particles(1:number_of_particles)%y
     zv = particles(1:number_of_particles)%z
-
-    DO  nb = 0, 7
-!
-!--    Interpolate u velocity-component       
-       i = ip
-       j = jp + block_offset(nb)%j_off
-       k = kp + block_offset(nb)%k_off
-
-       DO  n = start_index(nb), end_index(nb)
-!
-!--       Interpolation of the u velocity component onto particle position.  
-!--       Particles are interpolation bi-linearly in the horizontal and a 
-!--       linearly in the vertical. An exception is made for particles below
-!--       the first vertical grid level in case of a prandtl layer. In this 
-!--       case the horizontal particle velocity components are determined using
-!--       Monin-Obukhov relations (if branch). 
-!--       First, check if particle is located below first vertical grid level 
-!--       above topography (Prandtl-layer height)
-!--       Determine vertical index of topography top
-          k_wall = get_topography_top_index_ji( jp, ip, 's' )
-
-          IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
-!
-!--          Resolved-scale horizontal particle velocity is zero below z0. 
-             IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
-                u_int(n) = 0.0_wp
-             ELSE
-!
-!--             Determine the sublayer. Further used as index.
-                height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
-                                     * REAL( number_of_sublayers, KIND=wp )    &
-                                     * d_z_p_z0 
-!
-!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
-!--             interpolate linearly between precalculated logarithm. 
-                log_z_z0_int = log_z_z0(INT(height_p))                         &
-                                 + ( height_p - INT(height_p) )                &
-                                 * ( log_z_z0(INT(height_p)+1)                 &
-                                      - log_z_z0(INT(height_p))                &
-                                   ) 
-!
-!--             Get friction velocity and momentum flux from new surface data 
-!--             types.
-                IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
-                     surf_def_h(0)%end_index(jp,ip) )  THEN
-                   surf_start = surf_def_h(0)%start_index(jp,ip)
-!--                Limit friction velocity. In narrow canyons or holes the 
-!--                friction velocity can become very small, resulting in a too
-!--                large particle speed.
-                   us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )  
-                   usws_int  = surf_def_h(0)%usws(surf_start)
-                ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
-                         surf_lsm_h%end_index(jp,ip) )  THEN
-                   surf_start = surf_lsm_h%start_index(jp,ip)
-                   us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )  
-                   usws_int  = surf_lsm_h%usws(surf_start)
-                ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
-                         surf_usm_h%end_index(jp,ip) )  THEN
-                   surf_start = surf_usm_h%start_index(jp,ip)
-                   us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )  
-                   usws_int  = surf_usm_h%usws(surf_start)
+    dt_particle = dt_3d
+
+
+    SELECT CASE ( interpolation_method )
+
+!
+!--    This case uses a simple interpolation method for the particle velocites,
+!--    and applying a predictor-corrector method. @attention: for the corrector
+!--    step the velocities of t(n+1) are required. However, at this moment of
+!--    the time integration they are not free of divergence. This interpolation
+!--    method is described in more detail in Grabowski et al., 2018 (GMD).
+       CASE ( 'simple_corrector' )
+!
+!--       Predictor step
+          kkw = kp - 1
+          DO n = 1, number_of_particles
+
+             alpha = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
+             u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha
+
+             beta  = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
+             v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta
+
+             gamma = MAX( MIN( ( particles(n)%z - zw(kkw) ) /                   &
+                               ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp )
+             w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
+
+          ENDDO
+
+!
+!--       Corrector step
+          DO n = 1, number_of_particles
+
+             IF ( .NOT. particles(n)%particle_mask ) CYCLE
+
+             DO nn = 1, iteration_steps
+
+!
+!--             Guess new position
+                xp = particles(n)%x + u_int(n) * dt_particle(n)
+                yp = particles(n)%y + v_int(n) * dt_particle(n)
+                zp = particles(n)%z + w_int(n) * dt_particle(n)
+!
+!--             x direction
+                i_next = FLOOR( xp * ddx , KIND=iwp)
+                alpha  = MAX( MIN( ( xp - i_next * dx ) * ddx, 1.0_wp ), 0.0_wp )
+!
+!--             y direction
+                j_next = FLOOR( yp * ddy )
+                beta   = MAX( MIN( ( yp - j_next * dy ) * ddy, 1.0_wp ), 0.0_wp )
+!
+!--             z_direction
+                k_next = MAX( MIN( FLOOR( zp / (zw(kkw+1)-zw(kkw)) ), nzt ), 0)
+                gamma = MAX( MIN( ( zp - zw(k_next) ) /                      &
+                                  ( zw(k_next+1) - zw(k_next) ), 1.0_wp ), 0.0_wp )
+!
+!--             Calculate part of the corrector step
+                unext = u_p(k_next+1, j_next, i_next) * ( 1.0_wp - alpha ) +    &
+                        u_p(k_next+1, j_next,   i_next+1) * alpha
+
+                vnext = v_p(k_next+1, j_next, i_next) * ( 1.0_wp - beta  ) +    &
+                        v_p(k_next+1, j_next+1, i_next  ) * beta
+
+                wnext = w_p(k_next,   j_next, i_next) * ( 1.0_wp - gamma ) +    &
+                        w_p(k_next+1, j_next, i_next  ) * gamma
+
+!
+!--             Calculate interpolated particle velocity with predictor
+!--             corrector step. u_int, v_int and w_int describes the part of
+!--             the predictor step. unext, vnext and wnext is the part of the
+!--             corrector step. The resulting new position is set below. The
+!--             implementation is based on Grabowski et al., 2018 (GMD).
+                u_int(n) = 0.5_wp * ( u_int(n) + unext )
+                v_int(n) = 0.5_wp * ( v_int(n) + vnext )
+                w_int(n) = 0.5_wp * ( w_int(n) + wnext )
+
+             ENDDO
+          ENDDO
+
+
+!
+!--    This case uses a simple interpolation method for the particle velocites,
+!--    and applying a predictor.
+       CASE ( 'simple_predictor' )
+!
+!--       The particle position for the w velociy is based on the value of kp and kp-1
+          kkw = kp - 1
+          DO n = 1, number_of_particles
+             IF ( .NOT. particles(n)%particle_mask ) CYCLE
+
+             alpha    = MAX( MIN( ( particles(n)%x - ip * dx ) * ddx, 1.0_wp ), 0.0_wp )
+             u_int(n) = u(kp,jp,ip) * ( 1.0_wp - alpha ) + u(kp,jp,ip+1) * alpha
+
+             beta     = MAX( MIN( ( particles(n)%y - jp * dy ) * ddy, 1.0_wp ), 0.0_wp )
+             v_int(n) = v(kp,jp,ip) * ( 1.0_wp - beta ) + v(kp,jp+1,ip) * beta
+
+             gamma    = MAX( MIN( ( particles(n)%z - zw(kkw) ) /                   &
+                                  ( zw(kkw+1) - zw(kkw) ), 1.0_wp ), 0.0_wp )
+             w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
+          ENDDO
+!
+!--    The trilinear interpolation.
+       CASE ( 'trilinear' )
+
+          start_index = grid_particles(kp,jp,ip)%start_index
+          end_index   = grid_particles(kp,jp,ip)%end_index
+
+          DO  nb = 0, 7
+!
+!--          Interpolate u velocity-component
+             i = ip
+             j = jp + block_offset(nb)%j_off
+             k = kp + block_offset(nb)%k_off
+
+             DO  n = start_index(nb), end_index(nb)
+!
+!--             Interpolation of the u velocity component onto particle position.
+!--             Particles are interpolation bi-linearly in the horizontal and a
+!--             linearly in the vertical. An exception is made for particles below
+!--             the first vertical grid level in case of a prandtl layer. In this
+!--             case the horizontal particle velocity components are determined using
+!--             Monin-Obukhov relations (if branch).
+!--             First, check if particle is located below first vertical grid level
+!--             above topography (Prandtl-layer height)
+!--             Determine vertical index of topography top
+                k_wall = get_topography_top_index_ji( jp, ip, 's' )
+
+                IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
+!
+!--                Resolved-scale horizontal particle velocity is zero below z0.
+                   IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
+                      u_int(n) = 0.0_wp
+                   ELSE
+!
+!--                   Determine the sublayer. Further used as index.
+                      height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
+                                           * REAL( number_of_sublayers, KIND=wp )    &
+                                           * d_z_p_z0
+!
+!--                   Calculate LOG(z/z0) for exact particle height. Therefore,
+!--                   interpolate linearly between precalculated logarithm.
+                      log_z_z0_int = log_z_z0(INT(height_p))                         &
+                                       + ( height_p - INT(height_p) )                &
+                                       * ( log_z_z0(INT(height_p)+1)                 &
+                                            - log_z_z0(INT(height_p))                &
+                                         )
+!
+!--                   Get friction velocity and momentum flux from new surface data
+!--                   types.
+                      IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
+                           surf_def_h(0)%end_index(jp,ip) )  THEN
+                         surf_start = surf_def_h(0)%start_index(jp,ip)
+!--                      Limit friction velocity. In narrow canyons or holes the
+!--                      friction velocity can become very small, resulting in a too
+!--                      large particle speed.
+                         us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
+                         usws_int  = surf_def_h(0)%usws(surf_start)
+                      ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
+                               surf_lsm_h%end_index(jp,ip) )  THEN
+                         surf_start = surf_lsm_h%start_index(jp,ip)
+                         us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
+                         usws_int  = surf_lsm_h%usws(surf_start)
+                      ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
+                               surf_usm_h%end_index(jp,ip) )  THEN
+                         surf_start = surf_usm_h%start_index(jp,ip)
+                         us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
+                         usws_int  = surf_usm_h%usws(surf_start)
+                      ENDIF
+!
+!--                   Neutral solution is applied for all situations, e.g. also for
+!--                   unstable and stable situations. Even though this is not exact
+!--                   this saves a lot of CPU time since several calls of intrinsic
+!--                   FORTRAN procedures (LOG, ATAN) are avoided, This is justified
+!--                   as sensitivity studies revealed no significant effect of
+!--                   using the neutral solution also for un/stable situations.
+                      u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )           &
+                                  * log_z_z0_int - u_gtrans
+                   ENDIF
+!
+!--             Particle above the first grid level. Bi-linear interpolation in the
+!--             horizontal and linear interpolation in the vertical direction.
+                ELSE
+                   x  = xv(n) - i * dx
+                   y  = yv(n) + ( 0.5_wp - j ) * dy
+                   aa = x**2          + y**2
+                   bb = ( dx - x )**2 + y**2
+                   cc = x**2          + ( dy - y )**2
+                   dd = ( dx - x )**2 + ( dy - y )**2
+                   gg = aa + bb + cc + dd
+
+                   u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
+                               + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
+                               u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
+
+                   IF ( k == nzt )  THEN
+                      u_int(n) = u_int_l
+                   ELSE
+                      u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
+                                  + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
+                                  u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
+                      u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *            &
+                                 ( u_int_u - u_int_l )
+                   ENDIF
                 ENDIF
-
-!
-!--             Neutral solution is applied for all situations, e.g. also for 
-!--             unstable and stable situations. Even though this is not exact 
-!--             this saves a lot of CPU time since several calls of intrinsic 
-!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified 
-!--             as sensitivity studies revealed no significant effect of 
-!--             using the neutral solution also for un/stable situations.
-                u_int(n) = -usws_int / ( us_int * kappa + 1E-10_wp )           & 
-                            * log_z_z0_int - u_gtrans
-
-             ENDIF
-!
-!--       Particle above the first grid level. Bi-linear interpolation in the 
-!--       horizontal and linear interpolation in the vertical direction.
-          ELSE
-
-             x  = xv(n) - i * dx
-             y  = yv(n) + ( 0.5_wp - j ) * dy
-             aa = x**2          + y**2
-             bb = ( dx - x )**2 + y**2
-             cc = x**2          + ( dy - y )**2
-             dd = ( dx - x )**2 + ( dy - y )**2
-             gg = aa + bb + cc + dd
-
-             u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
-                         + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) *            &
-                         u(k,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
-
-             IF ( k == nzt )  THEN
-                u_int(n) = u_int_l
-             ELSE
-                u_int_u = ( ( gg-aa ) * u(k+1,j,i) + ( gg-bb ) * u(k+1,j,i+1)  &
-                            + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) *           &
-                            u(k+1,j+1,i+1) ) / ( 3.0_wp * gg ) - u_gtrans
-                u_int(n) = u_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *               &
-                           ( u_int_u - u_int_l )
-             ENDIF
-
-          ENDIF
-
-       ENDDO
-!
-!--    Same procedure for interpolation of the v velocity-component 
-       i = ip + block_offset(nb)%i_off
-       j = jp
-       k = kp + block_offset(nb)%k_off
-
-       DO  n = start_index(nb), end_index(nb)
-
-!
-!--       Determine vertical index of topography top
-          k_wall = get_topography_top_index_ji( jp,ip, 's' )
-
-          IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
-             IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
-!
-!--             Resolved-scale horizontal particle velocity is zero below z0. 
-                v_int(n) = 0.0_wp
-             ELSE       
-!
-!--             Determine the sublayer. Further used as index. Please note, 
-!--             logarithmus can not be reused from above, as in in case of
-!--             topography particle on u-grid can be above surface-layer height,
-!--             whereas it can be below on v-grid. 
-                height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
-                                  * REAL( number_of_sublayers, KIND=wp )       &
-                                  * d_z_p_z0 
-!
-!--             Calculate LOG(z/z0) for exact particle height. Therefore,   
-!--             interpolate linearly between precalculated logarithm. 
-                log_z_z0_int = log_z_z0(INT(height_p))                         &
-                                 + ( height_p - INT(height_p) )                &
-                                 * ( log_z_z0(INT(height_p)+1)                 &
-                                      - log_z_z0(INT(height_p))                &
-                                   ) 
-!
-!--             Get friction velocity and momentum flux from new surface data 
-!--             types.
-                IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
-                     surf_def_h(0)%end_index(jp,ip) )  THEN
-                   surf_start = surf_def_h(0)%start_index(jp,ip)
-!--                Limit friction velocity. In narrow canyons or holes the 
-!--                friction velocity can become very small, resulting in a too
-!--                large particle speed.
-                   us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )  
-                   vsws_int  = surf_def_h(0)%vsws(surf_start)
-                ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
-                         surf_lsm_h%end_index(jp,ip) )  THEN
-                   surf_start = surf_lsm_h%start_index(jp,ip)
-                   us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )  
-                   vsws_int  = surf_lsm_h%vsws(surf_start)
-                ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
-                         surf_usm_h%end_index(jp,ip) )  THEN
-                   surf_start = surf_usm_h%start_index(jp,ip)
-                   us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )  
-                   vsws_int  = surf_usm_h%vsws(surf_start)
-                ENDIF 
-!
-!--             Neutral solution is applied for all situations, e.g. also for 
-!--             unstable and stable situations. Even though this is not exact 
-!--             this saves a lot of CPU time since several calls of intrinsic 
-!--             FORTRAN procedures (LOG, ATAN) are avoided, This is justified 
-!--             as sensitivity studies revealed no significant effect of 
-!--             using the neutral solution also for un/stable situations.
-                v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )           &
-                         * log_z_z0_int - v_gtrans
-
-             ENDIF
-
-          ELSE
-             x  = xv(n) + ( 0.5_wp - i ) * dx
-             y  = yv(n) - j * dy
-             aa = x**2          + y**2
-             bb = ( dx - x )**2 + y**2
-             cc = x**2          + ( dy - y )**2
-             dd = ( dx - x )**2 + ( dy - y )**2
-             gg = aa + bb + cc + dd
-
-             v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
-                       + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
-                       ) / ( 3.0_wp * gg ) - v_gtrans
-
-             IF ( k == nzt )  THEN
-                v_int(n) = v_int_l
-             ELSE
-                v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
-                          + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
-                          ) / ( 3.0_wp * gg ) - v_gtrans
-                v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *               &
-                                  ( v_int_u - v_int_l )
-             ENDIF
-
-          ENDIF
-
-       ENDDO
-!
-!--    Same procedure for interpolation of the w velocity-component
-       i = ip + block_offset(nb)%i_off
-       j = jp + block_offset(nb)%j_off
-       k = kp - 1
-
-       DO  n = start_index(nb), end_index(nb)
-
-          IF ( vertical_particle_advection(particles(n)%group) )  THEN
-
-             x  = xv(n) + ( 0.5_wp - i ) * dx
-             y  = yv(n) + ( 0.5_wp - j ) * dy
-             aa = x**2          + y**2
-             bb = ( dx - x )**2 + y**2
-             cc = x**2          + ( dy - y )**2
-             dd = ( dx - x )**2 + ( dy - y )**2
-             gg = aa + bb + cc + dd
-
-             w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
-                       + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
-                       ) / ( 3.0_wp * gg )
-
-             IF ( k == nzt )  THEN
-                w_int(n) = w_int_l
-             ELSE
-                w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
-                            ( gg-bb ) * w(k+1,j,i+1) + &
-                            ( gg-cc ) * w(k+1,j+1,i) + &
-                            ( gg-dd ) * w(k+1,j+1,i+1) &
-                          ) / ( 3.0_wp * gg )
-                w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) *               &
-                           ( w_int_u - w_int_l )
-             ENDIF
-
-          ELSE
-
-             w_int(n) = 0.0_wp
-
-          ENDIF
-
-       ENDDO
-
-    ENDDO
+             ENDDO
+!
+!--          Same procedure for interpolation of the v velocity-component
+             i = ip + block_offset(nb)%i_off
+             j = jp
+             k = kp + block_offset(nb)%k_off
+
+             DO  n = start_index(nb), end_index(nb)
+!
+!--             Determine vertical index of topography top
+                k_wall = get_topography_top_index_ji( jp,ip, 's' )
+
+                IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
+                   IF ( zv(n) - zw(k_wall) < z0_av_global )  THEN
+!
+!--                   Resolved-scale horizontal particle velocity is zero below z0.
+                      v_int(n) = 0.0_wp
+                   ELSE
+!
+!--                   Determine the sublayer. Further used as index. Please note,
+!--                   logarithmus can not be reused from above, as in in case of
+!--                   topography particle on u-grid can be above surface-layer height,
+!--                   whereas it can be below on v-grid.
+                      height_p = ( zv(n) - zw(k_wall) - z0_av_global ) &
+                                        * REAL( number_of_sublayers, KIND=wp )       &
+                                        * d_z_p_z0
+!
+!--                   Calculate LOG(z/z0) for exact particle height. Therefore,
+!--                   interpolate linearly between precalculated logarithm.
+                      log_z_z0_int = log_z_z0(INT(height_p))                         &
+                                       + ( height_p - INT(height_p) )                &
+                                       * ( log_z_z0(INT(height_p)+1)                 &
+                                            - log_z_z0(INT(height_p))                &
+                                         )
+!
+!--                   Get friction velocity and momentum flux from new surface data
+!--                   types.
+                      IF ( surf_def_h(0)%start_index(jp,ip) <=                   &
+                           surf_def_h(0)%end_index(jp,ip) )  THEN
+                         surf_start = surf_def_h(0)%start_index(jp,ip)
+!--                      Limit friction velocity. In narrow canyons or holes the
+!--                      friction velocity can become very small, resulting in a too
+!--                      large particle speed.
+                         us_int    = MAX( surf_def_h(0)%us(surf_start), 0.01_wp )
+                         vsws_int  = surf_def_h(0)%vsws(surf_start)
+                      ELSEIF ( surf_lsm_h%start_index(jp,ip) <=                  &
+                               surf_lsm_h%end_index(jp,ip) )  THEN
+                         surf_start = surf_lsm_h%start_index(jp,ip)
+                         us_int    = MAX( surf_lsm_h%us(surf_start), 0.01_wp )
+                         vsws_int  = surf_lsm_h%vsws(surf_start)
+                      ELSEIF ( surf_usm_h%start_index(jp,ip) <=                  &
+                               surf_usm_h%end_index(jp,ip) )  THEN
+                         surf_start = surf_usm_h%start_index(jp,ip)
+                         us_int    = MAX( surf_usm_h%us(surf_start), 0.01_wp )
+                         vsws_int  = surf_usm_h%vsws(surf_start)
+                      ENDIF
+!
+!--                   Neutral solution is applied for all situations, e.g. also for
+!--                   unstable and stable situations. Even though this is not exact
+!--                   this saves a lot of CPU time since several calls of intrinsic
+!--                   FORTRAN procedures (LOG, ATAN) are avoided, This is justified
+!--                   as sensitivity studies revealed no significant effect of
+!--                   using the neutral solution also for un/stable situations.
+                      v_int(n) = -vsws_int / ( us_int * kappa + 1E-10_wp )           &
+                               * log_z_z0_int - v_gtrans
+
+                   ENDIF
+                ELSE
+                   x  = xv(n) + ( 0.5_wp - i ) * dx
+                   y  = yv(n) - j * dy
+                   aa = x**2          + y**2
+                   bb = ( dx - x )**2 + y**2
+                   cc = x**2          + ( dy - y )**2
+                   dd = ( dx - x )**2 + ( dy - y )**2
+                   gg = aa + bb + cc + dd
+
+                   v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
+                             + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
+                             ) / ( 3.0_wp * gg ) - v_gtrans
+
+                   IF ( k == nzt )  THEN
+                      v_int(n) = v_int_l
+                   ELSE
+                      v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
+                                + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
+                                ) / ( 3.0_wp * gg ) - v_gtrans
+                      v_int(n) = v_int_l + ( zv(n) - zu(k) ) / dzw(k+1) *               &
+                                        ( v_int_u - v_int_l )
+                   ENDIF
+                ENDIF
+             ENDDO
+!
+!--          Same procedure for interpolation of the w velocity-component
+             i = ip + block_offset(nb)%i_off
+             j = jp + block_offset(nb)%j_off
+             k = kp - 1
+
+             DO  n = start_index(nb), end_index(nb)
+                IF ( vertical_particle_advection(particles(n)%group) )  THEN
+                   x  = xv(n) + ( 0.5_wp - i ) * dx
+                   y  = yv(n) + ( 0.5_wp - j ) * dy
+                   aa = x**2          + y**2
+                   bb = ( dx - x )**2 + y**2
+                   cc = x**2          + ( dy - y )**2
+                   dd = ( dx - x )**2 + ( dy - y )**2
+                   gg = aa + bb + cc + dd
+
+                   w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
+                             + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
+                             ) / ( 3.0_wp * gg )
+
+                   IF ( k == nzt )  THEN
+                      w_int(n) = w_int_l
+                   ELSE
+                      w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + &
+                                  ( gg-bb ) * w(k+1,j,i+1) + &
+                                  ( gg-cc ) * w(k+1,j+1,i) + &
+                                  ( gg-dd ) * w(k+1,j+1,i+1) &
+                                ) / ( 3.0_wp * gg )
+                      w_int(n) = w_int_l + ( zv(n) - zw(k) ) / dzw(k+1) *               &
+                                 ( w_int_u - w_int_l )
+                   ENDIF
+                ELSE
+                   w_int(n) = 0.0_wp
+                ENDIF
+             ENDDO
+          ENDDO
+
+       CASE DEFAULT
+          WRITE( message_string, * )  'unknown particle velocity interpolation method = "',  &
+                                       TRIM( interpolation_method ), '"'
+          CALL message( 'lpm_advec', 'PA0660', 1, 2, 0, 6, 0 )
+
+    END SELECT
 
 !-- Interpolate and calculate quantities needed for calculating the SGS
@@ -3819,5 +3951,5 @@
 !--    and calculate SGS velocities again
        dz_temp = zw(kp)-zw(kp-1)
-       
+
        DO  nb = 0, 7
           DO  n = start_index(nb), end_index(nb)
@@ -3826,5 +3958,5 @@
                  ABS( v_int(n) + rvar2_temp(n) ) > (dy/dt_particle(n))  .OR.   &
                  ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n)))) THEN
-                
+
                 dt_particle(n) = 0.9_wp * MIN(                                 &
                                  ( dx / ABS( u_int(n) + rvar1_temp(n) ) ),     &
@@ -3837,5 +3969,5 @@
                 rvar2_temp(n) = particles(n)%rvar2
                 rvar3_temp(n) = particles(n)%rvar3
-                
+
                 de_dt_min = - e_int(n) / dt_particle(n)
 
@@ -3857,6 +3989,6 @@
                                         de_dz_int(n), de_dt, diss_int(n),       &
                                         dt_particle(n), rg(n,3), term_1_2(n) )
-             ENDIF                           
-             
+             ENDIF
+
 !
 !--          Update particle velocites 
@@ -3873,5 +4005,5 @@
           ENDDO
        ENDDO
-       
+
     ELSE
 !
@@ -3883,8 +4015,35 @@
 
     dens_ratio = particle_groups(particles(1:number_of_particles)%group)%density_ratio
-
     IF ( ANY( dens_ratio == 0.0_wp ) )  THEN
-       DO  nb = 0, 7
-          DO  n = start_index(nb), end_index(nb)
+!
+!--    Decide whether the particle loop runs over the subboxes or only over 1,
+!--    number_of_particles. This depends on the selected interpolation method.
+!--    If particle interpolation method is not trilinear, then the sorting within
+!--    subboxes is not required. However, therefore the index start_index(nb) and
+!--    end_index(nb) are not defined and the loops are still over
+!--    number_of_particles. @todo find a more generic way to write this loop or
+!--    delete trilinear interpolation
+       IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
+          subbox_start = 0
+          subbox_end   = 7
+       ELSE
+          subbox_start = 1
+          subbox_end   = 1
+       ENDIF
+!
+!--    loop over subboxes. In case of simple interpolation scheme no subboxes
+!--    are introduced, as they are not required. Accordingly, this loops goes
+!--    from 1 to 1.
+       DO  nb = subbox_start, subbox_end
+          IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
+             particle_start = start_index(nb)
+             particle_end   = end_index(nb)
+          ELSE
+             particle_start = 1
+             particle_end   = number_of_particles
+          ENDIF
+!
+!--         Loop from particle start to particle end
+            DO  n = particle_start, particle_end
 
 !
@@ -3915,5 +4074,5 @@
                 IF ( cloud_droplets )  THEN
 !
-!--                Terminal velocity is computed for vertical direction (Rogers et 
+!--                Terminal velocity is computed for vertical direction (Rogers et
 !--                al., 1993, J. Appl. Meteorol.)
                    diameter = particles(n)%radius * 2000.0_wp !diameter in mm
@@ -3974,9 +4133,31 @@
           ENDDO
        ENDDO
-    
+
     ELSE
-
-       DO  nb = 0, 7
-          DO  n = start_index(nb), end_index(nb)
+!
+!--    Decide whether the particle loop runs over the subboxes or only over 1,
+!--    number_of_particles. This depends on the selected interpolation method.
+       IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
+          subbox_start = 0
+          subbox_end   = 7
+       ELSE
+          subbox_start = 1
+          subbox_end   = 1
+       ENDIF
+!--    loop over subboxes. In case of simple interpolation scheme no subboxes
+!--    are introduced, as they are not required. Accordingly, this loops goes
+!--    from 1 to 1.
+       DO  nb = subbox_start, subbox_end
+          IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
+             particle_start = start_index(nb)
+             particle_end   = end_index(nb)
+          ELSE
+             particle_start = 1
+             particle_end   = number_of_particles
+          ENDIF
+!
+!--         Loop from particle start to particle end
+            DO  n = particle_start, particle_end
+
 !
 !--          Transport of particles with inertia
@@ -3988,5 +4169,5 @@
              IF ( cloud_droplets )  THEN
 !
-!--             Terminal velocity is computed for vertical direction (Rogers et al., 
+!--             Terminal velocity is computed for vertical direction (Rogers et al.,
 !--             1993, J. Appl. Meteorol.)
                 diameter = particles(n)%radius * 2000.0_wp !diameter in mm
@@ -4054,6 +4235,29 @@
     particles(1:number_of_particles)%age_m = particles(1:number_of_particles)%age
 
-    DO  nb = 0, 7
-       DO  n = start_index(nb), end_index(nb)
+!
+!--    loop over subboxes. In case of simple interpolation scheme no subboxes
+!--    are introduced, as they are not required. Accordingly, this loops goes
+!--    from 1 to 1.
+!
+!-- Decide whether the particle loop runs over the subboxes or only over 1,
+!-- number_of_particles. This depends on the selected interpolation method.
+    IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
+       subbox_start = 0
+       subbox_end   = 7
+    ELSE
+       subbox_start = 1
+       subbox_end   = 1
+    ENDIF
+    DO  nb = subbox_start, subbox_end
+       IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
+          particle_start = start_index(nb)
+          particle_end   = end_index(nb)
+       ELSE
+          particle_start = 1
+          particle_end   = number_of_particles
+       ENDIF
+!
+!--    Loop from particle start to particle end
+       DO  n = particle_start, particle_end
 !
 !--       Increment the particle age and the total time that the particle
@@ -4063,5 +4267,5 @@
 
 !
-!--       Check whether there is still a particle that has not yet completed 
+!--       Check whether there is still a particle that has not yet completed
 !--       the total LES timestep
           IF ( ( dt_3d - particles(n)%dt_sum ) > 1E-8_wp )  THEN
@@ -4086,5 +4290,5 @@
  SUBROUTINE weil_stochastic_eq( v_sgs, fs_n, e_n, dedxi_n, dedt_n, diss_n,     &
                                 dt_n, rg_n, fac )
-                                
+
     REAL(wp) ::  a1      !< dummy argument
     REAL(wp) ::  dedt_n  !< time derivative of TKE at particle position 
@@ -4161,5 +4365,5 @@
     INTEGER(iwp), INTENT(IN) ::  j !< grid index of particle box along y
     INTEGER(iwp), INTENT(IN) ::  k !< grid index of particle box along z
-    
+
     INTEGER(iwp) ::  inc            !< dummy for sorting algorithmus
     INTEGER(iwp) ::  ir             !< dummy for sorting algorithmus
@@ -4184,5 +4388,5 @@
     INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions
     INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions
-    
+
     LOGICAL  ::  cross_wall_x    !< flag to check if particle reflection along x is necessary
     LOGICAL  ::  cross_wall_y    !< flag to check if particle reflection along y is necessary
@@ -4225,6 +4429,4 @@
        CASE ( 'bottom/top' )
 
-!     IF ( location_bc == 'bottom/top' )  THEN
-
 !
 !--    Apply boundary conditions to those particles that have crossed the top or
@@ -4270,5 +4472,5 @@
              ENDIF
           ENDIF
-          
+
           IF ( particles(n)%z < zw(0)  .AND.  particles(n)%particle_mask )  THEN
              IF ( ibc_par_b == 1 )  THEN
@@ -4290,5 +4492,4 @@
        ENDDO
 
-!     ELSEIF ( location_bc == 'walls' )  THEN
       CASE ( 'walls' )
 
@@ -4347,5 +4548,5 @@
           ELSE
              zwall = zw(k-1)
-          ENDIF     
+          ENDIF
 !
 !--       Initialize flags to check if particle reflection is necessary
@@ -4418,5 +4619,5 @@
                               MIN( prt_z - pos_z_old, -1E-30_wp ),             &
                               prt_z > pos_z_old )
-                     
+
           x_ind(t_index)   = i1
           y_ind(t_index)   = j1
@@ -4429,5 +4630,5 @@
              cross_wall_z = .TRUE.
           ENDIF
-          
+
           t_index_number = t_index - 1
 !
@@ -4484,5 +4685,5 @@
              t_old = 0.0_wp
              DO t_index = 1, t_index_number 
-!           
+!
 !--             Calculate intermediate particle position according to the 
 !--             timesteps a particle reaches any wall.
@@ -4583,5 +4784,5 @@
                     y_ind(t_index:t_index_number) = j2
                 ENDIF !particle reflection in y direction done
-                
+
 !
 !--             Check if a particle needs to be reflected at any xy-wall. If 
@@ -4623,5 +4824,5 @@
                     z_ind(t_index:t_index_number) = k2
                 ENDIF !particle reflection in z direction done                
-                
+
 !
 !--             Swap time
@@ -4655,6 +4856,16 @@
 
  END SUBROUTINE lpm_boundary_conds 
- 
- 
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculates change in droplet radius by condensation/evaporation, using
+!> either an analytic formula or by numerically integrating the radius growth
+!> equation including curvature and solution effects using Rosenbrocks method
+!> (see Numerical recipes in FORTRAN, 2nd edition, p. 731).
+!> The analytical formula and growth equation follow those given in
+!> Rogers and Yau (A short course in cloud physics, 3rd edition, p. 102/103).
+!------------------------------------------------------------------------------!
  SUBROUTINE lpm_droplet_condensation (i,j,k)
 
@@ -4812,5 +5023,5 @@
              DO
 
-                drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -    &
+                drdt = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - &
                                                           afactor / r_ros +    &
                                                           bfactor / r_ros**3   &
@@ -4818,26 +5029,26 @@
 
                 d2rdtdr = -ddenom * ventilation_effect(n) * (                  &
-                                                (e_a / e_s - 1.0) * r_ros**4 - &
-                                                afactor * r0 * r_ros**2 -      &
-                                                2.0 * afactor * r_ros**3 +     &
-                                                3.0 * bfactor * r0 +           &
-                                                4.0 * bfactor * r_ros          &
+                                            (e_a / e_s - 1.0_wp ) * r_ros**4 - &
+                                            afactor * r0 * r_ros**2 -          &
+                                            2.0_wp * afactor * r_ros**3 +      &
+                                            3.0_wp * bfactor * r0 +            &
+                                            4.0_wp * bfactor * r_ros           &
                                                             )                  &
                           / ( r_ros**4 * ( r_ros + r0 )**2 )
 
-                k1    = drdt / ( 1.0 - gamma * dt_ros * d2rdtdr )
+                k1    = drdt / ( 1.0_wp - gamma * dt_ros * d2rdtdr )
 
                 r_ros = MAX(r_ros_ini + k1 * dt_ros, particles(n)%aux1)
                 r_err = r_ros
 
-                drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0 -   &
-                                                           afactor / r_ros +   &
-                                                           bfactor / r_ros**3  &
+                drdt  = ddenom * ventilation_effect(n) * ( e_a / e_s - 1.0_wp - &
+                                                           afactor / r_ros +    &
+                                                           bfactor / r_ros**3   &
                                                          ) / ( r_ros + r0 )
 
                 k2 = ( drdt - dt_ros * 2.0 * gamma * d2rdtdr * k1 ) / &
-                     ( 1.0 - dt_ros * gamma * d2rdtdr )
-
-                r_ros = MAX(r_ros_ini + dt_ros * ( 1.5 * k1 + 0.5 * k2), particles(n)%aux1)
+                     ( 1.0_wp - dt_ros * gamma * d2rdtdr )
+
+                r_ros = MAX(r_ros_ini + dt_ros * ( 1.5_wp * k1 + 0.5_wp * k2), particles(n)%aux1)
    !
    !--          Check error of the solution, and reduce dt_ros if necessary.
@@ -4909,6 +5120,6 @@
           particles(n)%class = MAX( particles(n)%class, 1 )
        ENDIF
- !
- !--   Store new radius to particle features
+!
+!--    Store new radius to particle features
        particles(n)%radius = new_r(n)
 
@@ -5134,5 +5345,5 @@
           aero_mass(1:number_of_particles) = particles(1:number_of_particles)%weight_factor * &
                                              particles(1:number_of_particles)%aux1**3       * &
-                                             4.0 / 3.0 * pi * rho_s
+                                             4.0_wp / 3.0_wp * pi * rho_s
        ENDIF
 !
@@ -5144,5 +5355,5 @@
 !--          For collisions, the weighting factor of at least one super-droplet
 !--          needs to be larger or equal to one.
-             IF ( MIN( weight(n), weight(m) ) .LT. 1.0 )  CYCLE
+             IF ( MIN( weight(n), weight(m) ) .LT. 1.0_wp )  CYCLE
 !
 !--          Get mass of individual droplets (aerosols)
@@ -5175,5 +5386,5 @@
              ENDIF
 
-             IF ( collection_probability .GT. 0.0 )  THEN
+             IF ( collection_probability .GT. 0.0_wp )  THEN
 !
 !--             Super-droplet n collects droplets of super-droplet m
@@ -5363,5 +5574,5 @@
                                        ( hkernel(i,j), i = 1,radius_classes )
           ENDDO
-          
+
           DO  k = 1, dissipation_classes
              DO  i = 1, radius_classes
@@ -6154,5 +6365,5 @@
                      ql(k,j,i) < ql_crit )  CYCLE
                 particles => grid_particles(k,j,i)%particles(1:number_of_particles)
-!                
+!
 !--             Start splitting operations. Each particle is checked if it
 !--             fulfilled the splitting criterion's. In splitting mode 'const'   
@@ -6168,12 +6379,12 @@
                         particles(n)%radius >= radius_split  .AND.             & 
                         particles(n)%weight_factor >= weight_factor_split )    &
-                   THEN          
+                   THEN
 !
 !--                   Calculate the new number of particles.
-                      new_size = prt_count(k,j,i) + splitting_factor - 1                                    
+                      new_size = prt_count(k,j,i) + splitting_factor - 1
 !
 !--                   Cycle if maximum number of particles per grid box
 !--                   is greater than the allowed maximum number.
-                      IF ( new_size >= max_number_particles_per_gridbox )  CYCLE                      
+                      IF ( new_size >= max_number_particles_per_gridbox )  CYCLE
 !
 !--                   Reallocate particle array if necessary. 
@@ -6191,5 +6402,5 @@
                       DO  jpp = 1, splitting_factor-1
                          grid_particles(k,j,i)%particles(jpp+old_size) =       & 
-                            tmp_particle         
+                            tmp_particle
                       ENDDO  
                       new_particles_gb = new_particles_gb + splitting_factor - 1
@@ -6197,5 +6408,5 @@
 !--                   Save the new number of super droplets for every grid box.
                       prt_count(k,j,i) = prt_count(k,j,i) +                    &
-                                         splitting_factor - 1         
+                                         splitting_factor - 1
                    ENDIF
                 ENDDO
@@ -6216,11 +6427,11 @@
        m3           = 0.0_wp
        m3_total     = 0.0_wp
-       nr           = 0.0_wp   
+       nr           = 0.0_wp
        nrclgb       = 0.0_wp
-       nrclgb_total = 0.0_wp 
+       nrclgb_total = 0.0_wp
        nr_total     = 0.0_wp
        rm           = 0.0_wp
        rm_total     = 0.0_wp
-       
+
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -6231,6 +6442,6 @@
                 particles => grid_particles(k,j,i)%particles(1:number_of_particles)
                 nrclgb = nrclgb + 1.0_wp
-!                
-!--             Calculate moments of DSD.                
+!
+!--             Calculate moments of DSD.
                 DO  n = 1, number_of_particles
                    IF ( particles(n)%particle_mask  .AND.                      &
@@ -6242,13 +6453,13 @@
                                  particles(n)%weight_factor
                       IF ( isf == 1 )  THEN           
-                         diameter   = particles(n)%radius * 2.0_wp            
+                         diameter   = particles(n)%radius * 2.0_wp
                          lwc = lwc + factor_volume_to_mass *                   &
                                      particles(n)%radius**3 *                  & 
                                      particles(n)%weight_factor 
-                         m1  = m1  + particles(n)%weight_factor * diameter                                                
-                         m2  = m2  + particles(n)%weight_factor * diameter**2           
+                         m1  = m1  + particles(n)%weight_factor * diameter
+                         m2  = m2  + particles(n)%weight_factor * diameter**2
                          m3  = m3  + particles(n)%weight_factor * diameter**3
-                      ENDIF   
-                   ENDIF                        
+                      ENDIF
+                   ENDIF
                 ENDDO 
              ENDDO
@@ -6287,5 +6498,5 @@
                             ( nr_total * factor_volume_to_mass )               &
                           )**0.3333333_wp, 0.0_wp, nrclgb_total > 0.0_wp       &
-                       )                         
+                       )
 !
 !--    Check which function should be used to approximate the DSD.
@@ -6303,5 +6514,5 @@
                              0.0_wp, nrclgb_total > 0.0_wp                     &
                            )
-          zeta = m1_total * m3_total / m2_total**2                             
+          zeta = m1_total * m3_total / m2_total**2
           mu   = MAX( ( ( 1.0_wp - zeta ) * 2.0_wp + 1.0_wp ) /                &
                         ( zeta - 1.0_wp ), 0.0_wp                              &
@@ -6309,10 +6520,10 @@
 
           lambda = ( pirho_l * nr_total / lwc_total *                          &
-                     ( mu + 3.0_wp ) * ( mu + 2.0_wp ) * ( mu + 1.0_wp )       &                                          
+                     ( mu + 3.0_wp ) * ( mu + 2.0_wp ) * ( mu + 1.0_wp )       &
                    )**0.3333333_wp
           nr0 = nr_total / gamma( mu + 1.0_wp ) * lambda**( mu + 1.0_wp ) 
-          
+
           DO  n = 0, n_max-1
-             diameter  = r_bin_mid(n) * 2.0_wp           
+             diameter  = r_bin_mid(n) * 2.0_wp
              an_spl(n) = nr0 * diameter**mu * EXP( -lambda * diameter ) *      & 
                          ( r_bin(n+1) - r_bin(n) ) * 2.0_wp 
@@ -6338,27 +6549,27 @@
 !--    Criterion to avoid super droplets with a weighting factor < 1.0.
        an_spl = MAX(an_spl, 1.0_wp)
-  
+
        DO  i = nxl, nxr
           DO  j = nys, nyn
              DO  k = nzb+1, nzt
                 number_of_particles = prt_count(k,j,i)
-                IF ( number_of_particles <= 0  .OR.                            &  
+                IF ( number_of_particles <= 0  .OR.                            &
                      ql(k,j,i) < ql_crit )  CYCLE
                 particles => grid_particles(k,j,i)%particles(1:number_of_particles)
-                new_particles_gb = 0         
-!                
+                new_particles_gb = 0
+!
 !--             Start splitting operations. Each particle is checked if it
-!--             fulfilled the splitting criterion's. In splitting mode 'cl_av'   
-!--             a critical radius (radius_split) and a splitting function must 
+!--             fulfilled the splitting criterion's. In splitting mode 'cl_av'
+!--             a critical radius (radius_split) and a splitting function must
 !--             be prescribed (see particles_par). The critical weighting factor
 !--             is calculated while approximating a 'gamma', 'log' or 'exp'-
-!--             drop size distribution. In this mode the DSD is calculated as 
-!--             an average over all cloudy grid boxes. Super droplets which 
+!--             drop size distribution. In this mode the DSD is calculated as
+!--             an average over all cloudy grid boxes. Super droplets which
 !--             have a larger radius and larger weighting factor are split into
-!--             'splitting_factor' super droplets. In this case the splitting  
-!--             factor is calculated of weighting factor of the super droplet  
-!--             and the approximated number concentration for droplet of such 
-!--             a size. Due to the splitting, the weighting factor of the  
-!--             super droplet and all created clones is reduced by the factor  
+!--             'splitting_factor' super droplets. In this case the splitting
+!--             factor is calculated of weighting factor of the super droplet
+!--             and the approximated number concentration for droplet of such
+!--             a size. Due to the splitting, the weighting factor of the
+!--             super droplet and all created clones is reduced by the factor
 !--             of 'splitting_facor'.
                 DO  n = 1, number_of_particles
@@ -6379,9 +6590,9 @@
                          IF ( splitting_factor < 2 )  CYCLE
 !
-!--                      Calculate the new number of particles.                                                           
+!--                      Calculate the new number of particles.
                          new_size = prt_count(k,j,i) + splitting_factor - 1
 !
 !--                      Cycle if maximum number of particles per grid box
-!--                      is greater than the allowed maximum number.                         
+!--                      is greater than the allowed maximum number.
                          IF ( new_size >= max_number_particles_per_gridbox )   & 
                          CYCLE
@@ -6391,5 +6602,5 @@
                             CALL realloc_particles_array(i,j,k,new_size)
                          ENDIF
-                         old_size  = prt_count(k,j,i)                             
+                         old_size  = prt_count(k,j,i)
                          new_particles_gb = new_particles_gb +                 &
                                             splitting_factor - 1
@@ -6400,10 +6611,10 @@
                          tmp_particle = particles(n)
 !
-!--                      Create splitting_factor-1 new particles.                                                  
+!--                      Create splitting_factor-1 new particles.
                          DO  jpp = 1, splitting_factor-1
                             grid_particles(k,j,i)%particles(jpp+old_size) =    &
-                                                                    tmp_particle         
+                                                                    tmp_particle
                          ENDDO
-!   
+!
 !--                      Save the new number of super droplets. 
                          prt_count(k,j,i) = prt_count(k,j,i) +                 &
@@ -6412,5 +6623,5 @@
                    ENDDO
                 ENDDO 
-                
+
              ENDDO
           ENDDO
@@ -6422,7 +6633,7 @@
           DO  j = nys, nyn
              DO  k = nzb+1, nzt
-             
-!
-!--             Initialize summing variables.             
+
+!
+!--             Initialize summing variables.
                 lwc = 0.0_wp
                 m1  = 0.0_wp
@@ -6431,5 +6642,5 @@
                 nr  = 0.0_wp
                 rm  = 0.0_wp  
-                
+
                 new_particles_gb = 0
                 number_of_particles = prt_count(k,j,i)
@@ -6437,6 +6648,6 @@
                      ql(k,j,i) < ql_crit )  CYCLE
                 particles => grid_particles(k,j,i)%particles
-!                
-!--             Calculate moments of DSD.                
+!
+!--             Calculate moments of DSD.
                 DO  n = 1, number_of_particles
                    IF ( particles(n)%particle_mask  .AND.                      &
@@ -6447,22 +6658,22 @@
                                  particles(n)%radius**3  *                     &
                                  particles(n)%weight_factor
-                      IF ( isf == 1 )  THEN            
-                         diameter   = particles(n)%radius * 2.0_wp           
+                      IF ( isf == 1 )  THEN
+                         diameter   = particles(n)%radius * 2.0_wp
                          lwc = lwc + factor_volume_to_mass *                   &
-                                     particles(n)%radius**3 *                  &  
+                                     particles(n)%radius**3 *                  &
                                      particles(n)%weight_factor 
                          m1  = m1 + particles(n)%weight_factor * diameter
                          m2  = m2 + particles(n)%weight_factor * diameter**2
                          m3  = m3 + particles(n)%weight_factor * diameter**3
-                      ENDIF      
-                   ENDIF                                            
-                ENDDO  
-
-                IF ( nr <= 0.0  .OR.  rm <= 0.0_wp )  CYCLE 
-!
-!--             Calculate mean volume averaged radius.                
+                      ENDIF
+                   ENDIF
+                ENDDO
+
+                IF ( nr <= 0.0_wp  .OR.  rm <= 0.0_wp )  CYCLE
+!
+!--             Calculate mean volume averaged radius.
                 rm = ( rm / ( nr * factor_volume_to_mass ) )**0.3333333_wp
 !
-!--             Check which function should be used to approximate the DSD.              
+!--             Check which function should be used to approximate the DSD.
                 IF ( isf == 1 )  THEN
 !
@@ -6475,5 +6686,5 @@
                    lambda = ( pirho_l * nr / lwc *                             &
                               ( mu + 3.0_wp ) * ( mu + 2.0_wp ) *              &
-                              ( mu + 1.0_wp )                                  &                                  
+                              ( mu + 1.0_wp )                                  &
                             )**0.3333333_wp
                    nr0 =  ( nr / (gamma( mu + 1.0_wp ) ) ) *                   &
@@ -6481,5 +6692,5 @@
 
                    DO  n = 0, n_max-1
-                      diameter         = r_bin_mid(n) * 2.0_wp           
+                      diameter         = r_bin_mid(n) * 2.0_wp
                       an_spl(n) = nr0 * diameter**mu *                         &
                                   EXP( -lambda * diameter ) *                  & 
@@ -6489,5 +6700,5 @@
 !
 !--                Lognormal size distribution to calculate critical 
-!--                weight_factor (e.g. Levin, 1971, Bradley + Stow, 1974).   
+!--                weight_factor (e.g. Levin, 1971, Bradley + Stow, 1974).
                    DO  n = 0, n_max-1
                       an_spl(n) = nr / ( SQRT( 2.0_wp * pi ) *                 &
@@ -6495,6 +6706,6 @@
                                         ) *                                    &
                                   EXP( -( LOG( r_bin_mid(n) / rm )**2 ) /      &
-                                        ( 2.0_wp * LOG(sigma_log)**2 )         & 
-                                      ) *                                      & 
+                                        ( 2.0_wp * LOG(sigma_log)**2 )         &
+                                      ) *                                      &
                                   ( r_bin(n+1) - r_bin(n) )
                    ENDDO
@@ -6509,11 +6720,11 @@
                    ENDDO
                 ENDIF
-                
-!
-!--             Criterion to avoid super droplets with a weighting factor < 1.0.                                   
+
+!
+!--             Criterion to avoid super droplets with a weighting factor < 1.0.
                 an_spl = MAX(an_spl, 1.0_wp)
-!                
+!
 !--             Start splitting operations. Each particle is checked if it
-!--             fulfilled the splitting criterion's. In splitting mode 'gb_av'   
+!--             fulfilled the splitting criterion's. In splitting mode 'gb_av'
 !--             a critical radius (radius_split) and a splitting function must 
 !--             be prescribed (see particles_par). The critical weighting factor
@@ -6546,13 +6757,13 @@
 
 !
-!--                      Calculate the new number of particles.                                                                                         
+!--                      Calculate the new number of particles.
                          new_size = prt_count(k,j,i) + splitting_factor - 1
 !
 !--                      Cycle if maximum number of particles per grid box
-!--                      is greater than the allowed maximum number.                                                  
+!--                      is greater than the allowed maximum number.
                          IF ( new_size >= max_number_particles_per_gridbox )   &
                          CYCLE
 !
-!--                      Reallocate particle array if necessary.                          
+!--                      Reallocate particle array if necessary.
                          IF ( new_size > SIZE(particles) )  THEN 
                             CALL realloc_particles_array(i,j,k,new_size)
@@ -6575,5 +6786,5 @@
                                                splitting_factor - 1
                          new_particles_gb    = new_particles_gb +              &
-                                               splitting_factor - 1                                       
+                                               splitting_factor - 1
                       ENDIF
                    ENDDO  
@@ -6602,18 +6813,18 @@
  SUBROUTINE lpm_merging
 
-    INTEGER(iwp) ::  i         !<      
-    INTEGER(iwp) ::  j         !<       
-    INTEGER(iwp) ::  k         !<       
+    INTEGER(iwp) ::  i         !<
+    INTEGER(iwp) ::  j         !<
+    INTEGER(iwp) ::  k         !<
     INTEGER(iwp) ::  n         !<
     INTEGER(iwp) ::  merge_drp = 0     !< number of merged droplets
-    
-    
+
+
     REAL(wp) ::  ql_crit = 1.0E-5_wp  !< threshold lwc for cloudy grid cells 
                                       !< (e.g. Siebesma et al 2003, JAS, 60)
-    
+
     CALL cpu_log( log_point_s(81), 'lpm_merging', 'start' )
 
     merge_drp  = 0
-    
+
     IF ( weight_factor_merge == -1.0_wp )  THEN
        weight_factor_merge = 0.5_wp * initial_weighting_factor 
@@ -6623,5 +6834,5 @@
        DO  j = nys, nyn
           DO  k = nzb+1, nzt
-   
+
              number_of_particles = prt_count(k,j,i)
              IF ( number_of_particles <= 0  .OR.                               &
@@ -6635,5 +6846,5 @@
 !--          particle array. Tests showed that this simplified method can be 
 !--          used because it will only take place outside of cloudy grid 
-!--          boxes where ql <= 1.0E-5 kg/kg. Therefore, especially former cloned   
+!--          boxes where ql <= 1.0E-5 kg/kg. Therefore, especially former cloned
 !--          and subsequent evaporated super droplets will be merged.
              DO  n = 1, number_of_particles-1
@@ -6641,5 +6852,5 @@
                      particles(n+1)%particle_mask                  .AND.       &
                      particles(n)%radius        <= radius_merge    .AND.       &
-                     particles(n)%weight_factor <= weight_factor_merge )       &   
+                     particles(n)%weight_factor <= weight_factor_merge )       &
                 THEN
                    particles(n+1)%weight_factor  =                             &
@@ -6652,5 +6863,5 @@
                    deleted_particles          = deleted_particles + 1 
                    merge_drp                  = merge_drp + 1
-                
+
                 ENDIF
              ENDDO
@@ -6658,6 +6869,6 @@
        ENDDO
     ENDDO
-       
-       
+
+
     CALL cpu_log( log_point_s(81), 'lpm_merging', 'stop' )
 
@@ -7145,5 +7356,5 @@
 
        ALLOCATE(rvnp(MAX(1,trnp_count_recv)))
-!      
+!
 !--    This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit 
 !--    variables in structure particle_type (due to the calculation of par_size)
@@ -7164,5 +7375,5 @@
 
        ALLOCATE(rvsp(MAX(1,trsp_count_recv)))
-!      
+!
 !--    This MPI_SENDRECV should work even with odd mixture on 32 and 64 Bit 
 !--    variables in structure particle_type (due to the calculation of par_size)
@@ -7516,5 +7727,5 @@
              IF ( np_before_move <= 0 )  CYCLE
              particles_before_move => grid_particles(kp,jp,ip)%particles(1:np_before_move)
-             
+
              DO  n = 1, np_before_move
                 i = particles_before_move(n)%x * ddx
@@ -7578,7 +7789,7 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE lpm_check_cfl  
-       
+
     IMPLICIT NONE
-    
+
     INTEGER(iwp)  ::  i !< running index, x-direction
     INTEGER(iwp)  ::  j !< running index, y-direction
@@ -7759,7 +7970,8 @@
 ! -----------
 !> Routine for the whole processor
-!> Sort all particles into the 8 respective subgrid boxes
+!> Sort all particles into the 8 respective subgrid boxes and free space of
+!> particles which has been marked for deletion
 !------------------------------------------------------------------------------!
-   SUBROUTINE lpm_sort_in_subboxes
+   SUBROUTINE lpm_sort_and_delete
 
        INTEGER(iwp) ::  i  !<
@@ -7776,61 +7988,83 @@
        INTEGER(iwp), DIMENSION(0:7) ::  sort_count  !<
 
-       TYPE(particle_type), DIMENSION(:,:), ALLOCATABLE ::  sort_particles  !<
-
-       CALL cpu_log( log_point_s(51), 'lpm_sort_in_subboxes', 'start' )
-       DO  ip = nxl, nxr
-          DO  jp = nys, nyn
-             DO  kp = nzb+1, nzt
-                number_of_particles = prt_count(kp,jp,ip)
-                IF ( number_of_particles <= 0 )  CYCLE
-                particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
-                      
-                nn = 0
-                sort_count = 0
-                ALLOCATE( sort_particles(number_of_particles, 0:7) )
-                
-                DO  n = 1, number_of_particles
-                   sort_index = 0
-
-                   IF ( particles(n)%particle_mask )  THEN
-                      nn = nn + 1
-!
-!--                   Sorting particles with a binary scheme
-!--                   sort_index=111_2=7_10 -> particle at the left,south,bottom subgridbox
-!--                   sort_index=000_2=0_10 -> particle at the right,north,top subgridbox
-!--                   For this the center of the gridbox is calculated 
-                      i = (particles(n)%x + 0.5_wp * dx) * ddx
-                      j = (particles(n)%y + 0.5_wp * dy) * ddy
-                      
-                      IF ( i == ip )  sort_index = sort_index + 4 
-                      IF ( j == jp )  sort_index = sort_index + 2                      
-                      IF ( zu(kp) > particles(n)%z ) sort_index = sort_index + 1
-                      
-                      sort_count(sort_index) = sort_count(sort_index) + 1
-                      m = sort_count(sort_index)
-                      sort_particles(m,sort_index) = particles(n)
-                      sort_particles(m,sort_index)%block_nr = sort_index
-                   ENDIF
+       TYPE(particle_type), DIMENSION(:,:), ALLOCATABLE ::  sort_particles    !<
+
+       CALL cpu_log( log_point_s(51), 'lpm_sort_and_delete', 'start' )
+       IF ( TRIM(particle_interpolation)  == 'trilinear' ) THEN
+          DO  ip = nxl, nxr
+             DO  jp = nys, nyn
+                DO  kp = nzb+1, nzt
+                   number_of_particles = prt_count(kp,jp,ip)
+                   IF ( number_of_particles <= 0 )  CYCLE
+                   particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+                   nn = 0
+                   sort_count = 0
+                   ALLOCATE( sort_particles(number_of_particles, 0:7) )
+
+                   DO  n = 1, number_of_particles
+                      sort_index = 0
+
+                      IF ( particles(n)%particle_mask )  THEN
+                         nn = nn + 1
+!
+!--                      Sorting particles with a binary scheme
+!--                      sort_index=111_2=7_10 -> particle at the left,south,bottom subgridbox
+!--                      sort_index=000_2=0_10 -> particle at the right,north,top subgridbox
+!--                      For this the center of the gridbox is calculated
+                         i = (particles(n)%x + 0.5_wp * dx) * ddx
+                         j = (particles(n)%y + 0.5_wp * dy) * ddy
+
+                         IF ( i == ip )  sort_index = sort_index + 4
+                         IF ( j == jp )  sort_index = sort_index + 2
+                         IF ( zu(kp) > particles(n)%z ) sort_index = sort_index + 1
+
+                         sort_count(sort_index) = sort_count(sort_index) + 1
+                         m = sort_count(sort_index)
+                         sort_particles(m,sort_index) = particles(n)
+                         sort_particles(m,sort_index)%block_nr = sort_index
+                      ENDIF
+                   ENDDO
+
+                   nn = 0
+                   DO is = 0,7
+                      grid_particles(kp,jp,ip)%start_index(is) = nn + 1
+                      DO n = 1,sort_count(is)
+                         nn = nn + 1
+                         particles(nn) = sort_particles(n,is)
+                      ENDDO
+                      grid_particles(kp,jp,ip)%end_index(is) = nn
+                   ENDDO
+
+                   number_of_particles = nn
+                   prt_count(kp,jp,ip) = number_of_particles
+                   DEALLOCATE(sort_particles)
                 ENDDO
-
-                nn = 0
-                DO is = 0,7
-                   grid_particles(kp,jp,ip)%start_index(is) = nn + 1
-                   DO n = 1,sort_count(is)
-                      nn = nn + 1
-                      particles(nn) = sort_particles(n,is)
-                   ENDDO
-                   grid_particles(kp,jp,ip)%end_index(is) = nn
-                ENDDO
-
-                number_of_particles = nn    
-                prt_count(kp,jp,ip) = number_of_particles
-                DEALLOCATE(sort_particles)
              ENDDO
           ENDDO
-       ENDDO
-       CALL cpu_log( log_point_s(51), 'lpm_sort_in_subboxes', 'stop' )
-
-    END SUBROUTINE lpm_sort_in_subboxes
+
+!--    In case of the simple interpolation method the particles must not
+!--    be sorted in subboxes but particles marked for deletion must be
+!--    deleted and number of particles must be recalculated
+       ELSE
+
+          DO  ip = nxl, nxr
+             DO  jp = nys, nyn
+                DO  kp = nzb+1, nzt
+
+                   number_of_particles = prt_count(kp,jp,ip)
+                   IF ( number_of_particles <= 0 )  CYCLE
+                   particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
+!
+!--                repack particles array, i.e. delete particles and recalculate
+!--                number of particles
+                   CALL lpm_pack
+                   prt_count(kp,jp,ip) = number_of_particles
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+       CALL cpu_log( log_point_s(51), 'lpm_sort_and_delete', 'stop' )
+
+    END SUBROUTINE lpm_sort_and_delete
 
  
@@ -7841,5 +8075,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE lpm_pack
-    
+
        INTEGER(iwp) ::  n       !<
        INTEGER(iwp) ::  nn      !<
@@ -7874,6 +8108,6 @@
 
     END SUBROUTINE lpm_pack 
-                
-                
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -7904,5 +8138,4 @@
                 number_of_particles = prt_count(k,j,i)
                 IF ( number_of_particles <= 0 )  CYCLE
-
                 particles => grid_particles(k,j,i)%particles(1:number_of_particles)
 
