Index: palm/trunk/SOURCE/lpm_droplet_collision.f90
===================================================================
--- palm/trunk/SOURCE/lpm_droplet_collision.f90	(revision 1323)
+++ palm/trunk/SOURCE/lpm_droplet_collision.f90	(revision 1359)
@@ -1,3 +1,3 @@
- SUBROUTINE lpm_droplet_collision
+ SUBROUTINE lpm_droplet_collision (i,j,k)
 
 !--------------------------------------------------------------------------------!
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! ------------------
-!
+! New particle structure integrated. 
+! Kind definition added to all floating point numbers.
 !
 ! Former revisions:
@@ -75,4 +76,5 @@
 !------------------------------------------------------------------------------!
 
+
     USE arrays_3d,                                                             &
         ONLY:  diss, ql, ql_v, ql_vp, u, v, w, zu, zw
@@ -103,6 +105,8 @@
     USE particle_attributes,                                                   &
         ONLY:  deleted_particles, dissipation_classes, hall_kernel,            &
-               palm_kernel, particles, particle_mask, particle_type,           &
-               prt_count, prt_start_index, use_kernel_tables, wang_kernel
+               palm_kernel, particles, particle_type,                          &
+               prt_count, use_kernel_tables, wang_kernel
+
+    USE pegrid
 
     IMPLICIT NONE
@@ -124,5 +128,13 @@
     INTEGER(iwp) ::  rclass_s !:
 
+    INTEGER(iwp), DIMENSION(prt_count(k,j,i)) ::  rclass_v !:
+
+    LOGICAL, SAVE ::  first_flag = .TRUE. !:
+
+    TYPE(particle_type) :: tmp_particle   !:
+
     REAL(wp) ::  aa       !:
+    REAL(wp) ::  auxn     !: temporary variables
+    REAL(wp) ::  auxs     !: temporary variables
     REAL(wp) ::  bb       !:
     REAL(wp) ::  cc       !:
@@ -158,535 +170,522 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  weight !:
 
-
-    TYPE(particle_type) ::  tmp_particle           !:
-
-
+    REAL, DIMENSION(prt_count(k,j,i))    :: ck
+    REAL, DIMENSION(prt_count(k,j,i))    :: r3v
+    REAL, DIMENSION(prt_count(k,j,i))    :: sum1v
+    REAL, DIMENSION(prt_count(k,j,i))    :: sum2v
 
     CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'start' )
 
-    DO  i = nxl, nxr
-       DO  j = nys, nyn
-          DO  k = nzb+1, nzt
-!
-!--          Collision requires at least two particles in the box
-             IF ( prt_count(k,j,i) > 1 )  THEN
-!
-!--             First, sort particles within the gridbox by their size,
-!--             using Shell's method (see Numerical Recipes)
-!--             NOTE: In case of using particle tails, the re-sorting of
-!--             ----  tails would have to be included here!
-                psi = prt_start_index(k,j,i) - 1
-                inc = 1
-                DO WHILE ( inc <= prt_count(k,j,i) )
-                   inc = 3 * inc + 1
+!
+!-- Collision requires at least two particles in the box
+    IF ( prt_count(k,j,i) > 1 )  THEN
+!
+!--    First, sort particles within the gridbox by their size,
+!--    using Shell's method (see Numerical Recipes)
+!--    NOTE: In case of using particle tails, the re-sorting of
+!--    ----  tails would have to be included here!
+       IF ( .NOT. ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
+       use_kernel_tables ) )  THEN
+          psi = 0
+          inc = 1
+          DO WHILE ( inc <= prt_count(k,j,i) )
+             inc = 3 * inc + 1
+          ENDDO
+
+          DO WHILE ( inc > 1 )
+             inc = inc / 3
+             DO  is = inc+1, prt_count(k,j,i)
+                tmp_particle = particles(psi+is)
+                js = is
+                DO WHILE ( particles(psi+js-inc)%radius >             &
+                tmp_particle%radius )
+                   particles(psi+js) = particles(psi+js-inc)
+                   js = js - inc
+                   IF ( js <= inc )  EXIT
                 ENDDO
-
-                DO WHILE ( inc > 1 )
-                   inc = inc / 3
-                   DO  is = inc+1, prt_count(k,j,i)
-                      tmp_particle = particles(psi+is)
-                      js = is
-                      DO WHILE ( particles(psi+js-inc)%radius >             &
-                                 tmp_particle%radius )
-                         particles(psi+js) = particles(psi+js-inc)
-                         js = js - inc
-                         IF ( js <= inc )  EXIT
-                      ENDDO
-                      particles(psi+js) = tmp_particle
-                   ENDDO
+                particles(psi+js) = tmp_particle
+             ENDDO
+          ENDDO
+       ENDIF
+
+       psi = 1
+       pse = prt_count(k,j,i)
+
+!
+!--    Now apply the different kernels
+       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
+            use_kernel_tables )  THEN
+!
+!--       Fast method with pre-calculated efficiencies for
+!--       discrete radius- and dissipation-classes.
+!--
+!--       Determine dissipation class index of this gridbox
+          IF ( wang_kernel )  THEN
+             eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
+                           dissipation_classes ) + 1
+             epsilon = diss(k,j,i)
+          ELSE
+             epsilon = 0.0_wp
+          ENDIF
+          IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
+             eclass = 0   ! Hall kernel is used
+          ELSE
+             eclass = MIN( dissipation_classes, eclass )
+          ENDIF
+
+!
+!--       Droplet collision are calculated using collision-coalescence
+!--       formulation proposed by Wang (see PALM documentation)
+!--       Since new radii after collision are defined by radii of all
+!--       droplets before collision, temporary fields for new radii and
+!--       weighting factors are needed
+          ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
+
+          rad    = 0.0_wp
+          weight = 0.0_wp
+
+          sum1v(1:prt_count(k,j,i)) = 0.0_wp
+          sum2v(1:prt_count(k,j,i)) = 0.0_wp
+
+          DO  n = 1, prt_count(k,j,i)
+
+             rclass_l = particles(n)%class
+!
+!--          Mass added due to collisions with smaller droplets
+             DO  is = n+1, prt_count(k,j,i)
+                rclass_s = particles(is)%class
+                auxs = ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor
+                auxn = ckernel(rclass_s,rclass_l,eclass) * particles(n)%weight_factor
+                IF ( particles(is)%radius <  particles(n)%radius )  THEN
+                   sum1v(n) =  sum1v(n)  + particles(is)%radius**3 * auxs
+                   sum2v(is) = sum2v(is) + auxn
+                ELSE
+                   sum2v(n)  = sum2v(n)  + auxs
+                   sum1v(is) = sum1v(is) + particles(n)%radius**3 * auxn
+                ENDIF
+             ENDDO
+          ENDDO
+          rclass_v = particles(1:prt_count(k,j,i))%class
+          DO  n = 1, prt_count(k,j,i)
+            ck(n)       = ckernel(rclass_v(n),rclass_v(n),eclass)
+          ENDDO
+          r3v      = particles(1:prt_count(k,j,i))%radius**3
+          DO  n = 1, prt_count(k,j,i)
+             sum3 = 0.0_wp
+             ddV = ddx * ddy / dz
+!
+!--          Change of the current weighting factor
+             sum3 = 1 - dt_3d * ddV *                                 &
+                        ck(n) *                                       &
+                        ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
+                    dt_3d * ddV * sum2v(n)
+             weight(n) = particles(n)%weight_factor * sum3
+!
+!--          Change of the current droplet radius
+             rad(n) = ( (r3v(n) + dt_3d * ddV * (sum1v(n) - sum2v(n) * r3v(n)) )/&
+                              sum3 )**0.33333333333333_wp
+
+             ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n) &
+                                         * rad(n)**3
+
+          ENDDO
+          IF ( ANY(weight < 0.0_wp) )  THEN
+                WRITE( message_string, * ) 'negative weighting'
+                CALL message( 'lpm_droplet_collision', 'PA0028',      &
+                               2, 2, -1, 6, 1 )
+          ENDIF
+
+          particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
+          particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
+
+          DEALLOCATE(rad, weight)
+
+       ELSEIF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
+                .NOT. use_kernel_tables )  THEN
+!
+!--       Collision efficiencies are calculated for every new
+!--       grid box. First, allocate memory for kernel table.
+!--       Third dimension is 1, because table is re-calculated for
+!--       every new dissipation value.
+          ALLOCATE( ckernel(1:prt_count(k,j,i),1:prt_count(k,j,i),1:1) )
+!
+!--       Now calculate collision efficiencies for this box
+          CALL recalculate_kernel( i, j, k )
+
+!
+!--       Droplet collision are calculated using collision-coalescence
+!--       formulation proposed by Wang (see PALM documentation)
+!--       Since new radii after collision are defined by radii of all
+!--       droplets before collision, temporary fields for new radii and
+!--       weighting factors are needed
+          ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
+
+          rad = 0.0_wp
+          weight = 0.0_wp
+
+          DO  n = psi, pse, 1
+
+             sum1 = 0.0_wp
+             sum2 = 0.0_wp
+             sum3 = 0.0_wp
+!
+!--          Mass added due to collisions with smaller droplets
+             DO  is = psi, n-1
+                sum1 = sum1 + ( particles(is)%radius**3 *            &
+                                ckernel(n,is,1) *  &
+                                particles(is)%weight_factor )
+             ENDDO
+!
+!--          Rate of collisions with larger droplets
+             DO  is = n+1, pse
+                sum2 = sum2 + ( ckernel(n,is,1) *  &
+                                particles(is)%weight_factor )
+             ENDDO
+
+             r3 = particles(n)%radius**3
+             ddV = ddx * ddy / dz 
+             is = 1
+!
+!--                   Change of the current weighting factor
+             sum3 = 1 - dt_3d * ddV *                                 &
+                        ckernel(n,n,1) *           &
+                        ( particles(n)%weight_factor - 1 ) * 0.5_wp - &
+                    dt_3d * ddV * sum2
+             weight(n-is+1) = particles(n)%weight_factor * sum3
+!
+!--                   Change of the current droplet radius
+             rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
+                              sum3 )**0.33333333333333_wp
+
+             IF ( weight(n-is+1) < 0.0_wp )  THEN
+                WRITE( message_string, * ) 'negative weighting',      &
+                                           'factor: ', weight(n-is+1)
+                CALL message( 'lpm_droplet_collision', 'PA0037',      &
+                               2, 2, -1, 6, 1 )
+             ENDIF
+
+             ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
+                                         * rad(n-is+1)**3
+
+          ENDDO
+
+          particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
+          particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
+
+          DEALLOCATE( rad, weight, ckernel )
+
+       ELSEIF ( palm_kernel )  THEN
+!
+!--       PALM collision kernel
+!
+!--       Calculate the mean radius of all those particles which
+!--       are of smaller size than the current particle and
+!--       use this radius for calculating the collision efficiency
+          DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
+
+             sl_r3 = 0.0_wp
+             sl_r4 = 0.0_wp
+
+             DO is = n-1, psi, -1
+                IF ( particles(is)%radius < particles(n)%radius )  &
+                THEN
+                   sl_r3 = sl_r3 + particles(is)%weight_factor     &
+                                   * particles(is)%radius**3
+                   sl_r4 = sl_r4 + particles(is)%weight_factor     &
+                                   * particles(is)%radius**4
+                ENDIF
+             ENDDO
+
+             IF ( ( sl_r3 ) > 0.0_wp )  THEN
+                mean_r = ( sl_r4 ) / ( sl_r3 )
+
+                CALL collision_efficiency_rogers( mean_r,             &
+                                           particles(n)%radius,       &
+                                           effective_coll_efficiency )
+
+             ELSE
+                effective_coll_efficiency = 0.0_wp
+             ENDIF
+
+             IF ( effective_coll_efficiency > 1.0_wp  .OR.            &
+                  effective_coll_efficiency < 0.0_wp )                &
+             THEN
+                WRITE( message_string, * )  'collision_efficien' , &
+                          'cy out of range:' ,effective_coll_efficiency
+                CALL message( 'lpm_droplet_collision', 'PA0145', 2, &
+                              2, -1, 6, 1 )
+             ENDIF
+
+!
+!--          Interpolation of liquid water content
+             ii = particles(n)%x * ddx
+             jj = particles(n)%y * ddy
+             kk = ( particles(n)%z + 0.5_wp * dz ) / dz
+
+             x  = particles(n)%x - ii * dx
+             y  = particles(n)%y - jj * 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
+
+             ql_int_l = ( (gg-aa) * ql(kk,jj,ii)   + (gg-bb) *     &
+                                                  ql(kk,jj,ii+1)   &
+                        + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) *   &
+                                                  ql(kk,jj+1,ii+1) &
+                        ) / ( 3.0_wp * gg )
+
+             ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii)   + (gg-bb) *   &
+                                                ql(kk+1,jj,ii+1)   &
+                        + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) *   &
+                                                ql(kk+1,jj+1,ii+1) &
+                        ) / ( 3.0_wp * gg )
+
+             ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
+                                 ( ql_int_u - ql_int_l )
+
+!
+!--          Interpolate u velocity-component
+             ii = ( particles(n)%x + 0.5_wp * dx ) * ddx
+             jj =   particles(n)%y * ddy
+             kk = ( particles(n)%z + 0.5_wp * dz ) / dz ! only if equidistant
+
+             IF ( ( particles(n)%z - zu(kk) ) > ( 0.5_wp * dz ) ) kk = kk+1
+
+             x  = particles(n)%x + ( 0.5_wp - ii ) * dx
+             y  = particles(n)%y - jj * 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(kk,jj,ii)   + (gg-bb) *       &
+                                                  u(kk,jj,ii+1)    &
+                       + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) *       &
+                                                  u(kk,jj+1,ii+1)  &
+                       ) / ( 3.0_wp * gg ) - u_gtrans
+             IF ( kk+1 == nzt+1 )  THEN
+                u_int = u_int_l
+             ELSE
+                u_int_u = ( (gg-aa) * u(kk+1,jj,ii)   + (gg-bb) *  &
+                                                 u(kk+1,jj,ii+1)   &
+                          + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) *  &
+                                                 u(kk+1,jj+1,ii+1) &
+                          ) / ( 3.0_wp * gg ) - u_gtrans
+                u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
+                                  * ( u_int_u - u_int_l )
+             ENDIF
+
+!
+!--          Same procedure for interpolation of the v velocity-component 
+!--          (adopt index k from u velocity-component)
+             ii =   particles(n)%x * ddx
+             jj = ( particles(n)%y + 0.5_wp * dy ) * ddy
+
+             x  = particles(n)%x - ii * dx
+             y  = particles(n)%y + ( 0.5_wp - jj ) * 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(kk,jj,ii)   + ( gg-bb ) *   &
+                                                   v(kk,jj,ii+1)   &
+                       + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *   &
+                                                   v(kk,jj+1,ii+1) &
+                       ) / ( 3.0_wp * gg ) - v_gtrans
+             IF ( kk+1 == nzt+1 )  THEN
+                v_int = v_int_l
+             ELSE
+                v_int_u = ( (gg-aa) * v(kk+1,jj,ii)   + (gg-bb) *  &
+                                                   v(kk+1,jj,ii+1) &
+                          + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) *  &
+                                                 v(kk+1,jj+1,ii+1) &
+                          ) / ( 3.0_wp * gg ) - v_gtrans
+                v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
+                                  * ( v_int_u - v_int_l )
+             ENDIF
+
+!
+!--          Same procedure for interpolation of the w velocity-component 
+!--          (adopt index i from v velocity-component)
+             jj = particles(n)%y * ddy
+             kk = particles(n)%z / dz
+
+             x  = particles(n)%x - ii * dx
+             y  = particles(n)%y - jj * 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(kk,jj,ii)   + ( gg-bb ) *   &
+                                                     w(kk,jj,ii+1) &
+                       + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *   &
+                                                   w(kk,jj+1,ii+1) &
+                       ) / ( 3.0_wp * gg )
+             IF ( kk+1 == nzt+1 )  THEN
+                w_int = w_int_l
+             ELSE
+                w_int_u = ( (gg-aa) * w(kk+1,jj,ii)   + (gg-bb) *  &
+                                                   w(kk+1,jj,ii+1) &
+                          + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) *  &
+                                                 w(kk+1,jj+1,ii+1) &
+                          ) / ( 3.0_wp * gg )
+                w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
+                                  * ( w_int_u - w_int_l )
+             ENDIF
+
+!
+!--          Change in radius due to collision
+             delta_r = effective_coll_efficiency / 3.0_wp          &
+                       * pi * sl_r3 * ddx * ddy / dz               &
+                       * SQRT( ( u_int - particles(n)%speed_x )**2 &
+                             + ( v_int - particles(n)%speed_y )**2 &
+                             + ( w_int - particles(n)%speed_z )**2 &
+                             ) * dt_3d
+!
+!--          Change in volume due to collision
+             delta_v = particles(n)%weight_factor                  &
+                       * ( ( particles(n)%radius + delta_r )**3    &
+                           - particles(n)%radius**3 )
+
+!
+!--          Check if collected particles provide enough LWC for
+!--          volume change of collector particle
+             IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
+
+                delta_r = ( ( sl_r3/particles(n)%weight_factor )               &
+                            + particles(n)%radius**3 )**( 1.0_wp / 3.0_wp )    &
+                            - particles(n)%radius
+
+                DO  is = n-1, psi, -1
+                   IF ( particles(is)%radius < particles(n)%radius )  THEN
+                      particles(is)%weight_factor = 0.0_wp
+                      particles(is)%particle_mask = .FALSE.
+                      deleted_particles = deleted_particles + 1
+                   ENDIF
                 ENDDO
 
-                psi = prt_start_index(k,j,i)
-                pse = psi + prt_count(k,j,i)-1
-
-!
-!--             Now apply the different kernels
-                IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
-                     use_kernel_tables )  THEN
-!
-!--                Fast method with pre-calculated efficiencies for
-!--                discrete radius- and dissipation-classes.
-!
-!--                Determine dissipation class index of this gridbox
-                   IF ( wang_kernel )  THEN
-                      eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
-                                    dissipation_classes ) + 1
-                      epsilon = diss(k,j,i)
-                   ELSE
-                      epsilon = 0.0
+             ELSE IF ( delta_v < sl_r3  .AND.  sl_r3 > 0.0_wp )  THEN
+
+                DO  is = n-1, psi, -1
+                   IF ( particles(is)%radius < particles(n)%radius &
+                        .AND.  sl_r3 > 0.0_wp )  THEN
+                      particles(is)%weight_factor =                &
+                                   ( ( particles(is)%weight_factor &
+                                   * ( particles(is)%radius**3 ) ) &
+                                   - ( delta_v                     &
+                                   * particles(is)%weight_factor   &
+                                   * ( particles(is)%radius**3 )   &
+                                   / sl_r3 ) )                     &
+                                   / ( particles(is)%radius**3 )
+
+                      IF ( particles(is)%weight_factor < 0.0_wp )  THEN
+                         WRITE( message_string, * ) 'negative ',   &
+                                         'weighting factor: ',     &
+                                         particles(is)%weight_factor
+                         CALL message( 'lpm_droplet_collision', &
+                                       'PA0039',                &
+                                       2, 2, -1, 6, 1 )
+                      ENDIF
                    ENDIF
-                   IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001 )  THEN
-                      eclass = 0   ! Hall kernel is used
-                   ELSE
-                      eclass = MIN( dissipation_classes, eclass )
-                   ENDIF
-
- !
-!--                Droplet collision are calculated using collision-coalescence
-!--                formulation proposed by Wang (see PALM documentation)
-!--                Since new radii after collision are defined by radii of all
-!--                droplets before collision, temporary fields for new radii and
-!--                weighting factors are needed
-                   ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
-
-                   rad = 0.0
-                   weight = 0.0
-
-                   DO  n = psi, pse, 1
-
-                      sum1 = 0.0
-                      sum2 = 0.0
-                      sum3 = 0.0
-
-                      rclass_l = particles(n)%class
-!
-!--                   Mass added due to collisions with smaller droplets
-                      DO  is = psi, n-1
-
-                         rclass_s = particles(is)%class
-                         sum1 = sum1 + ( particles(is)%radius**3 *            &
-                                         ckernel(rclass_l,rclass_s,eclass) *  &
-                                         particles(is)%weight_factor )
-
-                      ENDDO
-!
-!--                   Rate of collisions with larger droplets
-                      DO  is = n+1, pse
-
-                         rclass_s = particles(is)%class
-                         sum2 = sum2 + ( ckernel(rclass_l,rclass_s,eclass) *  &
-                                         particles(is)%weight_factor )
-
-                      ENDDO
-
-                      r3 = particles(n)%radius**3
-                      ddV = ddx * ddy / dz 
-                      is = prt_start_index(k,j,i)
-!
-!--                   Change of the current weighting factor
-                      sum3 = 1 - dt_3d * ddV *                                 &
-                                 ckernel(rclass_l,rclass_l,eclass) *           &
-                                 ( particles(n)%weight_factor - 1 ) * 0.5 -    &
-                             dt_3d * ddV * sum2
-                      weight(n-is+1) = particles(n)%weight_factor * sum3
-!
-!--                   Change of the current droplet radius
-                      rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
-                                       sum3 )**0.33333333333333_wp
-
-                      IF ( weight(n-is+1) < 0.0 )  THEN
-                         WRITE( message_string, * ) 'negative weighting',      &
-                                                    'factor: ', weight(n-is+1)
-                         CALL message( 'lpm_droplet_collision', 'PA0028',      &
-                                        2, 2, -1, 6, 1 )
-                      ENDIF
-
-                      ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
-                                                  * rad(n-is+1)**3
-
-                   ENDDO
-
-                   particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
-                   particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
-
-                   DEALLOCATE(rad, weight)
-
-                ELSEIF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
-                         .NOT. use_kernel_tables )  THEN
-!
-!--                Collision efficiencies are calculated for every new
-!--                grid box. First, allocate memory for kernel table.
-!--                Third dimension is 1, because table is re-calculated for
-!--                every new dissipation value.
-                   ALLOCATE( ckernel(prt_start_index(k,j,i):                 &
-                             prt_start_index(k,j,i)+prt_count(k,j,i)-1,      &
-                             prt_start_index(k,j,i):                         &
-                             prt_start_index(k,j,i)+prt_count(k,j,i)-1,1:1) )
-!
-!--                Now calculate collision efficiencies for this box
-                   CALL recalculate_kernel( i, j, k )
-
-!
-!--                Droplet collision are calculated using collision-coalescence
-!--                formulation proposed by Wang (see PALM documentation)
-!--                Since new radii after collision are defined by radii of all
-!--                droplets before collision, temporary fields for new radii and
-!--                weighting factors are needed
-                   ALLOCATE(rad(1:prt_count(k,j,i)), weight(1:prt_count(k,j,i)))
-
-                   rad = 0.0
-                   weight = 0.0
-
-                   DO  n = psi, pse, 1
-
-                      sum1 = 0.0
-                      sum2 = 0.0
-                      sum3 = 0.0
-!
-!--                   Mass added due to collisions with smaller droplets
-                      DO  is = psi, n-1
-                         sum1 = sum1 + ( particles(is)%radius**3 *            &
-                                         ckernel(n,is,1) *  &
-                                         particles(is)%weight_factor )
-                      ENDDO
-!
-!--                   Rate of collisions with larger droplets
-                      DO  is = n+1, pse
-                         sum2 = sum2 + ( ckernel(n,is,1) *  &
-                                         particles(is)%weight_factor )
-                      ENDDO
-
-                      r3 = particles(n)%radius**3
-                      ddV = ddx * ddy / dz 
-                      is = prt_start_index(k,j,i)
-!
-!--                   Change of the current weighting factor
-                      sum3 = 1 - dt_3d * ddV *                                 &
-                                 ckernel(n,n,1) *           &
-                                 ( particles(n)%weight_factor - 1 ) * 0.5 -    &
-                             dt_3d * ddV * sum2
-                      weight(n-is+1) = particles(n)%weight_factor * sum3
-!
-!--                   Change of the current droplet radius
-                      rad(n-is+1) = ( (r3 + dt_3d * ddV * (sum1 - sum2 * r3) )/&
-                                       sum3 )**0.33333333333333_wp
-
-                      IF ( weight(n-is+1) < 0.0 )  THEN
-                         WRITE( message_string, * ) 'negative weighting',      &
-                                                    'factor: ', weight(n-is+1)
-                         CALL message( 'lpm_droplet_collision', 'PA0037',      &
-                                        2, 2, -1, 6, 1 )
-                      ENDIF
-
-                      ql_vp(k,j,i) = ql_vp(k,j,i) + weight(n-is+1) &
-                                                  * rad(n-is+1)**3
-
-                   ENDDO
-
-                   particles(psi:pse)%radius = rad(1:prt_count(k,j,i))
-                   particles(psi:pse)%weight_factor = weight(1:prt_count(k,j,i))
-
-                   DEALLOCATE( rad, weight, ckernel )
-
-                ELSEIF ( palm_kernel )  THEN
-!
-!--                PALM collision kernel
-!
-!--                Calculate the mean radius of all those particles which
-!--                are of smaller size than the current particle and
-!--                use this radius for calculating the collision efficiency
-                   DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
-
-                      sl_r3 = 0.0
-                      sl_r4 = 0.0
-
-                      DO is = n-1, psi, -1
-                         IF ( particles(is)%radius < particles(n)%radius )  &
-                         THEN
-                            sl_r3 = sl_r3 + particles(is)%weight_factor     &
-                                            * particles(is)%radius**3
-                            sl_r4 = sl_r4 + particles(is)%weight_factor     &
-                                            * particles(is)%radius**4
-                         ENDIF
-                      ENDDO
-
-                      IF ( ( sl_r3 ) > 0.0 )  THEN
-                         mean_r = ( sl_r4 ) / ( sl_r3 )
-
-                         CALL collision_efficiency_rogers( mean_r,             &
-                                                    particles(n)%radius,       &
-                                                    effective_coll_efficiency )
-
-                      ELSE
-                         effective_coll_efficiency = 0.0
-                      ENDIF
-
-                      IF ( effective_coll_efficiency > 1.0  .OR.            &
-                           effective_coll_efficiency < 0.0 )                &
-                      THEN
-                         WRITE( message_string, * )  'collision_efficien' , &
-                                   'cy out of range:' ,effective_coll_efficiency
-                         CALL message( 'lpm_droplet_collision', 'PA0145', 2, &
-                                       2, -1, 6, 1 )
-                      ENDIF
-
-!
-!--                   Interpolation of ...
-                      ii = particles(n)%x * ddx
-                      jj = particles(n)%y * ddy
-                      kk = ( particles(n)%z + 0.5 * dz ) / dz
-
-                      x  = particles(n)%x - ii * dx
-                      y  = particles(n)%y - jj * 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
-
-                      ql_int_l = ( (gg-aa) * ql(kk,jj,ii)   + (gg-bb) *     &
-                                                           ql(kk,jj,ii+1)   &
-                                 + (gg-cc) * ql(kk,jj+1,ii) + ( gg-dd ) *   &
-                                                           ql(kk,jj+1,ii+1) &
-                                 ) / ( 3.0 * gg )
-
-                      ql_int_u = ( (gg-aa) * ql(kk+1,jj,ii)   + (gg-bb) *   &
-                                                         ql(kk+1,jj,ii+1)   &
-                                 + (gg-cc) * ql(kk+1,jj+1,ii) + (gg-dd) *   &
-                                                         ql(kk+1,jj+1,ii+1) &
-                                 ) / ( 3.0 * gg )
-
-                      ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *&
-                                          ( ql_int_u - ql_int_l )
-
-!
-!--                   Interpolate u velocity-component
-                      ii = ( particles(n)%x + 0.5 * dx ) * ddx
-                      jj =   particles(n)%y * ddy
-                      kk = ( particles(n)%z + 0.5 * dz ) / dz ! only if eqist
-
-                      IF ( ( particles(n)%z - zu(kk) ) > (0.5*dz) ) kk = kk+1
-
-                      x  = particles(n)%x + ( 0.5 - ii ) * dx
-                      y  = particles(n)%y - jj * 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(kk,jj,ii)   + (gg-bb) *       &
-                                                           u(kk,jj,ii+1)    &
-                                + (gg-cc) * u(kk,jj+1,ii) + (gg-dd) *       &
-                                                           u(kk,jj+1,ii+1)  &
-                                ) / ( 3.0 * gg ) - u_gtrans
-                      IF ( kk+1 == nzt+1 )  THEN
-                         u_int = u_int_l
-                      ELSE
-                         u_int_u = ( (gg-aa) * u(kk+1,jj,ii)   + (gg-bb) *  &
-                                                          u(kk+1,jj,ii+1)   &
-                                   + (gg-cc) * u(kk+1,jj+1,ii) + (gg-dd) *  &
-                                                          u(kk+1,jj+1,ii+1) &
-                                   ) / ( 3.0 * gg ) - u_gtrans
-                         u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz &
-                                           * ( u_int_u - u_int_l )
-                      ENDIF
-
-!
-!--                   Same procedure for interpolation of the v velocity-com-
-!--                   ponent (adopt index k from u velocity-component)
-                      ii =   particles(n)%x * ddx
-                      jj = ( particles(n)%y + 0.5 * dy ) * ddy
-
-                      x  = particles(n)%x - ii * dx
-                      y  = particles(n)%y + ( 0.5 - jj ) * 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(kk,jj,ii)   + ( gg-bb ) *   &
-                                                            v(kk,jj,ii+1)   &
-                                + ( gg-cc ) * v(kk,jj+1,ii) + ( gg-dd ) *   &
-                                                            v(kk,jj+1,ii+1) &
-                                ) / ( 3.0 * gg ) - v_gtrans
-                      IF ( kk+1 == nzt+1 )  THEN
-                         v_int = v_int_l
-                      ELSE
-                         v_int_u = ( (gg-aa) * v(kk+1,jj,ii)   + (gg-bb) *  &
-                                                            v(kk+1,jj,ii+1) &
-                                   + (gg-cc) * v(kk+1,jj+1,ii) + (gg-dd) *  &
-                                                          v(kk+1,jj+1,ii+1) &
-                                   ) / ( 3.0 * gg ) - v_gtrans
-                         v_int = v_int_l + ( particles(n)%z - zu(kk) ) / dz &
-                                           * ( v_int_u - v_int_l )
-                      ENDIF
-
-!
-!--                   Same procedure for interpolation of the w velocity-com-
-!--                   ponent (adopt index i from v velocity-component)
-                      jj = particles(n)%y * ddy
-                      kk = particles(n)%z / dz
-
-                      x  = particles(n)%x - ii * dx
-                      y  = particles(n)%y - jj * 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(kk,jj,ii)   + ( gg-bb ) *   &
-                                                              w(kk,jj,ii+1) &
-                                + ( gg-cc ) * w(kk,jj+1,ii) + ( gg-dd ) *   &
-                                                            w(kk,jj+1,ii+1) &
-                                ) / ( 3.0 * gg )
-                      IF ( kk+1 == nzt+1 )  THEN
-                         w_int = w_int_l
-                      ELSE
-                         w_int_u = ( (gg-aa) * w(kk+1,jj,ii)   + (gg-bb) *  &
-                                                            w(kk+1,jj,ii+1) &
-                                   + (gg-cc) * w(kk+1,jj+1,ii) + (gg-dd) *  &
-                                                          w(kk+1,jj+1,ii+1) &
-                                   ) / ( 3.0 * gg )
-                         w_int = w_int_l + ( particles(n)%z - zw(kk) ) / dz &
-                                           * ( w_int_u - w_int_l )
-                      ENDIF
-
-!
-!--                   Change in radius due to collision
-                      delta_r = effective_coll_efficiency / 3.0_wp          &
-                                * pi * sl_r3 * ddx * ddy / dz               &
-                                * SQRT( ( u_int - particles(n)%speed_x )**2 &
-                                      + ( v_int - particles(n)%speed_y )**2 &
-                                      + ( w_int - particles(n)%speed_z )**2 &
-                                      ) * dt_3d
-!
-!--                   Change in volume due to collision
-                      delta_v = particles(n)%weight_factor                  &
-                                * ( ( particles(n)%radius + delta_r )**3    &
-                                    - particles(n)%radius**3 )
-
-!
-!--                   Check if collected particles provide enough LWC for
-!--                   volume change of collector particle
-                      IF ( delta_v >= sl_r3  .AND.  sl_r3 > 0.0 )  THEN
-
-                         delta_r = ( ( sl_r3/particles(n)%weight_factor )           &
-                                     + particles(n)%radius**3 )**( 1.0_wp/3.0_wp )  &
-                                     - particles(n)%radius
-
-                         DO  is = n-1, psi, -1
-                            IF ( particles(is)%radius < &
-                                 particles(n)%radius )  THEN
-                               particles(is)%weight_factor = 0.0
-                               particle_mask(is) = .FALSE.
-                               deleted_particles = deleted_particles + 1
-                            ENDIF
-                         ENDDO
-
-                      ELSE IF ( delta_v < sl_r3  .AND.  sl_r3 > 0.0 )  THEN
-
-                         DO  is = n-1, psi, -1
-                            IF ( particles(is)%radius < particles(n)%radius &
-                                 .AND.  sl_r3 > 0.0 )  THEN
-                               particles(is)%weight_factor =                &
-                                            ( ( particles(is)%weight_factor &
-                                            * ( particles(is)%radius**3 ) ) &
-                                            - ( delta_v                     &
-                                            * particles(is)%weight_factor   &
-                                            * ( particles(is)%radius**3 )   &
-                                            / sl_r3 ) )                     &
-                                            / ( particles(is)%radius**3 )
-
-                               IF ( particles(is)%weight_factor < 0.0 )  THEN
-                                  WRITE( message_string, * ) 'negative ',   &
-                                                  'weighting factor: ',     &
-                                                  particles(is)%weight_factor
-                                  CALL message( 'lpm_droplet_collision', &
-                                                'PA0039',                &
-                                                2, 2, -1, 6, 1 )
-                               ENDIF
-                            ENDIF
-                         ENDDO
-
-                      ENDIF
-
-                      particles(n)%radius = particles(n)%radius + delta_r
-                      ql_vp(k,j,i) = ql_vp(k,j,i) +               &
-                                     particles(n)%weight_factor * &
-                                     ( particles(n)%radius**3 )
-                   ENDDO
-
-                ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
-                                              * particles(psi)%radius**3
-
-                ENDIF  ! collision kernel
-
-             ELSE IF ( prt_count(k,j,i) == 1 )  THEN
-
-                psi = prt_start_index(k,j,i)
-
-!
-!--             Calculate change of weighting factor due to self collision
-                IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
-                     use_kernel_tables )  THEN
-
-                   IF ( wang_kernel )  THEN
-                      eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
-                                    dissipation_classes ) + 1
-                      epsilon = diss(k,j,i)
-                   ELSE
-                      epsilon = 0.0
-                   ENDIF
-                   IF ( hall_kernel .OR. epsilon * 1.0E4_wp < 0.001 )  THEN
-                      eclass = 0   ! Hall kernel is used
-                   ELSE
-                      eclass = MIN( dissipation_classes, eclass )
-                   ENDIF
-
-                   ddV = ddx * ddy / dz 
-                   rclass_l = particles(psi)%class
-                   sum3 = 1 - dt_3d * ddV *                                 &
-                              ( ckernel(rclass_l,rclass_l,eclass) *         &
-                                ( particles(psi)%weight_factor-1 ) * 0.5 )
-
-                   particles(psi)%radius = ( particles(psi)%radius**3 / &
-                                             sum3 )**0.33333333333333_wp
-                   particles(psi)%weight_factor = particles(psi)%weight_factor &
-                                                  * sum3
-
-                ELSE IF ( ( hall_kernel .OR. wang_kernel )  .AND.  &
-                           .NOT. use_kernel_tables )  THEN
-!
-!--                Collision efficiencies are calculated for every new
-!--                grid box. First, allocate memory for kernel table.
-!--                Third dimension is 1, because table is re-calculated for
-!--                every new dissipation value.
-                   ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) )
-!
-!--                Now calculate collision efficiencies for this box
-                   CALL recalculate_kernel( i, j, k )
-
-                   ddV = ddx * ddy / dz 
-                   sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) *  &
-                                ( particles(psi)%weight_factor - 1 ) * 0.5 ) 
-
-                   particles(psi)%radius = ( particles(psi)%radius**3 / &
-                                             sum3 )**0.33333333333333_wp
-                   particles(psi)%weight_factor = particles(psi)%weight_factor &
-                                                  * sum3
-
-                   DEALLOCATE( ckernel )
-                ENDIF 
-
-               ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
-                                particles(psi)%radius**3
+                ENDDO
+
              ENDIF
 
-!
-!--          Check if condensation of LWC was conserved during collision
-!--          process
-             IF ( ql_v(k,j,i) /= 0.0 )  THEN
-                IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001  .OR.             &
-                     ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999 )  THEN
-                   WRITE( message_string, * ) 'LWC is not conserved during',&
-                                              ' collision! ',               &
-                                              'LWC after condensation: ',   &
-                                              ql_v(k,j,i),                  &
-                                              ' LWC after collision: ',     &
-                                              ql_vp(k,j,i)
-                   CALL message( 'lpm_droplet_collision', 'PA0040',         &
-                                 2, 2, -1, 6, 1 )
-                ENDIF
-             ENDIF
-
+             particles(n)%radius = particles(n)%radius + delta_r
+             ql_vp(k,j,i) = ql_vp(k,j,i) +               &
+                            particles(n)%weight_factor * &
+                            ( particles(n)%radius**3 )
           ENDDO
-       ENDDO
-    ENDDO
+
+          ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
+                                        * particles(psi)%radius**3
+
+       ENDIF  ! collision kernel
+
+    ELSE IF ( prt_count(k,j,i) == 1 )  THEN
+
+       psi = 1
+
+!
+!--    Calculate change of weighting factor due to self collision
+       IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
+            use_kernel_tables )  THEN
+
+          IF ( wang_kernel )  THEN
+             eclass = INT( diss(k,j,i) * 1.0E4_wp / 1000.0_wp * &
+                           dissipation_classes ) + 1
+             epsilon = diss(k,j,i)
+          ELSE
+             epsilon = 0.0_wp
+          ENDIF
+          IF ( hall_kernel  .OR.  epsilon * 1.0E4_wp < 0.001_wp )  THEN
+             eclass = 0   ! Hall kernel is used
+          ELSE
+             eclass = MIN( dissipation_classes, eclass )
+          ENDIF
+
+          ddV = ddx * ddy / dz 
+          rclass_l = particles(psi)%class
+          sum3 = 1 - dt_3d * ddV *                                 &
+                     ( ckernel(rclass_l,rclass_l,eclass) *         &
+                       ( particles(psi)%weight_factor-1 ) * 0.5_wp )
+
+          particles(psi)%radius = ( particles(psi)%radius**3 / &
+                                    sum3 )**0.33333333333333_wp
+          particles(psi)%weight_factor = particles(psi)%weight_factor &
+                                         * sum3
+
+       ELSE IF ( ( hall_kernel  .OR.  wang_kernel )  .AND.  &
+                  .NOT. use_kernel_tables )  THEN
+!
+!--       Collision efficiencies are calculated for every new
+!--       grid box. First, allocate memory for kernel table.
+!--       Third dimension is 1, because table is re-calculated for
+!--       every new dissipation value.
+          ALLOCATE( ckernel(psi:psi, psi:psi, 1:1) )
+!
+!--       Now calculate collision efficiencies for this box
+          CALL recalculate_kernel( i, j, k )
+
+          ddV = ddx * ddy / dz 
+          sum3 = 1 - dt_3d * ddV * ( ckernel(psi,psi,1) *  &
+                       ( particles(psi)%weight_factor - 1 ) * 0.5_wp ) 
+
+          particles(psi)%radius = ( particles(psi)%radius**3 / &
+                                    sum3 )**0.33333333333333_wp
+          particles(psi)%weight_factor = particles(psi)%weight_factor &
+                                         * sum3
+
+          DEALLOCATE( ckernel )
+       ENDIF 
+
+      ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
+                       particles(psi)%radius**3
+    ENDIF
+
+!
+!-- Check if condensation of LWC was conserved during collision process
+    IF ( ql_v(k,j,i) /= 0.0_wp )  THEN
+       IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001_wp  .OR.             &
+            ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999_wp )  THEN
+          WRITE( message_string, * ) 'LWC is not conserved during',&
+                                     ' collision! ',               &
+                                     'LWC after condensation: ',   &
+                                     ql_v(k,j,i),                  &
+                                     ' LWC after collision: ',     &
+                                     ql_vp(k,j,i)
+          CALL message( 'lpm_droplet_collision', 'PA0040',         &
+                        2, 2, -1, 6, 1 )
+       ENDIF
+    ENDIF
 
     CALL cpu_log( log_point_s(43), 'lpm_droplet_coll', 'stop' )
 
-
  END SUBROUTINE lpm_droplet_collision
