Index: palm/trunk/SOURCE/lpm_collision_kernels.f90
===================================================================
--- palm/trunk/SOURCE/lpm_collision_kernels.f90	(revision 826)
+++ palm/trunk/SOURCE/lpm_collision_kernels.f90	(revision 828)
@@ -1,8 +1,13 @@
-MODULE lpm_collision_kernels_mod
+ MODULE lpm_collision_kernels_mod
 
 !------------------------------------------------------------------------------!
 ! Current revisions:
 ! -----------------
-! 
+! code has been completely reformatted, routine colker renamed
+! recalculate_kernel,
+! routine init_kernels added, radius is now communicated to the collision
+! routines by array radclass
+!
+! Bugfix: transformation factor for dissipation changed from 1E5 to 1E4
 !
 ! Former revisions:
@@ -24,7 +29,15 @@
 ! Description:
 ! ------------
-! This routine calculates the effect of (SGS) turbulence on the collision
-! efficiency of droplets.
-! It is based on the original kernel developed by Wang (...)
+! This module calculates collision efficiencies either due to pure gravitational
+! effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or
+! including the effects of (SGS) turbulence (Wang kernel, see Wang and
+! Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8). The original code has been
+! provided by L.-P. Wang but is substantially reformatted and speed optimized
+! here.
+!
+! ATTENTION:
+! Physical quantities (like g, densities, etc.) used in this module still
+! have to be adjusted to those values used in the main PALM code.
+! Also, quantities in CGS-units should be converted to SI-units eventually.
 !------------------------------------------------------------------------------!
 
@@ -33,4 +46,6 @@
     USE constants
     USE particle_attributes
+    USE pegrid
+
 
     IMPLICIT NONE
@@ -38,49 +53,164 @@
     PRIVATE
 
-    PUBLIC  colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi 
-
-    INTEGER, SAVE ::  ip, jp, kp, pend, pstart, psum
-
-    REAL, SAVE :: epsilon, eps2, urms, urms2
-
-    REAL, DIMENSION(:),   ALLOCATABLE, SAVE ::  winf
-    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  ec, ecf, gck
+    PUBLIC  ckernel, init_kernels,  rclass_lbound, rclass_ubound, &
+            recalculate_kernel            
+
+    REAL ::  epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2
+
+    REAL, DIMENSION(:),   ALLOCATABLE ::  epsclass, radclass, winf
+    REAL, DIMENSION(:,:), ALLOCATABLE ::  ec, ecf, gck, hkernel, hwratio
+    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  ckernel
+
+    SAVE
 
 !
 !-- Public interfaces
-    INTERFACE colker
-       MODULE PROCEDURE colker
-    END INTERFACE colker
-
-    INTERFACE effic
-       MODULE PROCEDURE effic
-    END INTERFACE effic
-
-    INTERFACE fallg
-       MODULE PROCEDURE fallg
-    END INTERFACE fallg
-
-    INTERFACE phi
-       MODULE PROCEDURE phi
-    END INTERFACE phi
-
-    INTERFACE turbsd
-       MODULE PROCEDURE turbsd
-    END INTERFACE turbsd
-
-    INTERFACE turb_enhance_eff
-       MODULE PROCEDURE turb_enhance_eff
-    END INTERFACE turb_enhance_eff
-
-    INTERFACE zhi
-       MODULE PROCEDURE zhi
-    END INTERFACE zhi
+    INTERFACE init_kernels
+       MODULE PROCEDURE init_kernels
+    END INTERFACE init_kernels
+
+    INTERFACE recalculate_kernel
+       MODULE PROCEDURE recalculate_kernel
+    END INTERFACE recalculate_kernel
+
 
     CONTAINS
 
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of collision kernel
-!------------------------------------------------------------------------------!
-    SUBROUTINE colker( i1, j1, k1, kernel )
+
+    SUBROUTINE init_kernels
+!------------------------------------------------------------------------------!
+! Initialization of the collision efficiency matrix with fixed radius and
+! dissipation classes, calculated at simulation start only.
+!------------------------------------------------------------------------------!
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, j, k
+
+
+!
+!--    Calculate collision efficiencies for fixed radius- and dissipation
+!--    classes
+       IF ( collision_kernel(6:9) == 'fast' )  THEN
+
+          ALLOCATE( ckernel(1:radius_classes,1:radius_classes,               &
+                    0:dissipation_classes), epsclass(1:dissipation_classes), &
+                    radclass(1:radius_classes) )
+
+!
+!--       Calculate the radius class bounds with logarithmic distances
+!--       in the interval [1.0E-6, 2.0E-4] m
+          rclass_lbound = LOG( 1.0E-6 )
+          rclass_ubound = LOG( 2.0E-4 )
+          radclass(1)   = 1.0E-6
+          DO  i = 2, radius_classes
+             radclass(i) = EXP( rclass_lbound +                                &
+                                ( rclass_ubound - rclass_lbound ) * ( i-1.0 ) /&
+                                ( radius_classes - 1.0 ) )
+!             IF ( myid == 0 )  THEN
+!                PRINT*, 'i=', i, ' r = ', radclass(i)*1.0E6
+!             ENDIF
+          ENDDO
+!
+!--       Collision routines expect radius to be in cm
+          radclass = radclass * 100.0
+
+!
+!--       Set the class bounds for dissipation in interval [0, 1000] cm**2/s**3
+          DO  i = 1, dissipation_classes
+             epsclass(i) = 1000.0 * REAL( i ) / dissipation_classes
+!             IF ( myid == 0 )  THEN
+!                PRINT*, 'i=', i, ' eps = ', epsclass(i)
+!             ENDIF
+          ENDDO
+!
+!--       Calculate collision efficiencies of the Wang/ayala kernel
+          ALLOCATE( ec(1:radius_classes,1:radius_classes),  &
+                    ecf(1:radius_classes,1:radius_classes), &
+                    gck(1:radius_classes,1:radius_classes), &
+                    winf(1:radius_classes) )
+
+          DO  k = 1, dissipation_classes
+
+             epsilon = epsclass(k)
+             urms    = 202.0 * ( epsilon / 400.0 )**( 1.0 / 3.0 )
+
+             CALL turbsd
+             CALL turb_enhance_eff
+             CALL effic
+
+             DO  j = 1, radius_classes
+                DO  i = 1, radius_classes
+                   ckernel(i,j,k) = ec(i,j) * gck(i,j) * ecf(i,j)
+                ENDDO
+             ENDDO
+
+          ENDDO
+
+!
+!--       Calculate collision efficiencies of the Hall kernel
+          ALLOCATE( hkernel(1:radius_classes,1:radius_classes), &
+                    hwratio(1:radius_classes,1:radius_classes) )
+
+          CALL fallg
+          CALL effic
+
+          DO  j = 1, radius_classes
+             DO  i =  1, radius_classes
+                hkernel(i,j) = pi * ( radclass(j) + radclass(i) )**2 &
+                                  * ec(i,j) * ABS( winf(j) - winf(i) )
+                ckernel(i,j,0) = hkernel(i,j)  ! hall kernel stored on index 0
+              ENDDO
+          ENDDO
+
+!
+!--       Test output of efficiencies
+          IF ( j == -1 )  THEN
+
+             PRINT*, '*** Hall kernel'
+             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
+             DO  j = 1, radius_classes
+                WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes )
+             ENDDO
+
+             DO  k = 1, dissipation_classes
+                DO  i = 1, radius_classes
+                   DO  j = 1, radius_classes
+                      IF ( hkernel(i,j) == 0.0 )  THEN
+                         hwratio(i,j) = 9999999.9
+                      ELSE
+                         hwratio(i,j) = ckernel(i,j,k) / hkernel(i,j)
+                      ENDIF
+                   ENDDO
+                ENDDO
+
+                PRINT*, '*** epsilon = ', epsclass(k)
+                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
+                DO  j = 1, radius_classes
+!                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( ckernel(i,j,k), i = 1,radius_classes )
+                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( hwratio(i,j), i = 1,radius_classes )
+                ENDDO
+             ENDDO
+
+          ENDIF
+
+          DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf )
+
+          ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
+
+       ELSEIF( collision_kernel == 'hall'  .OR.  collision_kernel == 'wang' ) &
+       THEN
+!
+!--       Initial settings for Hall- and Wang-Kernel
+!--       To be done: move here parts from turbsd, fallg, ecoll, etc.
+       ENDIF
+
+    END SUBROUTINE init_kernels
+
+
+!------------------------------------------------------------------------------!
+! Calculation of collision kernels during each timestep and for each grid box
+!------------------------------------------------------------------------------!
+    SUBROUTINE recalculate_kernel( i1, j1, k1 )
 
        USE arrays_3d
@@ -94,88 +224,71 @@
        IMPLICIT NONE
 
-       INTEGER ::  i, i1, j, j1, k1
-
-       REAL, DIMENSION(prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1,prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1) :: kernel
-
-!       CALL cpu_log( log_point_s(46), 'colker', 'start' )
-
-       ip = i1
-       jp = j1
-       kp = k1
-
-       pstart = prt_start_index(kp,jp,ip)
-       pend   = prt_start_index(kp,jp,ip) + prt_count(kp,jp,ip) - 1
-       psum   = prt_count(kp,jp,ip)
-
-       ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
-
-       IF ( wang_kernel )  THEN
-
-          ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
-
-          epsilon = diss(kp,jp,ip) * 1.0E5  !dissipation rate in cm**2/s**-3
-          urms     = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 )
-
-          IF ( epsilon <= 0.001 )  THEN
-
-             CALL fallg
-             CALL effic
-
-             DO  j = pstart, pend
-                DO  i =  pstart, pend
-                   kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
-                                        particles(i)%radius * 100.0 )**2 &
-                                    * ec(i,j) * ABS( winf(j) - winf(i) )
-                ENDDO
+       INTEGER ::  i, i1, j, j1, k1, pend, pstart
+
+
+       pstart = prt_start_index(k1,j1,i1)
+       pend   = prt_start_index(k1,j1,i1) + prt_count(k1,j1,i1) - 1
+       radius_classes = prt_count(k1,j1,i1)
+
+       ALLOCATE( ec(1:radius_classes,1:radius_classes), &
+                 radclass(1:radius_classes), winf(1:radius_classes) )
+
+!
+!--    Store particle radii on the radclass array. Collision routines
+!--    expect radii to be in cm.
+       radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0
+
+       epsilon = diss(k1,j1,i1) * 1.0E4   ! dissipation rate in cm**2/s**-3
+       urms    = 202.0 * ( epsilon / 400.0 )**( 0.33333333333 )
+
+       IF ( wang_kernel  .AND.  epsilon > 0.001 )  THEN
+!
+!--       Call routines to calculate efficiencies for the Wang kernel
+          ALLOCATE( gck(1:radius_classes,1:radius_classes), &
+                    ecf(1:radius_classes,1:radius_classes) )
+
+          CALL turbsd
+          CALL turb_enhance_eff
+          CALL effic
+
+          DO  j = 1, radius_classes
+             DO  i =  1, radius_classes
+                ckernel(pstart+i-1,pstart+j-1,1) = ec(i,j) * gck(i,j) * ecf(i,j)
              ENDDO
-
-          ELSE
-
-             CALL turbsd
-             CALL turb_enhance_eff
-             CALL effic
-
-             DO  j = pstart, pend, 1
-                DO  i =  pstart, pend, 1
-                   kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
-                ENDDO
+          ENDDO
+
+          DEALLOCATE( gck, ecf )
+
+       ELSE
+!
+!--       Call routines to calculate efficiencies for the Hall kernel
+          CALL fallg
+          CALL effic
+
+          DO  j = 1, radius_classes
+             DO  i =  1, radius_classes
+                ckernel(pstart+i-1,pstart+j-1,1) = pi *                       &
+                                          ( radclass(j) + radclass(i) )**2    &
+                                          * ec(i,j) * ABS( winf(j) - winf(i) )
              ENDDO
-
-          ENDIF
-
-          DEALLOCATE(gck, ecf)
-
-       ELSE
-
-!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )
-          CALL fallg
-!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )
-!          CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )
-          CALL effic
-!          CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )
-
-          DO  j = pstart, pend
-             DO  i =  pstart, pend
-                kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
-                                     particles(i)%radius * 100.0 )**2 &
-                                 * ec(i,j) * ABS( winf(j) - winf(i) )
-             ENDDO
           ENDDO
 
        ENDIF
 
-       DEALLOCATE( ec, winf )
-
-!       CALL cpu_log( log_point_s(46), 'colker', 'stop' )
-
-    END SUBROUTINE colker
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of w, g and gck
+       ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
+
+       DEALLOCATE( ec, radclass, winf )
+
+    END SUBROUTINE recalculate_kernel
+
+
+!------------------------------------------------------------------------------!
+! Calculation of gck
+! This is from Aayala 2008b, page 37ff.
+! Necessary input parameters: water density, radii of droplets, air density,
+! air viscosity, turbulent dissipation rate, taylor microscale reynolds number,
+! gravitational acceleration  --> to be replaced by PALM parameters
 !------------------------------------------------------------------------------!
     SUBROUTINE turbsd
-! from Aayala 2008b, page 37ff, necessary input parameter water density, radii
-! of droplets, air density, air viscosity, turbulent dissipation rate,
-! taylor microscale reynolds number, gravitational acceleration
 
        USE constants
@@ -186,97 +299,85 @@
        IMPLICIT NONE
 
-       REAL :: Relamda, &
-               Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
-               bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
-               vrms1xy, v2xysq, vrms2xy, v1v2xy, fR, wrtur2xy, wrgrav2, &
-               wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
-
-       REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens
-
-       REAL, DIMENSION(pstart:pend) :: St, tau
+       INTEGER ::  i, j
 
        LOGICAL, SAVE ::  first = .TRUE.
 
-       INTEGER :: i, j
-
-!
-!--   Initial assignment of constants
+       REAL ::  ao, ao_gr, bbb, be, b1, b2, ccc, c1, c1_gr, c2, d1, d2, eta, &
+                e1, e2, fao_gr, fr, grfin, lambda, lambda_re, lf, rc, rrp,   &
+                sst, tauk, tl, t2, tt, t1, vk, vrms1xy, vrms2xy, v1, v1v2xy, &
+                v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z
+
+       REAL, SAVE ::  airdens, airvisc, anu, gravity, waterdens
+
+       REAL, DIMENSION(1:radius_classes) ::  st, tau
+
+
+!
+!--    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
+          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/epsilon)**0.25  !in cm
-       vk   = eta/tauk                   !in cm/s
-
-       ao = (11.+7.*Relamda)/(205.+Relamda)
-
-       lambda = urms * sqrt(15.*anu/epsilon)     !in cm
-
-       tt = sqrt(2.*Relamda/(15.**0.5)/ao) * tauk !in s
-
-       CALL fallg !gives winf in cm/s
-
-       DO i = pstart, pend
-          tau(i) = winf(i)/gravity    !in s
-          St(i)  = tau(i)/tauk
+       lambda    = urms * SQRT( 15.0 * anu / epsilon )     ! in cm
+       lambda_re = urms**2 * SQRT( 15.0 / epsilon / anu )
+       tl        = urms**2 / epsilon                       ! in s
+       lf        = 0.5 * urms**3 / epsilon                 ! in cm
+       tauk      = SQRT( anu / epsilon )                   ! in s
+       eta       = ( anu**3 / epsilon )**0.25              ! in cm
+       vk        = eta / tauk                              ! in cm/s
+
+       ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
+       tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk   ! in s
+
+       CALL fallg    ! gives winf in cm/s
+
+       DO  i = 1, radius_classes
+          tau(i) = winf(i) / gravity    ! in s
+          st(i)  = tau(i) / tauk
        ENDDO
 
-!***** TO CALCULATE wr ********************************
-!from Aayala 2008b, page 38f
-
-       z  = tt/Tl
-       be = sqrt(2.0)*lambda/Lf
-
-       bbb = sqrt(1.0-2.0*be**2)
-       d1 = (1.+bbb)/2.0/bbb
-       e1 = Lf*(1.0+bbb)/2.0   !in cm
-       d2 = (1.0-bbb)/2.0/bbb
-       e2 = Lf*(1.0-bbb)/2.0   !in cm
-
-       ccc = sqrt(1.0-2.0*z**2)
-       b1 = (1.+ccc)/2./ccc
-       c1 = Tl*(1.+ccc)/2.   !in s
-       b2 = (1.-ccc)/2./ccc
-       c2 = Tl*(1.-ccc)/2.   !in s
-
-       DO i = pstart, pend
-          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  = winf(j)         !in cm/s
-             t2  = tau(j)         !in s
-
-             v1xysq =  b1*d1*PHI(c1,e1,v1,t1) &
-                     - b1*d2*PHI(c1,e2,v1,t1) &
-                     - b2*d1*PHI(c2,e1,v1,t1) &
-                     + b2*d2*PHI(c2,e2,v1,t1)
-             v1xysq = v1xysq * urms**2/t1    !in cm**2/s**2
-
-             vrms1xy= sqrt(v1xysq)             !in cm/s
-
-             v2xysq =  b1*d1*PHI(c1,e1,v2,t2) &
-                     - b1*d2*PHI(c1,e2,v2,t2) &
-                     - b2*d1*PHI(c2,e1,v2,t2) &
-                     + b2*d2*PHI(c2,e2,v2,t2)
-             v2xysq = v2xysq * urms**2/t2    !in cm**2/s**2
-
-              vrms2xy= sqrt(v2xysq)             !in cm/s
-
-             IF(winf(i).ge.winf(j)) THEN
+!
+!--    Calculate wr (from Aayala 2008b, page 38f)
+       z   = tt / tl
+       be  = SQRT( 2.0 ) * lambda / lf
+       bbb = SQRT( 1.0 - 2.0 * be**2 )
+       d1  = ( 1.0 + bbb ) / ( 2.0 * bbb )
+       e1  = lf * ( 1.0 + bbb ) * 0.5   ! in cm
+       d2  = ( 1.0 - bbb ) * 0.5 / bbb
+       e2  = lf * ( 1.0 - bbb ) * 0.5   ! in cm
+       ccc = SQRT( 1.0 - 2.0 * z**2 )
+       b1  = ( 1.0 + ccc ) * 0.5 / ccc
+       c1  = tl * ( 1.0 + ccc ) * 0.5   ! in s
+       b2  = ( 1.0 - ccc ) * 0.5 / ccc
+       c2  = tl * ( 1.0 - ccc ) * 0.5   ! in s
+
+       DO  i = 1, radius_classes
+
+          v1 = winf(i)        ! in cm/s
+          t1 = tau(i)         ! in s
+
+          DO  j = 1, i
+             rrp = radclass(i) + radclass(j)               ! radius in cm
+             v2  = winf(j)                                 ! in cm/s
+             t2  = tau(j)                                  ! in s
+
+             v1xysq  = b1 * d1 * phi(c1,e1,v1,t1) - b1 * d2 * phi(c1,e2,v1,t1) &
+                     - b2 * d1 * phi(c2,e1,v1,t1) + b2 * d2 * phi(c2,e2,v1,t1)
+             v1xysq  = v1xysq * urms**2 / t1                ! in cm**2/s**2
+             vrms1xy = SQRT( v1xysq )                       ! in cm/s
+
+             v2xysq  = b1 * d1 * phi(c1,e1,v2,t2) - b1 * d2 * phi(c1,e2,v2,t2) &
+                     - b2 * d1 * phi(c2,e1,v2,t2) + b2 * d2 * phi(c2,e2,v2,t2)
+             v2xysq  = v2xysq * urms**2 / t2                ! in cm**2/s**2
+             vrms2xy = SQRT( v2xysq )                       ! in cm/s
+
+             IF ( winf(i) >= winf(j) )  THEN
                 v1 = winf(i)
                 t1 = tau(i)
@@ -290,45 +391,38 @@
              ENDIF
 
-             v1v2xy =  b1*d1*ZHI(c1,e1,v1,t1,v2,t2) &
-                     - b1*d2*ZHI(c1,e2,v1,t1,v2,t2) &
-                     - b2*d1*ZHI(c2,e1,v1,t1,v2,t2) &
-                     + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
-             fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
-             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.*(winf(j)-winf(i))**2
-
-              wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
-
-
-!***** TO CALCULATE gr ********************************
-
-             IF(St(j).gt.St(i)) THEN
-                SSt = St(j)
+             v1v2xy   =  b1 * d1 * zhi(c1,e1,v1,t1,v2,t2) - &
+                         b1 * d2 * zhi(c1,e2,v1,t1,v2,t2) - &
+                         b2 * d1 * zhi(c2,e1,v1,t1,v2,t2) + &
+                         b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
+             fr       = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
+             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)  ! in cm**2/s**2
+             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in cm**2/s**2
+             IF ( wrtur2xy < 0.0 )  wrtur2xy = 0.0
+             wrgrav2  = pi / 8.0 * ( winf(j) - winf(i) )**2
+             wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in cm/s
+
+!
+!--          Calculate gr
+             IF ( st(j) > st(i) )  THEN
+                sst = st(j)
              ELSE
-                SSt = St(i)
+                sst = st(i)
              ENDIF
 
-             XX = -0.1988*SSt**4 + 1.5275*SSt**3 &
-                  -4.2942*SSt**2 + 5.3406*SSt
-
-             IF(XX.le.0.0) XX = 0.0
-
-             YY = 0.1886*exp(20.306/Relamda)
-
-             c1_gr = XX/(gravity/(vk/tauk))**YY
-
-             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+rc**2)/(rrp**2+rc**2))**(c1_gr/2.)
-             IF (grFIN.lt.1.0) grFIN = 1.0
-
-             gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN       ! in cm**3/s
+             xx = -0.1988 * sst**4 + 1.5275 * sst**3 - 4.2942 * sst**2 + &
+                   5.3406 * sst
+             IF ( xx < 0.0 )  xx = 0.0
+             yy = 0.1886 * EXP( 20.306 / lambda_re )
+
+             c1_gr  =  xx / ( gravity / vk * tauk )**yy
+
+             ao_gr  = ao + ( pi / 8.0) * ( gravity / vk * tauk )**2
+             fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
+             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
+
+             grfin  = ( ( eta**2 + rc**2 ) / ( rrp**2 + rc**2) )**( c1_gr*0.5 )
+             IF ( grfin < 1.0 )  grfin = 1.0
+
+             gck(i,j) = 2.0 * pi * rrp**2 * wrfin * grfin           ! in cm**3/s
              gck(j,i) = gck(i,j)
 
@@ -336,47 +430,49 @@
        ENDDO
 
-    END SUBROUTINE TurbSD
-
-!------------------------------------------------------------------------------!
-! PHI as a function
-!------------------------------------------------------------------------------!
-  REAL FUNCTION PHI(a,b,vsett,tau0)
+    END SUBROUTINE turbsd
+
+
+!------------------------------------------------------------------------------!
+! phi as a function
+!------------------------------------------------------------------------------!
+    REAL FUNCTION phi( a, b, vsett, tau0 )
 
        IMPLICIT NONE
 
-       REAL :: a, aa1, b, vsett, tau0
-
-       aa1 = 1./tau0 + 1./a + vsett/b
-
-       PHI = 1./aa1 - vsett/2.0/b/aa1**2  !in s
-
-    END FUNCTION PHI
-
-!------------------------------------------------------------------------------!
-! ZETA as a function
-!------------------------------------------------------------------------------!
-  REAL FUNCTION ZHI(a,b,vsett1,tau1,vsett2,tau2)
+       REAL ::  a, aa1, b, tau0, vsett
+
+       aa1 = 1.0 / tau0 + 1.0 / a + vsett / b
+       phi = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
+
+    END FUNCTION phi
+
+
+!------------------------------------------------------------------------------!
+! zeta as a function
+!------------------------------------------------------------------------------!
+    REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
 
        IMPLICIT NONE
 
-       REAL :: a, aa1, aa2, aa3, aa4, aa5, aa6, b, vsett1, tau1, vsett2, tau2
-
-        aa1 = vsett2/b - 1./tau2 - 1./a
-        aa2 = vsett1/b + 1./tau1 + 1./a
-        aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
-        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          &
-             + (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
+       REAL ::  a, aa1, aa2, aa3, aa4, aa5, aa6, b, tau1, tau2, vsett1, vsett2
+
+       aa1 = vsett2 / b - 1.0 / tau2 - 1.0 / a
+       aa2 = vsett1 / b + 1.0 / tau1 + 1.0 / a
+       aa3 = ( vsett1 - vsett2 ) / b + 1.0 / tau1 + 1.0 / tau2
+       aa4 = ( vsett2 / b )**2 - ( 1.0 / tau2 + 1.0 / a )**2
+       aa5 = vsett2 / b + 1.0 / tau2 + 1.0 / a
+       aa6 = 1.0 / tau1 - 1.0 / a + ( 1.0 / tau2 + 1.0 / a) * vsett1 / vsett2
+       zhi = (1.0 / aa1 - 1.0 / aa2 ) * ( vsett1 - vsett2 ) * 0.5 / b / aa3**2 &
+           + (4.0 / aa4 - 1.0 / aa5**2 - 1.0 / aa1**2 ) * vsett2 * 0.5 / b /aa6&
+           + (2.0 * ( b / aa2 - b / aa1 ) - vsett1 / aa2**2 + vsett2 / aa1**2 )&
+           * 0.5 / b / aa3      ! in s**2
+
+    END FUNCTION zhi
+
+
+!------------------------------------------------------------------------------!
+! Calculation of terminal velocity winf
+!------------------------------------------------------------------------------!
+    SUBROUTINE fallg
 
        USE constants
@@ -385,240 +481,230 @@
        USE arrays_3d
 
-      IMPLICIT NONE
-
-      INTEGER :: i, j
-
-      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
+       IMPLICIT NONE
+
+       INTEGER ::  i, j
+
+       LOGICAL, SAVE ::  first = .TRUE.
+
+       REAL, SAVE ::  cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, &
+                      t0, xlamb
+
+       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, &
+                 -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 * rhoa**2 / ( eta**4 * grav * ( rhow - rhoa ) )
+          py    = phy**( 1.0 / 6.0 )
+
+       ENDIF
+
+       DO  j = 1, radius_classes
+
+          IF ( radclass(j) <= 1.0E-3 ) THEN
+
+             winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) ! in cm/s
+
+          ELSEIF ( radclass(j) > 1.0E-3  .AND.  radclass(j) <= 5.35E-2 )  THEN
+
+             x = LOG( stb * radclass(j)**3 )
+             y = 0.0
+
+             DO  i = 1, 7
+                y = y + b(i) * x**(i-1)
+             ENDDO
+
+             xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
+             winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )      ! in cm/s
+
+          ELSEIF ( radclass(j) > 5.35E-2 )  THEN
+
+             IF ( radclass(j) > 0.35 )  THEN
+                bond = grav * ( rhow - rhoa ) * 0.35**2 / sigma
+             ELSE
+                bond = grav * ( rhow - rhoa ) * radclass(j)**2 / sigma
+             ENDIF
+
+             x = LOG( 16.0 * bond * py / 3.0 )
+             y = 0.0
+
+             DO  i = 1, 6
+                y = y + c(i) * x**(i-1)
+             ENDDO
+
+             xrey = py * EXP( y )
+
+             IF ( radclass(j) > 0.35 )  THEN
+                winf(j) = xrey * eta / ( 2.0 * rhoa * 0.35 )           ! in cm/s
+             ELSE
+                winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )    ! in cm/s
+             ENDIF
+
+          ENDIF
+
+       ENDDO
+
+    END SUBROUTINE fallg
+
+
+!------------------------------------------------------------------------------!
+! Calculation of collision efficencies for the Hall kernel
+!------------------------------------------------------------------------------!
+    SUBROUTINE effic
+
+       USE arrays_3d
+       USE cloud_parameters
+       USE constants
+       USE particle_attributes
+
+       IMPLICIT NONE
+
+       INTEGER ::  i, iq, ir, j, k, kk
+
+       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
+
+       LOGICAL, SAVE ::  first = .TRUE.
+
+       REAL ::  ek, particle_radius, pp, qq, rq
+
+       REAL, DIMENSION(1:21), SAVE ::  rat
+       REAL, DIMENSION(1:15), SAVE ::  r0
+       REAL, DIMENSION(1:15,1:21), SAVE ::  ecoll
+
+!
+!--    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,    &
-             -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) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa))
-         py   = phy**(1.0/6.0)
-
-      ENDIF
-
-!particle radius has to be in cm
-      DO j = pstart, pend
-
-         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
-
-         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)
-            y = 0.0
-
-            DO i = 1, 7
-               y = y + b(i) * (x**(i-1))
-            ENDDO
-
-            xrey = (1.0 + cunh/(particles(j)%radius*100.0)) * exp(y)
-            winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
-
-         ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
-
-            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
-
-            x = log(16.0*bond*py/3.0)
-            y = 0.0
-
-            DO i = 1, 6
-               y = y + c(i) * (x**(i-1))
-            ENDDO
-
-            xrey = py*exp(y)
-
-            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
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of collision efficencies
-!------------------------------------------------------------------------------!
-
-   SUBROUTINE effic
-
-      USE arrays_3d
-      USE constants
-      USE cloud_parameters
-      USE particle_attributes
-
-!collision efficiencies of hall kernel
-      IMPLICIT NONE
-
-      INTEGER :: i, ir, iq, j, k, kk
-
-      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
-
-      LOGICAL, SAVE ::  first = .TRUE.
-
-      REAL ::  ek, particle_radius, pp, qq, rq
-
-      REAL, DIMENSION(1:21), SAVE :: rat
-      REAL, DIMENSION(1:15), SAVE :: r0
-      REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
-
-!
-!--   Initial assignment of constants
-      IF ( first )  THEN
-
-         first = .FALSE.
-         r0  = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., &
-                 300./)
-         rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,  &
-                 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
-
-         ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, &
-                        0.001,0.001,0.001,0.001,0.001,0.001/)
-         ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, &
-                        0.050,0.200,0.500,0.770,0.870,0.970/)
-         ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, &
-                        0.430,0.580,0.790,0.930,0.960,1.000/)
-         ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, &
-                        0.640,0.750,0.910,0.970,0.980,1.000/)
-         ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, &
-                        0.770,0.840,0.950,0.970,1.000,1.000/)
-         ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, &
-                        0.840,0.880,0.950,1.000,1.000,1.000/)
-         ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, &
-                        0.870,0.900,0.950,1.000,1.000,1.000/)
-         ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, &
-                        0.890,0.920,1.000,1.000,1.000,1.000/)
-         ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, &
-                        0.900,0.940,1.000,1.000,1.000,1.000/)
-         ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, &
-                        0.910,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, &
-                        0.920,0.950,1.000,1.000,1.000,1.000/)
-         ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, &
-                        0.930,0.970,1.000,1.000,1.000,1.000/)
-         ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, &
-                        0.950,1.000,1.000,1.000,1.000,1.000/)
-         ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, &
-                        1.000,1.020,1.020,1.020,1.020,1.020/)
-         ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, &
-                        1.030,1.040,1.040,1.040,1.040,1.040/)
-         ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, &
-                        1.700,2.300,2.300,2.300,2.300,2.300/)
-         ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, &
-                        3.000,4.000,4.000,4.000,4.000,4.000/)
-      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, 15
-            IF ( particle_radius < r0(k) )  THEN
-               ira(j) = k
-               EXIT
-            ENDIF
-         ENDDO
-         IF ( particle_radius >= r0(15) )  ira(j) = 16
-      ENDDO
-
-!
-!--   Two-dimensional linear interpolation of the collision efficiency.
-!--   Radius has to be in Âµm
-      DO  j = pstart, pend
-         DO  i = pstart, j
-
-            ir = ira(j)
-
-            rq = particles(i)%radius / particles(j)%radius
-
-!            DO kk = 2, 21
-!               IF ( rq <= rat(kk) )  THEN
-!                  iq = kk
-!                  EXIT
-!               ENDIF
-!            ENDDO
-
-            iq = INT( rq * 20 ) + 1
-            iq = MAX(iq , 2)
-
-            IF ( ir < 16 )  THEN
-
-               IF ( ir >= 2 )  THEN
-                  pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / &
-                       ( r0(ir) - r0(ir-1) )
-                  qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
-                  ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)   &
-                            + pp * ( 1.0-qq ) * ecoll(ir,iq-1)           &
-                            + qq * ( 1.0-pp ) * ecoll(ir-1,iq)           &
-                            + pp * qq * ecoll(ir,iq)
-               ELSE
-                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-                  ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
-               ENDIF
-
-            ELSE
-               qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
-               ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
-               ec(j,i) = MIN( ek, 1.0 )
-            ENDIF
-
-            ec(i,j) = ec(j,i)
-            IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
-
-         ENDDO
-      ENDDO
-
-      DEALLOCATE( ira )
-
-   END SUBROUTINE effic
-
-!------------------------------------------------------------------------------!
-! SUBROUTINE for calculation of enhancement factor collision efficencies
-!------------------------------------------------------------------------------!
-   SUBROUTINE turb_enhance_eff
+         r0  = (/ 6.0, 8.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60., &
+                  70.0, 100.0, 150.0, 200.0, 300.0 /)
+         rat = (/ 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, &
+                  0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, &
+                  1.00 /)
+
+         ecoll(:,1) = (/0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, &
+                        0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001/)
+         ecoll(:,2) = (/0.003, 0.003, 0.003, 0.004, 0.005, 0.005, 0.005, &
+                        0.010, 0.100, 0.050, 0.200, 0.500, 0.770, 0.870, 0.970/)
+         ecoll(:,3) = (/0.007, 0.007, 0.007, 0.008, 0.009, 0.010, 0.010, &
+                        0.070, 0.400, 0.430, 0.580, 0.790, 0.930, 0.960, 1.000/)
+         ecoll(:,4) = (/0.009, 0.009, 0.009, 0.012, 0.015, 0.010, 0.020, &
+                        0.280, 0.600, 0.640, 0.750, 0.910, 0.970, 0.980, 1.000/)
+         ecoll(:,5) = (/0.014, 0.014, 0.014, 0.015, 0.016, 0.030, 0.060, &
+                        0.500, 0.700, 0.770, 0.840, 0.950, 0.970, 1.000, 1.000/)
+         ecoll(:,6) = (/0.017, 0.017, 0.017, 0.020, 0.022, 0.060, 0.100, &
+                        0.620, 0.780, 0.840, 0.880, 0.950, 1.000, 1.000, 1.000/)
+         ecoll(:,7) = (/0.030, 0.030, 0.024, 0.022, 0.032, 0.062, 0.200, &
+                        0.680, 0.830, 0.870, 0.900, 0.950, 1.000, 1.000, 1.000/)
+         ecoll(:,8) = (/0.025, 0.025, 0.025, 0.036, 0.043, 0.130, 0.270, &
+                        0.740, 0.860, 0.890, 0.920, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,9) = (/0.027, 0.027, 0.027, 0.040, 0.052, 0.200, 0.400, &
+                        0.780, 0.880, 0.900, 0.940, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,10)= (/0.030, 0.030, 0.030, 0.047, 0.064, 0.250, 0.500, &
+                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,11)= (/0.040, 0.040, 0.033, 0.037, 0.068, 0.240, 0.550, &
+                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,12)= (/0.035, 0.035, 0.035, 0.055, 0.079, 0.290, 0.580, &
+                        0.800, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,13)= (/0.037, 0.037, 0.037, 0.062, 0.082, 0.290, 0.590, &
+                        0.780, 0.900, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,14)= (/0.037, 0.037, 0.037, 0.060, 0.080, 0.290, 0.580, &
+                        0.770, 0.890, 0.910, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,15)= (/0.037, 0.037, 0.037, 0.041, 0.075, 0.250, 0.540, &
+                        0.760, 0.880, 0.920, 0.950, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,16)= (/0.037, 0.037, 0.037, 0.052, 0.067, 0.250, 0.510, &
+                        0.770, 0.880, 0.930, 0.970, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,17)= (/0.037, 0.037, 0.037, 0.047, 0.057, 0.250, 0.490, &
+                        0.770, 0.890, 0.950, 1.000, 1.000, 1.000, 1.000, 1.000/)
+         ecoll(:,18)= (/0.036, 0.036, 0.036, 0.042, 0.048, 0.230, 0.470, &
+                        0.780, 0.920, 1.000, 1.020, 1.020, 1.020, 1.020, 1.020/)
+         ecoll(:,19)= (/0.040, 0.040, 0.035, 0.033, 0.040, 0.112, 0.450, &
+                        0.790, 1.010, 1.030, 1.040, 1.040, 1.040, 1.040, 1.040/)
+         ecoll(:,20)= (/0.033, 0.033, 0.033, 0.033, 0.033, 0.119, 0.470, &
+                        0.950, 1.300, 1.700, 2.300, 2.300, 2.300, 2.300, 2.300/)
+         ecoll(:,21)= (/0.027, 0.027, 0.027, 0.027, 0.027, 0.125, 0.520, &
+                        1.400, 2.300, 3.000, 4.000, 4.000, 4.000, 4.000, 4.000/)
+       ENDIF
+
+!
+!--    Calculate the radius class index of particles with respect to array r
+       ALLOCATE( ira(1:radius_classes) )
+       DO  j = 1, radius_classes
+          particle_radius = radclass(j) * 1.0E4   ! in microm
+          DO  k = 1, 15
+             IF ( particle_radius < r0(k) )  THEN
+                ira(j) = k
+                EXIT
+             ENDIF
+          ENDDO
+          IF ( particle_radius >= r0(15) )  ira(j) = 16
+       ENDDO
+
+!
+!--    Two-dimensional linear interpolation of the collision efficiency.
+!--    Radius has to be in Âµm
+       DO  j = 1, radius_classes
+          DO  i = 1, j
+
+             ir = ira(j)
+             rq = radclass(i) / radclass(j)
+             iq = INT( rq * 20 ) + 1
+             iq = MAX( iq , 2)
+
+             IF ( ir < 16 )  THEN
+                IF ( ir >= 2 )  THEN
+                   pp = ( ( radclass(j) * 1.0E04 ) - r0(ir-1) ) / &
+                        ( r0(ir) - r0(ir-1) )
+                   qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                   ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)  &
+                             + pp * ( 1.0-qq ) * ecoll(ir,iq-1)          &
+                             + qq * ( 1.0-pp ) * ecoll(ir-1,iq)          &
+                             + pp * qq * ecoll(ir,iq)
+                ELSE
+                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                   ec(j,i) = (1.0-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
+                ENDIF
+             ELSE
+                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
+                ec(j,i) = MIN( ek, 1.0 )
+             ENDIF 
+
+             ec(i,j) = ec(j,i)
+             IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
+
+          ENDDO
+       ENDDO
+
+       DEALLOCATE( ira )
+
+    END SUBROUTINE effic
+
+
+!------------------------------------------------------------------------------!
+! Calculation of enhancement factor for collision efficencies due to turbulence
+!------------------------------------------------------------------------------!
+    SUBROUTINE turb_enhance_eff
 
        USE constants
@@ -627,136 +713,131 @@
        USE arrays_3d
 
-      IMPLICIT NONE
-
-      INTEGER :: i, ik, ir, iq, j, k, kk
-
-      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 /)
-
-! 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 /)
-
-      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
-      DO j =  pstart, pend
-         DO i = pstart, j
-
-            ir = ira(j)
-
-            rq = particles(i)%radius/particles(j)%radius
-
-            DO kk = 2, 11
-               IF ( rq <= rat(kk) )  THEN
-                  iq = kk
-                  EXIT
-               ENDIF
-            ENDDO
-
-! 0  cm2/s3
-            y1 = 1.0
-! 100 cm2/s3, 400 cm2/s3
-            IF (ir.lt.8) THEN
-               IF (ir.ge.2) THEN
-                  pp = ((particles(j)%radius*1.0E06)-r0(ir-1))/(r0(ir)-r0(ir-1))
-                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-                  y2= (1.-pp)*(1.-qq)*ecoll_100(ir-1,iq-1)+  &
-                            pp*(1.-qq)*ecoll_100(ir,iq-1)+    &
-                            qq*(1.-pp)*ecoll_100(ir-1,iq)+    &
-                            pp*qq*ecoll_100(ir,iq)
-                  y3= (1.-pp)*(1.-qq)*ecoll_400(ir-1,iq-1)+  &
-                            pp*(1.-qq)*ecoll_400(ir,iq-1)+    &
-                            qq*(1.-pp)*ecoll_400(ir-1,iq)+    &
-                            pp*qq*ecoll_400(ir,iq)
-               ELSE
-                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-                  y2= (1.-qq)*ecoll_100(1,iq-1)+qq*ecoll_100(1,iq)
-                  y3= (1.-qq)*ecoll_400(1,iq-1)+qq*ecoll_400(1,iq)
-               ENDIF
-            ELSE
-            qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
-            y2 = (1.-qq) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
-            y3 = (1.-qq) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
-         ENDIF
-! linear interpolation
-! dissipation rate in cm2/s3
-         x1 = 0.0
-         x2 = 100.0
-         x3 = 400.0
-
-         IF (epsilon.le.100.) THEN
-            ecf(j,i) = (epsilon-100.)/(0.-100.) * y1 &
-                    +  (epsilon-0.)/(100.-0.) * y2
-         ELSE IF(epsilon.le.600.)THEN
-            ecf(j,i) = (epsilon-400.)/(100.-400.) * y2 &
-                    +  (epsilon-100.)/(400.-100.) * y3
-
-         ELSE
-            ecf(j,i) = (600.-400.)/(100.-400.) * y2 &
-                    +  (600.-100.)/(400.-100.) * y3
-         ENDIF
-
-         IF (ecf(j,i).lt.1.0) ecf(j,i) = 1.0
-
-         ecf(i,j)=ecf(j,i)
-      ENDDO
-   ENDDO
-
-   RETURN
-   END SUBROUTINE turb_enhance_eff
+       IMPLICIT NONE
+
+       INTEGER :: i, ik, iq, ir, j, k, kk
+
+       INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
+
+       REAL ::  particle_radius, pp, qq, rq, x1, x2, x3, y1, y2, y3
+
+       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.0, 20.0, 30.0, 40.0, 50.0, 60.0, 100.0 /)
+          rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /)
+!
+!--       In 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 /)
+!
+!--       In 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 /)
+
+       ENDIF
+
+!
+!--    Calculate the radius class index of particles with respect to array r0
+       ALLOCATE( ira(1:radius_classes) )
+
+       DO  j = 1, radius_classes
+          particle_radius = radclass(j) * 1.0E4    ! in microm
+          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 efficiencies
+       DO  j =  1, radius_classes
+          DO  i = 1, j
+
+             ir = ira(j)
+             rq = radclass(i) / radclass(j)
+
+             DO  kk = 2, 11
+                IF ( rq <= rat(kk) )  THEN
+                   iq = kk
+                   EXIT
+                ENDIF
+             ENDDO
+
+             y1 = 1.0      ! 0  cm2/s3
+!
+!--          100 cm2/s3, 400 cm2/s3
+             IF ( ir < 8 )  THEN
+                IF ( ir >= 2 )  THEN
+                   pp = ( radclass(j)*1.0E4 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
+                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                   y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) +  &
+                                pp * ( 1.0-qq ) * ecoll_100(ir,iq-1)   +  &
+                                qq * ( 10.-pp ) * ecoll_100(ir-1,iq)   +  &
+                                pp * qq         * ecoll_100(ir,iq)
+                   y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) +  &
+                                pp * ( 1.0-qq ) * ecoll_400(ir,iq-1)   +  &
+                                qq * ( 1.0-pp ) * ecoll_400(ir-1,iq)   +  &
+                                pp * qq         * ecoll_400(ir,iq)
+                ELSE
+                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                   y2 = ( 1.0-qq ) * ecoll_100(1,iq-1) + qq * ecoll_100(1,iq)
+                   y3 = ( 1.0-qq ) * ecoll_400(1,iq-1) + qq * ecoll_400(1,iq)
+                ENDIF
+             ELSE
+                qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
+                y2 = ( 1.0-qq ) * ecoll_100(7,iq-1) + qq * ecoll_100(7,iq)
+                y3 = ( 1.0-qq ) * ecoll_400(7,iq-1) + qq * ecoll_400(7,iq)
+             ENDIF
+!
+!--          Linear interpolation of dissipation rate in cm2/s3
+             IF ( epsilon <= 100.0 )  THEN
+                ecf(j,i) = ( epsilon - 100.0 ) / (   0.0 - 100.0 ) * y1 &
+                         + ( epsilon -   0.0 ) / ( 100.0 -   0.0 ) * y2
+             ELSEIF ( epsilon <= 600.0 )  THEN
+                ecf(j,i) = ( epsilon - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
+                         + ( epsilon - 100.0 ) / ( 400.0 - 100.0 ) * y3
+             ELSE
+                ecf(j,i) = (   600.0 - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
+                         + (   600.0 - 100.0 ) / ( 400.0 - 100.0 ) * y3
+             ENDIF
+
+             IF ( ecf(j,i) < 1.0 )  ecf(j,i) = 1.0
+
+             ecf(i,j) = ecf(j,i)
+
+          ENDDO
+       ENDDO
+
+    END SUBROUTINE turb_enhance_eff
 
  END MODULE lpm_collision_kernels_mod
