Index: /palm/trunk/SCRIPTS/subjob
===================================================================
--- /palm/trunk/SCRIPTS/subjob	(revision 798)
+++ /palm/trunk/SCRIPTS/subjob	(revision 799)
@@ -140,5 +140,5 @@
     # VARIABLENVEREINBARUNGEN + DEFAULTWERTE
  delete_dayfile=false
- email_notification=none 
+ email_notification=none
  group_number=none
  locat=normal
@@ -250,4 +250,5 @@
      (rte*)                  local_addres=133.5.185.60;   local_host=lcrte;;
      (scirocco)              local_addres=172.20.25.41;   local_host=lcmuk;;
+     (solano)                local_addres=130.75.105.110; local_host=lcmuk;;
      (sun1|sun2)             local_addres=130.75.6.1;     local_host=unics;;
      (sx-*)                  local_addres=172.16.1.131;   local_host=necriam;;
@@ -885,5 +886,5 @@
        feature_directive=""
     fi
-    
+
     if [[ $queue = dataq || $queue = permq ]]
     then
@@ -952,5 +953,5 @@
     then
        cat > $job_to_send << %%END%%
-#!/bin/ksh 
+#!/bin/ksh
 #PBS -S /bin/ksh
 #PBS -N $job_name
@@ -994,9 +995,9 @@
 
        cat > $job_to_send << %%END%%
-#!/bin/ksh 
+#!/bin/ksh
 # @\$-o $remote_dayfile
 # @\$-eo -oi
 # @\$-lP 16
-# @\$-lp 1    
+# @\$-lp 1
 # @\$-lm 28gb  -llm unlimited -ls unlimited
 # @\$-q $queue
@@ -1021,5 +1022,5 @@
     then
        cat > $job_to_send << %%END%%
-#!/bin/ksh 
+#!/bin/ksh
 #PBS -S /bin/ksh
 #PBS -N $job_name
@@ -1442,7 +1443,7 @@
           if [[ $queue = default ]]
           then
-             eval  $submcom  $job_on_remhost 
+             eval  $submcom  $job_on_remhost
           else
-             eval  $submcom  -q $queue  $job_on_remhost 
+             eval  $submcom  -q $queue  $job_on_remhost
           fi
        else
Index: /palm/trunk/SOURCE/advec_particles.f90
===================================================================
--- /palm/trunk/SOURCE/advec_particles.f90	(revision 798)
+++ /palm/trunk/SOURCE/advec_particles.f90	(revision 799)
@@ -4,6 +4,6 @@
 ! Current revisions:
 ! ------------------
-! particle arrays (particles, particles_temp) implemented as pointers in
-! order to speed up sorting (see routine sort_particles)
+! Implementation of Wang collision kernel and corresponding new parameter 
+! wang_collision_kernel
 !
 ! Former revisions:
@@ -11,4 +11,8 @@
 ! $Id$
 !
+! 792 2011-12-01 raasch
+! particle arrays (particles, particles_temp) implemented as pointers in
+! order to speed up sorting (see routine sort_particles)
+!
 ! 759 2011-09-15 13:58:31Z raasch
 ! Splitting of parallel I/O (routine write_particles)
@@ -17,5 +21,4 @@
 ! Declaration of de_dx, de_dy, de_dz adapted to additional 
 ! ghost points. Furthermore the calls of exchange_horiz were modified.
-! 
 !
 ! 622 2010-12-10 08:08:13Z raasch
@@ -127,10 +130,11 @@
     USE random_function_mod
     USE statistics
+    USE wang_kernel_mod
 
     IMPLICIT NONE
 
     INTEGER ::  agp, deleted_particles, deleted_tails, i, ie, ii, inc, is, j,  &
-                jj, js, k, kk, kw, m, n, nc, nn, num_gp, psi, tlength,         &
-                trlp_count, trlp_count_sum, trlp_count_recv,                   &
+                jj, js, k, kk, kw, m, n, nc, nd, nn, num_gp, pse, psi,         &
+                tlength, trlp_count, trlp_count_sum, trlp_count_recv,          &
                 trlp_count_recv_sum, trlpt_count, trlpt_count_recv,            &
                 trnp_count, trnp_count_sum, trnp_count_recv,                   &
@@ -139,5 +143,5 @@
                 trrp_count_recv_sum, trrpt_count, trrpt_count_recv,            &
                 trsp_count, trsp_count_sum, trsp_count_recv,                   &
-                trsp_count_recv_sum, trspt_count, trspt_count_recv, nd
+                trsp_count_recv_sum, trspt_count, trspt_count_recv 
 
     INTEGER ::  gp_outside_of_building(1:8)
@@ -151,11 +155,12 @@
                 dt_particle, dt_particle_m, d_radius, d_sum, e_a, e_int,       &
                 e_int_l, e_int_u, e_mean_int, e_s, exp_arg, exp_term, fs_int,  &
-                gg,lagr_timescale, mean_r, 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, &
-                random_gauss, sl_r3, sl_r4, t_int, u_int, u_int_l, u_int_u,    &
-                vv_int, v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u,      &
-                x, y
+                gg, integral, lagr_timescale, lw_max, mean_r, 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, random_gauss, sl_r3, sl_r4, t_int, u_int,  &
+                u_int_l, u_int_u,vv_int, v_int, v_int_l, v_int_u, w_int,       &
+                w_int_l, w_int_u, x, y 
 
     REAL, DIMENSION(1:30) ::  de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
+    REAL, DIMENSION(:,:), ALLOCATABLE ::  kern
 
     REAL    ::  location(1:30,1:3)
@@ -399,4 +404,169 @@
        CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' )
 
+       IF ( wang_collision_kernel ) THEN
+
+       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
+                   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
+                   ENDDO
+
+                   psi = prt_start_index(k,j,i)
+                   pse = psi + prt_count(k,j,i)-1
+
+                   ALLOCATE( kern(psi:pse,psi:pse) )
+
+!
+!--                Calculate collision kernel for all particles in grid box
+                   CALL colker( i, j, k, kern )
+!
+!--                collison kernel is calculated in cm**3/s but needed in m**3/s
+                   kern = kern * 1.0E-6
+
+                   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 
+                      DO  is = psi, n-1
+
+                         integral = integral +                               &
+                                    ( particles(is)%radius**3 * kern(n,is) * &
+                                      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 )
+
+                      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
+
+                         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
+
+                         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 - &
+                               ( ( kern(n,is) * 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( 'advec_particles', '', 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
+
+                   ENDDO
+
+                   ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
+                                                 * particles(psi)%radius**3
+
+                   DEALLOCATE( kern )
+
+                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 *              &
+                                   particles(psi)%radius**3
+                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( 'advec_particles', '', 2, 2, -1, 6,  &
+                                    1 )
+                   ENDIF
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDDO
+
+       ELSE
+
        DO  i = nxl, nxr
           DO  j = nys, nyn
@@ -683,4 +853,6 @@
           ENDDO
        ENDDO
+
+       ENDIF
 
        CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' )
@@ -1200,5 +1372,5 @@
                    de_dzi(num_gp) = de_dz(k+1,j+1,i)
                 ENDIF
- 
+
                 IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
                 THEN
@@ -1620,5 +1792,5 @@
                       de_dzi(num_gp) = 0.0
                    ENDIF
- 
+
 !
 !--                If wall between gridpoint 5 and gridpoint 7, then
@@ -1971,5 +2143,5 @@
 #else
        dt_3d_reached = dt_3d_reached_l
-#endif     
+#endif
 
        CALL cpu_log( log_point_s(44), 'advec_part_advec', 'stop' )
@@ -2162,8 +2334,8 @@
 !--       data
           DO  n = 1, number_of_particles
-	     i = ( particles(n)%x + 0.5 * dx ) * ddx
+             i = ( particles(n)%x + 0.5 * dx ) * ddx
 !
 !--          Above calculation does not work for indices less than zero
-	     IF ( particles(n)%x < -0.5 * dx )  i = -1
+             IF ( particles(n)%x < -0.5 * dx )  i = -1
 
              IF ( i < nxl )  THEN
@@ -2225,8 +2397,8 @@
           nn = particles(n)%tail_id
 
-	  i = ( particles(n)%x + 0.5 * dx ) * ddx
+          i = ( particles(n)%x + 0.5 * dx ) * ddx
 !
 !--       Above calculation does not work for indices less than zero
-	  IF ( particles(n)%x < - 0.5 * dx )  i = -1
+          IF ( particles(n)%x < - 0.5 * dx )  i = -1
 
           IF ( i <  nxl )  THEN
@@ -2602,8 +2774,8 @@
           DO  n = 1, number_of_particles
              IF ( particle_mask(n) )  THEN
-	        j = ( particles(n)%y + 0.5 * dy ) * ddy
+                j = ( particles(n)%y + 0.5 * dy ) * ddy
 !
 !--             Above calculation does not work for indices less than zero
-	        IF ( particles(n)%y < -0.5 * dy )  j = -1
+                IF ( particles(n)%y < -0.5 * dy )  j = -1
 
                 IF ( j < nys )  THEN
@@ -2669,8 +2841,8 @@
 !--       moved.
           IF ( particle_mask(n) )  THEN
-	     j = ( particles(n)%y + 0.5 * dy ) * ddy
+             j = ( particles(n)%y + 0.5 * dy ) * ddy
 !
 !--          Above calculation does not work for indices less than zero
-	     IF ( particles(n)%y < -0.5 * dy )  j = -1
+             IF ( particles(n)%y < -0.5 * dy )  j = -1
 
              IF ( j < nys )  THEN
@@ -3815,5 +3987,4 @@
 
     CHARACTER (LEN=10) ::  particle_binary_version
-
     INTEGER ::  i
 
@@ -3879,4 +4050,5 @@
 
  END SUBROUTINE write_particles
+
 
 
Index: /palm/trunk/SOURCE/compute_vpt.f90
===================================================================
--- /palm/trunk/SOURCE/compute_vpt.f90	(revision 798)
+++ /palm/trunk/SOURCE/compute_vpt.f90	(revision 799)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: ql is now included in calculation of vpt in case of 
+!         cloud droplets 
 !
 ! Former revisions:
@@ -32,11 +33,13 @@
     INTEGER :: k
 
-    IF ( .NOT. cloud_physics ) THEN
+    IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
        vpt = pt * ( 1.0 + 0.61 * q )
-    ELSE
+    ELSE IF (cloud_physics) THEN
        DO  k = nzb, nzt+1
           vpt(k,:,:) = ( pt(k,:,:) + pt_d_t(k) * l_d_cp * ql(k,:,:) ) * &
                        ( 1.0 + 0.61 * q(k,:,:) - 1.61 * ql(k,:,:) ) 
        ENDDO
+    ELSE
+       vpt = pt * ( 1.0 + 0.61 * q - 1.61 * ql ) 
     ENDIF
 
Index: /palm/trunk/SOURCE/interaction_droplets_ptq.f90
===================================================================
--- /palm/trunk/SOURCE/interaction_droplets_ptq.f90	(revision 798)
+++ /palm/trunk/SOURCE/interaction_droplets_ptq.f90	(revision 799)
@@ -4,5 +4,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: pt_d_t(k) was missing in calculation of pt_p
 !
 ! Former revisions:
@@ -53,5 +53,5 @@
              DO  k = nzb_2d(j,i)+1, nzt
                 q_p(k,j,i)  = q_p(k,j,i)  - ql_c(k,j,i)
-                pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i)
+                pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i) * pt_d_t(k)
              ENDDO
           ENDDO
@@ -80,5 +80,5 @@
        DO  k = nzb_2d(j,i)+1, nzt
           q_p(k,j,i)  = q_p(k,j,i)  - ql_c(k,j,i)
-          pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i)
+          pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i) * pt_d_t(k)
        ENDDO
 
Index: /palm/trunk/SOURCE/wang_kernel.f90
===================================================================
--- /palm/trunk/SOURCE/wang_kernel.f90	(revision 798)
+++ /palm/trunk/SOURCE/wang_kernel.f90	(revision 799)
@@ -4,5 +4,7 @@
 ! Current revisions:
 ! -----------------
-! speed optimizations + formatting
+! speed optimizations and formatting
+! Bugfix: iq=1 is not allowed (routine effic) 
+! Bugfix: replaced stop by ec=0.0 in case of very small ec (routine effic)  
 !
 ! Former revisions:
@@ -169,4 +171,5 @@
 ! of droplets, air density, air viscosity, turbulent dissipation rate,
 ! taylor microscale reynolds number, gravitational acceleration
+
        USE constants
        USE cloud_parameters
@@ -176,5 +179,5 @@
        IMPLICIT NONE
 
-       REAL :: airvisc, airdens, waterdens, gravity, anu, Relamda, &
+       REAL :: Relamda, &
                Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
                bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
@@ -182,22 +185,32 @@
                wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
 
-       REAL, DIMENSION(pstart:pend) :: vST, tau, St
+       REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens
+
+       REAL, DIMENSION(pstart:pend) :: St, tau
+
+       LOGICAL, SAVE ::  first = .TRUE.
 
        INTEGER :: i, j
 
-       airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
-       airdens   = 1.2250   !air density in mg/cm**3
-       waterdens = 1000.0   !water density in mg/cm**3
-       gravity   = 980.6650 !in cm/s**2
-
-       anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
-
-       Relamda = urms**2.0*sqrt(15.0/epsilon/anu)
-
-       Tl = urms**2.0/epsilon       !in s
-       Lf = 0.5 * (urms**3.0)/epsilon !in cm
+!
+!--   Initial assignment of constants
+       IF ( first )  THEN
+
+          first = .FALSE.
+          airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
+          airdens   = 1.2250   !air density in mg/cm**3
+          waterdens = 1000.0   !water density in mg/cm**3
+          gravity   = 980.6650 !in cm/s**2
+          anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
+
+       ENDIF   
+
+       Relamda = urms**2*sqrt(15.0/epsilon/anu)
+
+       Tl = urms**2/epsilon       !in s
+       Lf = 0.5 * (urms**3)/epsilon !in cm
 
        tauk = (anu/epsilon)**0.5       !in s
-       eta  = (anu**3.0/epsilon)**0.25  !in cm
+       eta  = (anu**3/epsilon)**0.25  !in cm
        vk   = eta/tauk                   !in cm/s
 
@@ -211,6 +224,5 @@
 
        DO i = pstart, pend
-          vST(i) = winf(i)           !in cm/s
-          tau(i) = vST(i)/gravity    !in s
+          tau(i) = winf(i)/gravity    !in s
           St(i)  = tau(i)/tauk
        ENDDO
@@ -222,5 +234,5 @@
        be = sqrt(2.0)*lambda/Lf
 
-       bbb = sqrt(1.0-2.0*be**2.0)
+       bbb = sqrt(1.0-2.0*be**2)
        d1 = (1.+bbb)/2.0/bbb
        e1 = Lf*(1.0+bbb)/2.0   !in cm
@@ -228,5 +240,5 @@
        e2 = Lf*(1.0-bbb)/2.0   !in cm
 
-       ccc = sqrt(1.0-2.0*z**2.0)
+       ccc = sqrt(1.0-2.0*z**2)
        b1 = (1.+ccc)/2./ccc
        c1 = Tl*(1.+ccc)/2.   !in s
@@ -235,10 +247,10 @@
 
        DO i = pstart, pend
-          v1 = vST(i)         !in cm/s
+          v1 = winf(i)         !in cm/s
           t1 = tau(i)         !in s
 
           DO j = pstart,i
              rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm
-             v2  = vST(j)         !in cm/s
+             v2  = winf(j)         !in cm/s
              t2  = tau(j)         !in s
 
@@ -247,5 +259,5 @@
                      - b2*d1*PHI(c2,e1,v1,t1) &
                      + b2*d2*PHI(c2,e2,v1,t1)
-             v1xysq = v1xysq * urms**2.0/t1    !in cm**2/s**2
+             v1xysq = v1xysq * urms**2/t1    !in cm**2/s**2
 
              vrms1xy= sqrt(v1xysq)             !in cm/s
@@ -255,17 +267,17 @@
                      - b2*d1*PHI(c2,e1,v2,t2) &
                      + b2*d2*PHI(c2,e2,v2,t2)
-             v2xysq = v2xysq * urms**2.0/t2    !in cm**2/s**2
+             v2xysq = v2xysq * urms**2/t2    !in cm**2/s**2
 
               vrms2xy= sqrt(v2xysq)             !in cm/s
 
-             IF(vST(i).ge.vST(j)) THEN
-                v1 = vST(i)
+             IF(winf(i).ge.winf(j)) THEN
+                v1 = winf(i)
                 t1 = tau(i)
-                v2 = vST(j)
+                v2 = winf(j)
                 t2 = tau(j)
              ELSE
-                v1 = vST(j)
+                v1 = winf(j)
                 t1 = tau(j)
-                v2 = vST(i)
+                v2 = winf(i)
                 t2 = tau(i)
              ENDIF
@@ -276,11 +288,11 @@
                      + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
              fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
-             v1v2xy = v1v2xy * fR * urms**2.0/tau(i)/tau(j)    !in cm**2/s**2
-
-             wrtur2xy=vrms1xy**2.0 + vrms2xy**2.0 - 2.*v1v2xy  !in cm**2/s**2
+             v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j)    !in cm**2/s**2
+
+             wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy  !in cm**2/s**2
 
               IF (wrtur2xy.le.0.0) wrtur2xy=0.0
 
-              wrgrav2=pi/8.*(vST(j)-vST(i))**2.0
+              wrgrav2=pi/8.*(winf(j)-winf(i))**2
 
               wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
@@ -295,6 +307,6 @@
              ENDIF
 
-             XX = -0.1988*SSt**4.0 + 1.5275*SSt**3.0 &
-                  -4.2942*SSt**2.0 + 5.3406*SSt
+             XX = -0.1988*SSt**4 + 1.5275*SSt**3 &
+                  -4.2942*SSt**2 + 5.3406*SSt
 
              IF(XX.le.0.0) XX = 0.0
@@ -304,12 +316,12 @@
              c1_gr = XX/(gravity/(vk/tauk))**YY
 
-             ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2.0
+             ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2
              fao_gr = 20.115 * (ao_gr/Relamda)**0.5
              rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
 
-             grFIN = ((eta**2.0+rc**2.0)/(rrp**2.0+rc**2.0))**(c1_gr/2.)
+             grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.)
              IF (grFIN.lt.1.0) grFIN = 1.0
 
-             gck(i,j) = 2. * pi * rrp**2.0 * wrFIN * grFIN       ! in cm**3/s
+             gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN       ! in cm**3/s
              gck(j,i) = gck(i,j)
 
@@ -330,5 +342,5 @@
        aa1 = 1./tau0 + 1./a + vsett/b
 
-       PHI = 1./aa1 - vsett/2.0/b/aa1**2.0  !in s
+       PHI = 1./aa1 - vsett/2.0/b/aa1**2  !in s
 
     END FUNCTION PHI
@@ -346,18 +358,17 @@
         aa2 = vsett1/b + 1./tau1 + 1./a
         aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
-        aa4 = (vsett2/b)**2.0 - (1./tau2 + 1./a)**2.0
+        aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2
         aa5 = vsett2/b + 1./tau2 + 1./a
         aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2
-        ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2.0          &
-             + (4./aa4 - 1./aa5**2.0 - 1./aa1**2.0) * vsett2/2./b/aa6      &
-             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2.0 + vsett2/aa1**2.0)   &
+        ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2          &
+             + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6      &
+             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2)   &
                                                            * 1./2./b/aa3             ! in s**2
+
     END FUNCTION ZHI
 
-
 !------------------------------------------------------------------------------!
 ! SUBROUTINE for calculation of terminal velocity winf
 !------------------------------------------------------------------------------!
-
    SUBROUTINE fallg
 
@@ -371,32 +382,39 @@
       INTEGER :: i, j
 
-      REAL :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py,  &
-              x, y, xrey, bond
-
-      REAL, DIMENSION(1:7)  :: b
-      REAL, DIMENSION(1:6)  :: c
-      REAL, DIMENSION(1:20) :: rat
-      REAL, DIMENSION(1:15) :: r0
-      REAL, DIMENSION(1:15,1:20) :: ecoll
-
-      b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
+      LOGICAL, SAVE ::  first = .TRUE.
+
+      REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py
+
+      REAL :: bond, x, xrey, y
+
+      REAL, DIMENSION(1:7), SAVE  :: b
+      REAL, DIMENSION(1:6), SAVE  :: c
+
+!
+!--   Initial assignment of constants
+      IF ( first )  THEN
+
+         first = .FALSE.
+         b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
              -0.578878e-3,0.855176e-4,-0.327815e-5/)
-      c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
+         c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
              -0.542819e-1, 0.238449e-2/)
 
-      eta   = 1.818e-4     !in poise = g/(cm s)
-      xlamb = 6.62e-6    !in cm
-
-      rhow = 1.0       !in g/cm**3
-      rhoa = 1.225e-3    !in g/cm**3
-
-      grav = 980.665     !in cm/s**2
-      cunh = 1.257 * xlamb !in cm
-      t0   = 273.15        !in K
-      sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
-      stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
-      stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
-      phy  = (sigma**3.0) * (rhoa**2.0)/((eta**4.0) * grav * (rhow-rhoa))
-      py   = phy**(1.0/6.0)
+         eta   = 1.818e-4     !in poise = g/(cm s)
+         xlamb = 6.62e-6    !in cm
+
+         rhow = 1.0       !in g/cm**3
+         rhoa = 1.225e-3    !in g/cm**3
+
+         grav = 980.665     !in cm/s**2
+         cunh = 1.257 * xlamb !in cm
+         t0   = 273.15        !in K
+         sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
+         stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
+         stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
+         phy  = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa))
+         py   = phy**(1.0/6.0)
+
+      ENDIF
 
 !particle radius has to be in cm
@@ -405,9 +423,9 @@
          IF (particles(j)%radius*100.0 .le. 1.e-3) THEN
 
-            winf(j)=stok*((particles(j)%radius*100.0)**2.+cunh* particles(j)%radius*100.0) !in cm/s
+            winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s
 
          ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN
 
-            x = log(stb*(particles(j)%radius*100.0)**3.0)
+            x = log(stb*(particles(j)%radius*100.0)**3)
             y = 0.0
 
@@ -421,8 +439,8 @@
          ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
 
-            bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2.0/sigma
-
             IF (particles(j)%radius*100.0.gt.0.35)  THEN
                bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma
+            ELSE
+               bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma
             ENDIF
 
@@ -435,13 +453,14 @@
 
             xrey = py*exp(y)
-            winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
 
             IF (particles(j)%radius*100.0 .gt.0.35)  THEN
               winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
+            ELSE
+               winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
             ENDIF
 
          ENDIF
       ENDDO
-
+      RETURN
    END SUBROUTINE fallg
 
@@ -457,4 +476,5 @@
       USE particle_attributes
 
+!collision efficiencies of hall kernel
       IMPLICIT NONE
 
@@ -545,23 +565,17 @@
          DO  i = pstart, j
 
-!            DO k = 2, 15
-!               IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN
-!                  ir = k
-!               ELSEIF ((particles(j)%radius*1.0E06).gt.r0(15)) THEN
-!                  ir = 16
-!               ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN
-!                  ir = 1
+            ir = ira(j)
+
+            rq = particles(i)%radius / particles(j)%radius
+
+!            DO kk = 2, 21
+!               IF ( rq <= rat(kk) )  THEN
+!                  iq = kk
+!                  EXIT
 !               ENDIF
 !            ENDDO
 
-            ir = ira(j)
-
-            rq = particles(i)%radius / particles(j)%radius
-
-!            DO kk = 2, 21
-!               IF ( rq .le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
-!            ENDDO
-
             iq = INT( rq * 20 ) + 1
+            iq = MAX(iq , 2)
 
             IF ( ir < 16 )  THEN
@@ -587,5 +601,5 @@
 
             ec(i,j) = ec(j,i)
-            IF ( ec(i,j) < 1.0E-20 )  STOP 99
+            IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
 
          ENDDO
@@ -599,5 +613,4 @@
 ! SUBROUTINE for calculation of enhancement factor collision efficencies
 !------------------------------------------------------------------------------!
-
    SUBROUTINE turb_enhance_eff
 
@@ -607,42 +620,69 @@
        USE arrays_3d
 
-!collision efficiencies of hall kernel
       IMPLICIT NONE
 
       INTEGER :: i, ik, ir, iq, j, k, kk
-      REAL    :: rq, y1, pp, qq, y2, y3, x1, x2, x3
-
-      REAL, DIMENSION(1:11) :: rat
-      REAL, DIMENSION(1:7)  :: r0
-      REAL, DIMENSION(1:7,1:11) :: ecoll_100, ecoll_400
-
-      r0  = (/10., 20., 30.,40., 50., 60.,100./)
-      rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
+
+      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
+
+      REAL    :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3
+
+      LOGICAL, SAVE ::  first = .TRUE.
+
+      REAL, DIMENSION(1:11), SAVE :: rat
+      REAL, DIMENSION(1:7), SAVE  :: r0
+      REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400
+
+!
+!--   Initial assignment of constants
+      IF ( first )  THEN
+
+         first = .FALSE.
+
+         r0  = (/10., 20., 30.,40., 50., 60.,100./)
+         rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
 
 ! 100 cm^2/s^3
-      ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
-      ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
-      ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
-      ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
-      ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
-      ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
-      ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
-      ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
-      ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
-      ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
-      ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
+         ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
+         ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
+         ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
+         ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
+         ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
+         ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
+         ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
+         ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
+         ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
+         ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
+         ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
 
 ! 400 cm^2/s^3
-      ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
-      ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
-      ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
-      ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
-      ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
-      ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
-      ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
-      ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
-      ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
-      ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
-      ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
+         ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
+         ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
+         ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
+         ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
+         ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
+         ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
+         ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
+         ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
+         ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
+         ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
+         ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
+
+      ENDIF
+
+!
+!--   Calculate the radius class index of particles with respect to array r
+      ALLOCATE( ira(pstart:pend) )
+
+      DO  j = pstart, pend
+         particle_radius = particles(j)%radius * 1.0E6
+         DO  k = 1, 7
+            IF ( particle_radius < r0(k) )  THEN
+               ira(j) = k
+               EXIT
+            ENDIF
+         ENDDO
+         IF ( particle_radius >= r0(7) )  ira(j) = 8
+      ENDDO
 
 ! two-dimensional linear interpolation of the collision efficiency
@@ -650,18 +690,13 @@
          DO i = pstart, j
 
-            DO k = 2, 7
-               IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN
-                  ir = k
-               ELSEIF ((particles(j)%radius*1.0E06).gt.r0(7)) THEN
-                  ir = 8
-               ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN
-                  ir = 1
+            ir = ira(j)
+
+            rq = particles(i)%radius/particles(j)%radius
+
+            DO kk = 2, 11
+               IF ( rq <= rat(kk) )  THEN
+                  iq = kk
+                  EXIT
                ENDIF
-            ENDDO
-
-            rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06)
-
-            DO kk = 2, 11
-               IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
             ENDDO
 
