Index: palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4122)
+++ palm/trunk/SOURCE/lagrangian_particle_model_mod.f90	(revision 4143)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Rename variable and change select case to if statement
+! 
+! 4122 2019-07-26 13:11:56Z schwenkel
 ! Implement reset method as bottom boundary condition
 ! 
@@ -278,5 +281,5 @@
     CHARACTER(LEN=5)  ::  splitting_mode     = 'const'            !< splitting mode
 
-    CHARACTER(LEN=25) ::  particle_interpolation = 'trilinear'    !< interpolation method for calculatin the particle
+    CHARACTER(LEN=25) ::  particle_advection_interpolation = 'trilinear' !< interpolation method for calculatin the particle
 
     INTEGER(iwp) ::  deleted_particles = 0                        !< number of deleted particles per time step    
@@ -315,4 +318,7 @@
     LOGICAL ::  use_kernel_tables = .FALSE.               !< parameter, which turns on the use of precalculated collision kernels
     LOGICAL ::  write_particle_statistics = .FALSE.       !< namelist parameter (see documentation)
+    LOGICAL ::  interpolation_simple_predictor = .FALSE.  !< flag for simple particle advection interpolation with predictor step
+    LOGICAL ::  interpolation_simple_corrector = .FALSE.  !< flag for simple particle advection interpolation with corrector step
+    LOGICAL ::  interpolation_trilinear = .FALSE.         !< flag for trilinear particle advection interpolation
 
     LOGICAL, DIMENSION(max_number_of_particle_groups) ::   vertical_particle_advection = .TRUE. !< Switch for vertical particle transport
@@ -555,5 +561,5 @@
        particles_per_point, &
        particle_advection_start, &
-       particle_interpolation, &
+       particle_advection_interpolation, &
        particle_maximum_age, &
        pdx, &
@@ -615,5 +621,5 @@
        particles_per_point, &
        particle_advection_start, &
-       particle_interpolation, &
+       particle_advection_interpolation, &
        particle_maximum_age, &
        pdx, &
@@ -862,5 +868,5 @@
 !-- 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.              &
+    IF ( .NOT. TRIM(particle_advection_interpolation) == 'trilinear'  .AND.              &
        use_sgs_for_particles .AND.  .NOT. cloud_droplets )  THEN
           message_string = 'subrgrid scale velocities in combination with ' // &
@@ -1012,6 +1018,6 @@
        de_dx = 0.0_wp
        de_dy = 0.0_wp
-       de_dz = 0.0_wp             
-                 
+       de_dz = 0.0_wp
+
        sgs_wf_part = 1.0_wp / 3.0_wp
     ENDIF
@@ -1063,4 +1069,20 @@
        ENDDO
 
+    ENDIF
+
+!
+!-- Check which particle interpolation method should be used
+    IF ( TRIM(particle_advection_interpolation)  ==  'trilinear' )  THEN
+       interpolation_simple_corrector = .FALSE.
+       interpolation_simple_predictor = .FALSE.
+       interpolation_trilinear        = .TRUE.
+    ELSEIF ( TRIM(particle_advection_interpolation)  ==  'simple_corrector' )  THEN
+       interpolation_simple_corrector = .TRUE.
+       interpolation_simple_predictor = .FALSE.
+       interpolation_trilinear        = .FALSE.
+    ELSEIF ( TRIM(particle_advection_interpolation)  ==  'simple_predictor' )  THEN
+       interpolation_simple_corrector = .FALSE.
+       interpolation_simple_predictor = .TRUE.
+       interpolation_trilinear        = .FALSE.
     ENDIF
 
@@ -1309,6 +1331,6 @@
 !
 !--    Calculate initial_weighting_factor diagnostically
-       IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp ) THEN
-          initial_weighting_factor =  number_concentration  *                        &
+       IF ( number_concentration /= -1.0_wp .AND. number_concentration > 0.0_wp )  THEN
+          initial_weighting_factor =  number_concentration  *                         &
                                       pdx(1) * pdy(1) * pdz(1)
        END IF
@@ -1322,9 +1344,9 @@
                 DO WHILE ( pos_y <= psn(i) )
                    IF ( pos_y >= nys * dy  .AND.                  &
-                        pos_y <  ( nyn + 1 ) * dy  ) THEN
+                        pos_y <  ( nyn + 1 ) * dy  )  THEN
                       pos_x = psl(i)
                xloop: DO WHILE ( pos_x <= psr(i) )
                          IF ( pos_x >= nxl * dx  .AND.            &
-                              pos_x <  ( nxr + 1) * dx ) THEN
+                              pos_x <  ( nxr + 1) * dx )  THEN
                             DO  j = 1, particles_per_point
                                n = n + 1
@@ -1366,5 +1388,5 @@
 !--                            In case of stretching the actual k index is found iteratively
                                IF ( dz_stretch_level .NE. -9999999.9_wp  .OR.           &
-                                    dz_stretch_level_start(1) .NE. -9999999.9_wp ) THEN
+                                    dz_stretch_level_start(1) .NE. -9999999.9_wp )  THEN
                                   kp = MINLOC( ABS( tmp_particle%z - zu ), DIM = 1 ) - 1
                                ELSE
@@ -1626,13 +1648,13 @@
 !
 !-- Set constants for different aerosol species
-    IF ( TRIM(aero_species) .EQ. 'nacl' ) THEN
+    IF ( TRIM(aero_species) .EQ. 'nacl' )  THEN
        molecular_weight_of_solute = 0.05844_wp 
        rho_s                      = 2165.0_wp
        vanthoff                   = 2.0_wp
-    ELSEIF ( TRIM(aero_species) .EQ. 'c3h4o4' ) THEN
+    ELSEIF ( TRIM(aero_species) .EQ. 'c3h4o4' )  THEN
        molecular_weight_of_solute = 0.10406_wp 
        rho_s                      = 1600.0_wp
        vanthoff                   = 1.37_wp
-    ELSEIF ( TRIM(aero_species) .EQ. 'nh4o3' ) THEN
+    ELSEIF ( TRIM(aero_species) .EQ. 'nh4o3' )  THEN
        molecular_weight_of_solute = 0.08004_wp 
        rho_s                      = 1720.0_wp
@@ -2153,5 +2175,5 @@
 
                       IF ( .NOT. first_loop_stride  .AND.  &
-                           grid_particles(k,j,i)%time_loop_done ) CYCLE
+                           grid_particles(k,j,i)%time_loop_done )  CYCLE
 
                       particles => grid_particles(k,j,i)%particles(1:number_of_particles)
@@ -2187,5 +2209,5 @@
 !
 !--                   Particle advection
-                      CALL lpm_advec(TRIM(particle_interpolation),i,j,k)
+                      CALL lpm_advec(i,j,k)
 !
 !--                   Particle reflection from walls. Only applied if the particles 
@@ -2243,8 +2265,8 @@
 !--          Apply splitting and merging algorithm
              IF ( cloud_droplets )  THEN
-                IF ( splitting ) THEN
+                IF ( splitting )  THEN
                    CALL lpm_splitting
                 ENDIF
-                IF ( merging ) THEN
+                IF ( merging )  THEN
                    CALL lpm_merging
                 ENDIF
@@ -2392,5 +2414,5 @@
     CALL close_file( 80 )
 
-    IF ( number_of_particles > 0 ) THEN
+    IF ( number_of_particles > 0 )  THEN
         WRITE(9,*) 'number_of_particles ', number_of_particles,                &
                     current_timestep_number + 1, simulated_time + dt_3d
@@ -3131,4 +3153,13 @@
     WRITE ( 14 )  curvature_solution_effects
 
+    CALL wrd_write_string( 'interpolation_simple_corrector' )
+    WRITE ( 14 )  interpolation_simple_corrector
+
+    CALL wrd_write_string( 'interpolation_simple_predictor' )
+    WRITE ( 14 )  interpolation_simple_predictor
+
+    CALL wrd_write_string( 'interpolation_trilinear' )
+    WRITE ( 14 )  interpolation_trilinear
+
  END SUBROUTINE lpm_wrd_global
  
@@ -3152,5 +3183,14 @@
        CASE ( 'curvature_solution_effects' )
           READ ( 13 )  curvature_solution_effects
-          
+
+       CASE ( 'interpolation_simple_corrector' )
+          READ ( 13 )  interpolation_simple_corrector
+
+       CASE ( 'interpolation_simple_predictor' )
+          READ ( 13 )  interpolation_simple_predictor
+
+       CASE ( 'interpolation_trilinear' )
+          READ ( 13 )  interpolation_trilinear
+
 !          CASE ( 'global_paramter' )
 !             READ ( 13 )  global_parameter
@@ -3176,7 +3216,6 @@
 !> this submodule should be excluded as an own file. 
 !------------------------------------------------------------------------------!
- SUBROUTINE lpm_advec (interpolation_method,ip,jp,kp)
-
-    CHARACTER (LEN=*), INTENT(IN) ::  interpolation_method !<
+ SUBROUTINE lpm_advec (ip,jp,kp)
+
     LOGICAL ::  subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall
 
@@ -3307,340 +3346,328 @@
     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) ) /                   &
+!
+!-- 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).
+    IF ( interpolation_simple_corrector )  THEN
+!
+!--    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.
+    ELSEIF ( interpolation_simple_predictor )  THEN
+!
+!--    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
-
+          w_int(n) = w(kkw,jp,ip) * ( 1.0_wp - gamma ) + w(kkw+1,jp,ip) * gamma
+       ENDDO
+!
+!-- The trilinear interpolation.
+    ELSEIF ( interpolation_trilinear )  THEN
+
+       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
           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
+!
+!--       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
-
-
-!
-!--    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
+!
+!--       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
-!
-!--    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
-             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
+       ENDDO
+    ENDIF
 
 !-- Interpolate and calculate quantities needed for calculating the SGS
@@ -3653,5 +3680,5 @@
 !
 !--       In case of topography check if subbox is adjacent to a wall
-          IF ( .NOT. topography == 'flat' ) THEN
+          IF ( .NOT. topography == 'flat' )  THEN
              i = ip + MERGE( -1_iwp , 1_iwp, BTEST( nb, 2 ) )
              j = jp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 1 ) )
@@ -3664,5 +3691,5 @@
              ENDIF
           ENDIF
-          IF ( subbox_at_wall ) THEN
+          IF ( subbox_at_wall )  THEN
              e_int(start_index(nb):end_index(nb))     = e(kp,jp,ip) 
              diss_int(start_index(nb):end_index(nb))  = diss(kp,jp,ip)
@@ -3960,5 +3987,5 @@
                 (ABS( u_int(n) + rvar1_temp(n) ) > (dx/dt_particle(n))  .OR.   &
                  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
+                 ABS( w_int(n) + rvar3_temp(n) ) > (dz_temp/dt_particle(n))))  THEN
 
                 dt_particle(n) = 0.9_wp * MIN(                                 &
@@ -4027,5 +4054,5 @@
 !--    number_of_particles. @todo find a more generic way to write this loop or
 !--    delete trilinear interpolation
-       IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
+       IF ( interpolation_trilinear )  THEN
           subbox_start = 0
           subbox_end   = 7
@@ -4039,5 +4066,5 @@
 !--    from 1 to 1.
        DO  nb = subbox_start, subbox_end
-          IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
+          IF ( interpolation_trilinear )  THEN
              particle_start = start_index(nb)
              particle_end   = end_index(nb)
@@ -4141,5 +4168,5 @@
 !--    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
+       IF ( interpolation_trilinear )  THEN
           subbox_start = 0
           subbox_end   = 7
@@ -4152,5 +4179,5 @@
 !--    from 1 to 1.
        DO  nb = subbox_start, subbox_end
-          IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
+          IF ( interpolation_trilinear )  THEN
              particle_start = start_index(nb)
              particle_end   = end_index(nb)
@@ -4245,5 +4272,5 @@
 !-- 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
+    IF ( interpolation_trilinear )  THEN
        subbox_start = 0
        subbox_end   = 7
@@ -4253,5 +4280,5 @@
     ENDIF
     DO  nb = subbox_start, subbox_end
-       IF ( TRIM(particle_interpolation)  ==  'trilinear' )  THEN
+       IF ( interpolation_trilinear )  THEN
           particle_start = start_index(nb)
           particle_end   = end_index(nb)
@@ -4505,5 +4532,5 @@
                 particles(n)%z       = reset_top *  ( 1.0  + ( ran_val / 10.0_wp) )
                 particles(n)%speed_z = 0.0_wp
-                IF ( curvature_solution_effects ) THEN
+                IF ( curvature_solution_effects )  THEN
                    particles(n)%radius = particles(n)%aux1
                 ELSE
@@ -4566,5 +4593,5 @@
           ENDIF
 
-          IF ( prt_z > pos_z_old ) THEN
+          IF ( prt_z > pos_z_old )  THEN
              zwall = zw(k)
           ELSE
@@ -4648,5 +4675,5 @@
           reach_y(t_index) = .FALSE.
           reach_z(t_index) = .TRUE.
-          IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp) THEN
+          IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp)  THEN
              t_index      = t_index + 1
              cross_wall_z = .TRUE.
@@ -5001,5 +5028,5 @@
     ENDDO
 
-    IF( .NOT. curvature_solution_effects ) then
+    IF( .NOT. curvature_solution_effects )  THEN
 !
 !--    Use analytic model for diffusional growth including gas-kinetic
@@ -6136,5 +6163,5 @@
 
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ira !<
-    
+
     LOGICAL, SAVE ::  first = .TRUE. !<
 
@@ -6149,5 +6176,5 @@
     REAL(wp), DIMENSION(1:11), SAVE ::  rat           !<
     REAL(wp), DIMENSION(1:7), SAVE  ::  r0            !<
-    
+
     REAL(wp), DIMENSION(1:7,1:11), SAVE ::  ecoll_100 !<
     REAL(wp), DIMENSION(1:7,1:11), SAVE ::  ecoll_400 !<
@@ -7010,5 +7037,5 @@
           DO  kp = nzb+1, nzt
              number_of_particles = prt_count(kp,jp,ip)
-             IF ( number_of_particles <= 0 ) CYCLE
+             IF ( number_of_particles <= 0 )  CYCLE
              particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
              DO  n = 1, number_of_particles
@@ -7589,5 +7616,5 @@
 !--    In case of stretching the actual k index must be found
        IF ( dz_stretch_level .NE. -9999999.9_wp  .OR.         &
-            dz_stretch_level_start(1) .NE. -9999999.9_wp ) THEN
+            dz_stretch_level_start(1) .NE. -9999999.9_wp )  THEN
           kp = MINLOC( ABS( particle_array(n)%z - zu ), DIM = 1 ) - 1
        ELSE
@@ -7992,6 +8019,7 @@
 ! -----------
 !> Routine for the whole processor
-!> Sort all particles into the 8 respective subgrid boxes and free space of
-!> particles which has been marked for deletion
+!> Sort all particles into the 8 respective subgrid boxes (in case of trilinear
+!> interpolation method) and free space of particles which has been marked for
+!> deletion.
 !------------------------------------------------------------------------------!
    SUBROUTINE lpm_sort_and_delete
@@ -8013,5 +8041,5 @@
 
        CALL cpu_log( log_point_s(51), 'lpm_sort_and_delete', 'start' )
-       IF ( TRIM(particle_interpolation)  == 'trilinear' ) THEN
+       IF ( interpolation_trilinear )  THEN
           DO  ip = nxl, nxr
              DO  jp = nys, nyn
@@ -8047,5 +8075,7 @@
                       ENDIF
                    ENDDO
-
+!
+!--                Delete and resort particles by overwritting and set
+!--                the number_of_particles to the actual value.
                    nn = 0
                    DO is = 0,7
@@ -8066,6 +8096,7 @@
 
 !--    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
+!--    be sorted in subboxes. Particles marked for deletion however, must be
+!--    deleted and number of particles must be recalculated as it is also
+!--    done for the trilinear particle advection interpolation method.
        ELSE
 
@@ -8078,5 +8109,5 @@
                    particles => grid_particles(kp,jp,ip)%particles(1:number_of_particles)
 !
-!--                repack particles array, i.e. delete particles and recalculate
+!--                Repack particles array, i.e. delete particles and recalculate
 !--                number of particles
                    CALL lpm_pack
