Index: palm/trunk/SOURCE/lpm_droplet_collision.f90
===================================================================
--- palm/trunk/SOURCE/lpm_droplet_collision.f90	(revision 1037)
+++ palm/trunk/SOURCE/lpm_droplet_collision.f90	(revision 1071)
@@ -20,5 +20,8 @@
 ! Current revisions:
 ! ------------------
-! 
+! Calculation of Hall and Wang kernel now uses collision-coalescence formulation
+! proposed by Wang instead of the continuous collection equation (for more
+! information about new method see PALM documentation)
+! Bugfix: message identifiers added
 !
 ! Former revisions:
@@ -35,5 +38,5 @@
 ! Description:
 ! ------------
-! Calculates chang in droplet radius by collision. Droplet collision is
+! Calculates change in droplet radius by collision. Droplet collision is
 ! calculated for each grid box seperately. Collision is parameterized by
 ! using collision kernels. Three different kernels are available:
@@ -74,7 +77,8 @@
              mean_r, ql_int, ql_int_l, ql_int_u, u_int, u_int_l, u_int_u,     &
              v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u, sl_r3, sl_r4,  &
-             x, y
+             x, y, sum1, sum2, sum3, r3, ddV
 
     TYPE(particle_type) ::  tmp_particle
+    REAL, DIMENSION(:), ALLOCATABLE :: rad, weight
 
 
@@ -138,91 +142,73 @@
                    ENDIF
 
-                   DO  n = pse, psi+1, -1
-
-                      integral = 0.0
-                      lw_max   = 0.0
+ !
+!--                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
 !
-!--                   Calculate growth of collector particle radius using all
-!--                   droplets smaller than current droplet 
+!--                   Mass added due to collisions with smaller droplets
                       DO  is = psi, n-1
 
                          rclass_s = particles(is)%class
-                         integral = integral +                             &
-                                    ( particles(is)%radius**3 *            &
-                                      ckernel(rclass_l,rclass_s,eclass) *  &
-                                      particles(is)%weight_factor )
-!
-!--                      Calculate volume of liquid water of the collected
-!--                      droplets which is the maximum liquid water available
-!--                      for droplet growth   
-                         lw_max = lw_max + ( particles(is)%radius**3 * &
-                                             particles(is)%weight_factor )
+                         sum1 = sum1 + ( particles(is)%radius**3 *            &
+                                         ckernel(rclass_l,rclass_s,eclass) *  &
+                                         particles(is)%weight_factor )
 
                       ENDDO
-
-!
-!--                   Change in radius of collector droplet due to collision
-                      delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
-                                integral * dt_3d * ddx * ddy / dz 
-
-!
-!--                   Change in volume of collector droplet due to collision
-                      delta_v = particles(n)%weight_factor                  &
-                                * ( ( particles(n)%radius + delta_r )**3    &
-                                    - particles(n)%radius**3 )
-
-                      IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
-!-- replace by message call
-                         PRINT*, 'Particle has grown to much because the',  &
-                                 ' change of volume of particle is larger', &
-                                 ' than liquid water available!'
-
-                      ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
-!-- can this case really happen??
-                         DO  is = psi, n-1
-
-                            particles(is)%weight_factor = 0.0
-                            particle_mask(is)  = .FALSE.
-                            deleted_particles = deleted_particles + 1
-
-                         ENDDO
-
-                      ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
-!
-!--                      Calculate new weighting factor of collected droplets
-                         DO  is = psi, n-1
-
-                            rclass_s = particles(is)%class
-                            particles(is)%weight_factor = &
-                               particles(is)%weight_factor - &
-                               ( ( ckernel(rclass_l,rclass_s,eclass) * particles(is)%weight_factor &
-                                 / integral ) * delta_v )
-
-                            IF ( particles(is)%weight_factor < 0.0 )  THEN
-                                  WRITE( message_string, * ) 'negative ',   &
-                                                     'weighting factor: ',  &
-                                                  particles(is)%weight_factor
-                                  CALL message( 'lpm_droplet_collision', '', &
-                                                2, 2, -1, 6, 1 )
-
-                            ELSEIF ( particles(is)%weight_factor == 0.0 )  &
-                            THEN
-
-                                particles(is)%weight_factor = 0.0
-                                particle_mask(is)  = .FALSE.
-                                deleted_particles = deleted_particles + 1
-
-                            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
+!
+!--                   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
+
+                      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.  &
@@ -241,89 +227,65 @@
                    CALL recalculate_kernel( i, j, k )
 
-                   DO  n = pse, psi+1, -1
-
-                      integral = 0.0
-                      lw_max   = 0.0
-!
-!--                   Calculate growth of collector particle radius using all
-!--                   droplets smaller than current droplet 
+!
+!--                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
-
-                         integral = integral + particles(is)%radius**3 * &
-                                               ckernel(n,is,1) *         &
-                                               particles(is)%weight_factor
-!
-!--                      Calculate volume of liquid water of the collected
-!--                      droplets which is the maximum liquid water available
-!--                      for droplet growth   
-                         lw_max = lw_max + ( particles(is)%radius**3 * &
-                                             particles(is)%weight_factor )
-
+                         sum1 = sum1 + ( particles(is)%radius**3 *            &
+                                         ckernel(n,is,1) *  &
+                                         particles(is)%weight_factor )
                       ENDDO
-
-!
-!--                   Change in radius of collector droplet due to collision
-                      delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
-                                integral * dt_3d * ddx * ddy / dz 
-
-!
-!--                   Change in volume of collector droplet due to collision
-                      delta_v = particles(n)%weight_factor                  &
-                                * ( ( particles(n)%radius + delta_r )**3    &
-                                    - particles(n)%radius**3 )
-
-                      IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
-!-- replace by message call
-                         PRINT*, 'Particle has grown to much because the',  &
-                                 ' change of volume of particle is larger', &
-                                 ' than liquid water available!'
-
-                      ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 ) THEN
-!-- can this case really happen??
-                         DO  is = psi, n-1
-
-                            particles(is)%weight_factor = 0.0
-                            particle_mask(is)  = .FALSE.
-                            deleted_particles = deleted_particles + 1
-
-                         ENDDO
-
-                      ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
-!
-!--                      Calculate new weighting factor of collected droplets
-                         DO  is = psi, n-1
-
-                            particles(is)%weight_factor =                   &
-                                   particles(is)%weight_factor -            &
-                                   ( ckernel(n,is,1) / integral * delta_v * &
-                                     particles(is)%weight_factor )
-
-                            IF ( particles(is)%weight_factor < 0.0 )  THEN
-                                  WRITE( message_string, * ) 'negative ',   &
-                                                     'weighting factor: ',  &
-                                                  particles(is)%weight_factor
-                                  CALL message( 'lpm_droplet_collision', '', &
-                                                2, 2, -1, 6, 1 )
-
-                            ELSEIF ( particles(is)%weight_factor == 0.0 )  &
-                            THEN
-
-                                particles(is)%weight_factor = 0.0
-                                particle_mask(is)  = .FALSE.
-                                deleted_particles = deleted_particles + 1
-
-                            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
+!
+!--                   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
+
+                      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
 
-                   DEALLOCATE( ckernel )
+                   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
@@ -543,5 +505,6 @@
                                                   'weighting factor: ',     &
                                                   particles(is)%weight_factor
-                                  CALL message( 'lpm_droplet_collision', '', &
+                                  CALL message( 'lpm_droplet_collision', &
+                                                'PA0039',                &
                                                 2, 2, -1, 6, 1 )
                                ENDIF
@@ -557,14 +520,67 @@
                    ENDDO
 
-                ENDIF  ! collision kernel
-
                 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)
-                ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
+
+!
+!--             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 / 1000.0 * &
+                                    dissipation_classes ) + 1
+                      epsilon = diss(k,j,i)
+                   ELSE
+                      epsilon = 0.0
+                   ENDIF
+                   IF ( hall_kernel .OR. epsilon * 1.0E4 < 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
+                   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
+                   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
@@ -582,5 +598,6 @@
                                               ' LWC after collision: ',     &
                                               ql_vp(k,j,i)
-                   CALL message( 'lpm_droplet_collision', '', 2, 2, -1, 6, 1 )
+                   CALL message( 'lpm_droplet_collision', 'PA0040',         &
+                                 2, 2, -1, 6, 1 )
                 ENDIF
              ENDIF
