Index: palm/trunk/SOURCE/lpm_collision_kernels.f90
===================================================================
--- palm/trunk/SOURCE/lpm_collision_kernels.f90	(revision 1070)
+++ palm/trunk/SOURCE/lpm_collision_kernels.f90	(revision 1071)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-!
+! Bugfix: collision efficiencies for Hall kernel should not be < 1.0E-20
 !
 ! Former revisions:
@@ -607,5 +607,5 @@
 
 !------------------------------------------------------------------------------!
-! Calculation of collision efficencies for the Hall kernel
+! Calculation of collision efficiencies for the Hall kernel
 !------------------------------------------------------------------------------!
     SUBROUTINE effic
@@ -728,8 +728,9 @@
                 ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
                 ec(j,i) = MIN( ek, 1.0 )
-             ENDIF 
+             ENDIF
+
+             IF ( ec(j,i) < 1.0E-20 )  ec(j,i) = 0.0
 
              ec(i,j) = ec(j,i)
-             IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
 
           ENDDO
Index: palm/trunk/SOURCE/lpm_droplet_collision.f90
===================================================================
--- palm/trunk/SOURCE/lpm_droplet_collision.f90	(revision 1070)
+++ 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
Index: palm/trunk/SOURCE/lpm_droplet_condensation.f90
===================================================================
--- palm/trunk/SOURCE/lpm_droplet_condensation.f90	(revision 1070)
+++ palm/trunk/SOURCE/lpm_droplet_condensation.f90	(revision 1071)
@@ -20,5 +20,14 @@
 ! Current revisions:
 ! ------------------
-! 
+! Ventilation effect for evaporation of large droplets included
+! Check for unreasonable results included in calculation of Rosenbrock method
+! since physically unlikely results were observed and for the same
+! reason the first internal time step in Rosenbrock method should be < 1.0E02 in
+! case of evaporation
+! Unnecessary calculation of ql_int removed
+! Unnecessary calculations in Rosenbrock method (d2rdt2, drdt_m, dt_ros_last)
+! removed
+! Bugfix: factor in calculation of surface tension changed from 0.00155 to
+! 0.000155
 !
 ! Former revisions:
@@ -55,24 +64,23 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, internal_timestep_count, j, jtry, k, n
+    INTEGER ::  i, internal_timestep_count, j, jtry, k, n, ros_count
 
     INTEGER, PARAMETER ::  maxtry = 40
 
+    LOGICAL ::  repeat
+
     REAL ::  aa, afactor, arg, bb, cc, dd, ddenom, delta_r, drdt, drdt_ini,    &
-             drdt_m, dt_ros, dt_ros_last, dt_ros_next, dt_ros_sum,             &
-             dt_ros_sum_ini, d2rdtdr, d2rdt2, errmax, err_ros, g1, g2, g3, g4, &
-             e_a, e_s, gg, new_r, p_int, pt_int, pt_int_l, pt_int_u, q_int,    &
-             q_int_l, q_int_u, ql_int, ql_int_l, ql_int_u, r_ros, r_ros_ini,   &
-             sigma, t_int, x, y
+             dt_ros, dt_ros_next, dt_ros_sum, dt_ros_sum_ini, d2rdtdr, errmax, &
+             err_ros, g1, g2, g3, g4, e_a, e_s, gg, new_r, p_int, pt_int,      &
+             pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l,    &
+             ql_int_u, r_ros, r_ros_ini, sigma, t_int, x, y, re_p
 
 !
 !-- Parameters for Rosenbrock method
     REAL, PARAMETER ::  a21 = 2.0, a31 = 48.0/25.0, a32 = 6.0/25.0,        &
-                        a2x = 1.0, a3x = 3.0/5.0, b1 = 19.0/9.0, b2 = 0.5, &
-                        b3 = 25.0/108.0, b4 = 125.0/108.0, c21 = -8.0,     &
-                        c31 = 372.0/25.0, c32 = 12.0/5.0,                  &
-                        c41 = -112.0/125.0, c42 = -54.0/125.0,             &
-                        c43 = -2.0/5.0, c1x = 0.5, c2x= -3.0/2.0,          &
-                        c3x = 121.0/50.0, c4x = 29.0/250.0,                &
+                        b1 = 19.0/9.0, b2 = 0.5, b3 = 25.0/108.0,          &
+                        b4 = 125.0/108.0, c21 = -8.0, c31 = 372.0/25.0,    &
+                        c32 = 12.0/5.0, c41 = -112.0/125.0,                &
+                        c42 = -54.0/125.0, c43 = -2.0/5.0,                 &
                         errcon = 0.1296, e1 = 17.0/54.0, e2 = 7.0/36.0,    &
                         e3 = 0.0, e4 = 125.0/108.0, gam = 0.5, grow = 1.5, &
@@ -121,15 +129,4 @@
                            ( q_int_u - q_int_l )
 
-       ql_int_l = ( ( gg - aa ) * ql(k,j,i)   + ( gg - bb ) * ql(k,j,i+1)   &
-                  + ( gg - cc ) * ql(k,j+1,i) + ( gg - dd ) * ql(k,j+1,i+1) &
-                  ) / ( 3.0 * gg )
-
-       ql_int_u = ( ( gg-aa ) * ql(k+1,j,i)   + ( gg-bb ) * ql(k+1,j,i+1)   &
-                  + ( gg-cc ) * ql(k+1,j+1,i) + ( gg-dd ) * ql(k+1,j+1,i+1) &
-                  ) / ( 3.0 * gg )
-
-       ql_int = ql_int_l + ( particles(n)%z - zu(k) ) / dz * &
-                             ( ql_int_u - ql_int_l )
-
 !
 !--    Calculate real temperature and saturation vapor pressure
@@ -151,9 +148,33 @@
 !
 !--    Change in radius by condensation/evaporation
-       IF ( particles(n)%radius >= 1.0E-6  .OR.  &
-            .NOT. curvature_solution_effects )  THEN
-!
-!--       Approximation for large radii, where curvature and solution
-!--       effects can be neglected
+       IF ( particles(n)%radius >= 4.0E-5  .AND.  e_a/e_s < 1.0 )  THEN
+!
+!--       Approximation for large radii, where curvature and solution effects
+!--       can be neglected but ventilation effect has to be included in case of
+!--       evaporation.
+!--       First calculate the droplet's Reynolds number
+          re_p = 2.0 * particles(n)%radius * ABS( particles(n)%speed_z ) / &
+                 molecular_viscosity
+!
+!--       Ventilation coefficient after Rogers and Yau, 1989
+          IF ( re_p > 2.5 )  THEN
+             afactor = 0.78 + 0.28 * SQRT( re_p )
+          ELSE
+             afactor = 1.0 + 0.09 * re_p
+          ENDIF
+
+          arg = particles(n)%radius**2 + 2.0 * dt_3d * afactor *              &
+                           ( e_a / e_s - 1.0 ) /                              &
+                           ( ( l_d_rv / t_int - 1.0 ) * l_v * rho_l / t_int / &
+                             thermal_conductivity_l +                         &
+                             rho_l * r_v * t_int / diff_coeff_l / e_s )
+
+          new_r = SQRT( arg )
+
+       ELSEIF ( particles(n)%radius >= 1.0E-6  .OR.  &
+                .NOT. curvature_solution_effects )  THEN
+!
+!--       Approximation for larger radii in case that curvature and solution
+!--       effects are neglected and ventilation effects does not play a role
           arg = particles(n)%radius**2 + 2.0 * dt_3d *                        &
                            ( e_a / e_s - 1.0 ) /                              &
@@ -178,145 +199,175 @@
 !--       For larger radii the simple analytic method (see ELSE) gives
 !--       almost the same results.
-!
-!--       Surface tension after (Straka, 2009)
-          sigma = 0.0761 - 0.00155 * ( t_int - 273.15 )
-
-          r_ros = particles(n)%radius
-          dt_ros_sum  = 0.0      ! internal integrated time (s)
-          internal_timestep_count = 0
-
-          ddenom  = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) +  &
-                            ( l_v / ( r_v * t_int ) - 1.0 ) *               &
-                            rho_l * l_v / ( thermal_conductivity_l * t_int )&
-                          )
-              
-          afactor = 2.0 * sigma / ( rho_l * r_v * t_int )
-
-          IF ( particles(n)%rvar3 == -9999999.9 )  THEN
-!
-!--          First particle timestep. Derivative has to be calculated.
-             drdt_m  = ddenom / r_ros * ( e_a / e_s - 1.0 -    &
-                                          afactor / r_ros +    &
+
+          ros_count = 0
+          repeat = .TRUE.
+!
+!--       Carry out the Rosenbrock algorithm. In case of unreasonable results
+!--       the switch "repeat" will be set true and the algorithm will be carried
+!--       out again with the internal time step set to its initial (small) value.
+!--       Unreasonable results may occur if the external conditions, especially the
+!--       supersaturation, has significantly changed compared to the last PALM
+!--       timestep.
+          DO WHILE ( repeat )
+
+             repeat = .FALSE.
+!
+!--          Surface tension after (Straka, 2009)
+             sigma = 0.0761 - 0.000155 * ( t_int - 273.15 )
+
+             r_ros = particles(n)%radius
+             dt_ros_sum  = 0.0      ! internal integrated time (s)
+             internal_timestep_count = 0
+
+             ddenom  = 1.0 / ( rho_l * r_v * t_int / ( e_s * diff_coeff_l ) +  &
+                               ( l_v / ( r_v * t_int ) - 1.0 ) *               &
+                               rho_l * l_v / ( thermal_conductivity_l * t_int )&
+                             )
+
+             afactor = 2.0 * sigma / ( rho_l * r_v * t_int )
+
+!
+!--          Take internal time step values from the end of last PALM time step
+             dt_ros_next = particles(n)%rvar1
+
+!
+!--          Internal time step should not be > 1.0E-2 in case of evaporation
+!--          because larger values may lead to secondary solutions which are
+!--          physically unlikely
+             IF ( dt_ros_next > 1.0E-2  .AND.  e_a/e_s < 1.0 )  THEN
+                dt_ros_next = 1.0E-3
+             ENDIF
+!
+!--          If calculation of Rosenbrock method is repeated due to unreasonalble
+!--          results during previous try the initial internal time step has to be
+!--          reduced
+             IF ( ros_count > 1 )  THEN
+                dt_ros_next = dt_ros_next - ( 0.2 * dt_ros_next )
+             ELSEIF ( ros_count > 5 )  THEN
+!
+!--             Prevent creation of infinite loop
+                message_string = 'ros_count > 5 in Rosenbrock method'
+                CALL message( 'lpm_droplet_condensation', 'PA0018', 2, 2, &
+                               0, 6, 0 )
+             ENDIF
+
+!
+!--          Internal time step must not be larger than PALM time step
+             dt_ros = MIN( dt_ros_next, dt_3d )
+!
+!--          Integrate growth equation in time unless PALM time step is reached
+             DO WHILE ( dt_ros_sum < dt_3d )
+
+                internal_timestep_count = internal_timestep_count + 1
+
+!
+!--             Derivative at starting value
+                drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + &
                                           bfactor / r_ros**3 )
-          ELSE
-!
-!--          Take value from last PALM timestep
-             drdt_m = particles(n)%rvar3
-          ENDIF
-!
-!--       Take internal timestep values from the end of last PALM timestep
-          dt_ros_last = particles(n)%rvar1
-          dt_ros_next = particles(n)%rvar2
-!
-!--       Internal timestep must not be larger than PALM timestep
-          dt_ros = MIN( dt_ros_next, dt_3d )
-!
-!--       Integrate growth equation in time unless PALM timestep is reached
-          DO WHILE ( dt_ros_sum < dt_3d )
-
-             internal_timestep_count = internal_timestep_count + 1
-
-!
-!--          Derivative at starting value
-             drdt = ddenom / r_ros * ( e_a / e_s - 1.0 - afactor / r_ros + &
-                                       bfactor / r_ros**3 )
-             drdt_ini       = drdt
-             dt_ros_sum_ini = dt_ros_sum
-             r_ros_ini      = r_ros
-
-!
-!--          Calculate time derivative of dr/dt
-             d2rdt2 = ( drdt - drdt_m ) / dt_ros_last 
-
-!
-!--          Calculate radial derivative of dr/dt
-             d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + &
-                                  2.0 * afactor / r_ros**3 -     &
-                                  4.0 * bfactor / r_ros**5 )
-!
-!--          Adjust stepsize unless required accuracy is reached
-             DO  jtry = 1, maxtry+1
-
-                IF ( jtry == maxtry+1 )  THEN
-                   message_string = 'maxtry > 40 in Rosenbrock method'
-                   CALL message( 'lpm_droplet_condensation', 'PA0347', 2, 2, &
-                                 0, 6, 0 )
+                drdt_ini       = drdt
+                dt_ros_sum_ini = dt_ros_sum
+                r_ros_ini      = r_ros
+
+!
+!--             Calculate radial derivative of dr/dt
+                d2rdtdr = ddenom * ( ( 1.0 - e_a/e_s ) / r_ros**2 + &
+                                     2.0 * afactor / r_ros**3 -     &
+                                     4.0 * bfactor / r_ros**5 )
+!
+!--             Adjust stepsize unless required accuracy is reached
+                DO  jtry = 1, maxtry+1
+
+                   IF ( jtry == maxtry+1 )  THEN
+                      message_string = 'maxtry > 40 in Rosenbrock method'
+                      CALL message( 'lpm_droplet_condensation', 'PA0347', 2, 2, &
+                                    0, 6, 0 )
+                   ENDIF
+
+                   aa    = 1.0 / ( gam * dt_ros ) - d2rdtdr
+                   g1    = drdt_ini / aa
+                   r_ros = r_ros_ini + a21 * g1
+                   drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
+                                              afactor / r_ros + &
+                                              bfactor / r_ros**3 )
+
+                   g2    = ( drdt + c21 * g1 / dt_ros )&
+                           / aa
+                   r_ros = r_ros_ini + a31 * g1 + a32 * g2
+                   drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
+                                              afactor / r_ros + &
+                                              bfactor / r_ros**3 )
+
+                   g3    = ( drdt +  &
+                             ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
+                   g4    = ( drdt +  &
+                             ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
+                   r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
+
+                   dt_ros_sum = dt_ros_sum_ini + dt_ros
+
+                   IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
+                      message_string = 'zero stepsize in Rosenbrock method'
+                      CALL message( 'lpm_droplet_condensation', 'PA0348', 2, 2, &
+                                    0, 6, 0 )
+                   ENDIF
+!
+!--                Calculate error
+                   err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4
+                   errmax  = 0.0
+                   errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
+!
+!--                Leave loop if accuracy is sufficient, otherwise try again
+!--                with a reduced stepsize
+                   IF ( errmax <= 1.0 )  THEN
+                      EXIT
+                   ELSE
+                      dt_ros = SIGN( MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), &
+                                          shrnk * ABS( dt_ros ) ), dt_ros )
+                   ENDIF
+
+                ENDDO  ! loop for stepsize adjustment
+
+!
+!--             Calculate next internal time step
+                IF ( errmax > errcon )  THEN
+                   dt_ros_next = 0.9 * dt_ros * errmax**pgrow
+                ELSE
+                   dt_ros_next = grow * dt_ros
                 ENDIF
 
-                aa    = 1.0 / ( gam * dt_ros ) - d2rdtdr
-                g1    = ( drdt_ini + dt_ros * c1x * d2rdt2 ) / aa
-                r_ros = r_ros_ini + a21 * g1
-                drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
-                                           afactor / r_ros + &
-                                           bfactor / r_ros**3 )
-
-                g2    = ( drdt + dt_ros * c2x * d2rdt2 + c21 * g1 / dt_ros )&
-                        / aa
-                r_ros = r_ros_ini + a31 * g1 + a32 * g2
-                drdt  = ddenom / r_ros * ( e_a / e_s - 1.0 - &
-                                           afactor / r_ros + &
-                                           bfactor / r_ros**3 )
-
-                g3    = ( drdt + dt_ros * c3x * d2rdt2 + &
-                          ( c31 * g1 + c32 * g2 ) / dt_ros ) / aa
-                g4    = ( drdt + dt_ros * c4x * d2rdt2 + &
-                          ( c41 * g1 + c42 * g2 + c43 * g3 ) / dt_ros ) / aa
-                r_ros = r_ros_ini + b1 * g1 + b2 * g2 + b3 * g3 + b4 * g4
-
-                dt_ros_sum = dt_ros_sum_ini + dt_ros
-
-                IF ( dt_ros_sum == dt_ros_sum_ini )  THEN
-                   message_string = 'zero stepsize in Rosenbrock method'
-                   CALL message( 'lpm_droplet_condensation', 'PA0348', 2, 2, &
-                                 0, 6, 0 )
+!
+!--             Estimated time step is reduced if the PALM time step is exceeded
+                IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
+                   dt_ros = dt_3d - dt_ros_sum
+                ELSE
+                   dt_ros = dt_ros_next
                 ENDIF
-!
-!--             Calculate error
-                err_ros = e1*g1 + e2*g2 + e3*g3 + e4*g4 
-                errmax  = 0.0
-                errmax  = MAX( errmax, ABS( err_ros / r_ros_ini ) ) / eps_ros
-!
-!--             Leave loop if accuracy is sufficient, otherwise try again
-!--             with a reduced stepsize
-                IF ( errmax <= 1.0 )  THEN
-                   EXIT
-                ELSE
-                   dt_ros = SIGN( MAX( ABS( 0.9 * dt_ros * errmax**pshrnk ), &
-                                       shrnk * ABS( dt_ros ) ), dt_ros )
-                ENDIF
-
-             ENDDO  ! loop for stepsize adjustment
-
-!
-!--          Calculate next internal timestep
-             IF ( errmax > errcon )  THEN
-                dt_ros_next = 0.9 * dt_ros * errmax**pgrow
-             ELSE
-                dt_ros_next = grow * dt_ros
+
+             ENDDO
+!
+!--          Store internal time step value for next PALM step
+             particles(n)%rvar1 = dt_ros_next
+
+             new_r = r_ros
+!
+!--          Radius should not fall below 1E-8 because Rosenbrock method may
+!--          lead to errors otherwise
+             new_r = MAX( new_r, 1.0E-8 )
+!
+!--          Check if calculated droplet radius change is reasonable since in
+!--          case of droplet evaporation the Rosenbrock method may lead to
+!--          secondary solutions which are physically unlikely.
+!--          Due to the solution effect the droplets may grow for relative
+!--          humidities below 100%, but change of radius should not be too large.
+!--          In case of unreasonable droplet growth the Rosenbrock method is
+!--          recalculated using a smaller initial time step.
+!--          Limiting values are tested for droplets down to 1.0E-7
+             IF ( new_r - particles(n)%radius >= 3.0E-7  .AND.  &
+                  e_a/e_s < 0.97 )  THEN
+                ros_count = ros_count + 1
+                repeat = .TRUE.
              ENDIF
 
-!
-!--          Estimated timestep is reduced if the PALM time step is exceeded
-             dt_ros_last = dt_ros
-             IF ( ( dt_ros_next + dt_ros_sum ) >= dt_3d )  THEN
-                dt_ros = dt_3d - dt_ros_sum
-             ELSE
-                dt_ros = dt_ros_next
-             ENDIF
-
-             drdt_m = drdt
-
-          ENDDO
-!
-!--       Store derivative and internal timestep values for next PALM step
-          particles(n)%rvar1 = dt_ros_last
-          particles(n)%rvar2 = dt_ros_next
-          particles(n)%rvar3 = drdt_m
-
-          new_r = r_ros
-!
-!--       Radius should not fall below 1E-8 because Rosenbrock method may
-!--       lead to errors otherwise
-          new_r = MAX( new_r, 1.0E-8 )
+          ENDDO    ! Rosenbrock method
 
        ENDIF
