Index: palm/trunk/SOURCE/wang_kernel.f90
===================================================================
--- palm/trunk/SOURCE/wang_kernel.f90	(revision 792)
+++ 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
 
