source: palm/trunk/SOURCE/wang_kernel.f90 @ 790

MODULE wang_kernel_mod

!------------------------------------------------------------------------------!
! Current revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: wang_kernel.f90 790 2011-11-29 03:11:20Z raasch $
!
!
! 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 (...)
!------------------------------------------------------------------------------!

    USE arrays_3d
    USE cloud_parameters
    USE constants
    USE particle_attributes

    IMPLICIT NONE

    PRIVATE

    PUBLIC TurbSD, turb_enhance_eff, effic, fallg, PHI, ZHI, colker

    INTEGER, SAVE :: ip, jp, kp, pstart, pend, psum
    REAL, SAVE :: epsilon, urms
    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: gck, ec, ecf
    REAL, DIMENSION(:),   ALLOCATABLE, SAVE :: winf
    REAL, SAVE :: eps2, urms2

!
!-- Public interfaces
    INTERFACE TurbSD
       MODULE PROCEDURE TurbSD
    END INTERFACE TurbSD

    INTERFACE PHI
       MODULE PROCEDURE PHI
    END INTERFACE PHI

    INTERFACE ZHI
       MODULE PROCEDURE ZHI
    END INTERFACE ZHI

    INTERFACE turb_enhance_eff
       MODULE PROCEDURE turb_enhance_eff
    END INTERFACE turb_enhance_eff

    INTERFACE effic
       MODULE PROCEDURE effic
    END INTERFACE effic

    INTERFACE fallg
       MODULE PROCEDURE fallg
    END INTERFACE fallg

    INTERFACE colker
       MODULE PROCEDURE colker
    END INTERFACE colker

    CONTAINS

!------------------------------------------------------------------------------!
! SUBROUTINE for calculation of collision kernel
!------------------------------------------------------------------------------!
    SUBROUTINE colker(i1,j1,k1,kernel)

       USE arrays_3d
       USE cloud_parameters
       USE constants
       USE indices
       USE particle_attributes

       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

       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(winf(pstart:pend), ec(pstart:pend,pstart:pend))

       IF ( turbulence_effects_on_collision ) THEN

          ALLOCATE(gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend))

          epsilon = diss(kp,jp,ip) * 10000.0  !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, 1
                DO i =  pstart, pend, 1
                   kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * ec(i,j) * abs(winf(j)-winf(i))
                ENDDO
             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
          ENDIF

          DEALLOCATE(gck, ecf)

       ELSE

          CALL fallg
          CALL effic

          DO j = pstart, pend, 1
             DO i =  pstart, pend, 1
                kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * ec(i,j) * abs(winf(j)-winf(i))
             ENDDO
          ENDDO

       ENDIF

       DEALLOCATE(winf, ec)

       RETURN

    END SUBROUTINE colker

!------------------------------------------------------------------------------!
! SUBROUTINE for calculation of w, g and gck
!------------------------------------------------------------------------------!
    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
       USE cloud_parameters
       USE particle_attributes
       USE arrays_3d

       IMPLICIT NONE

       REAL :: airvisc, airdens, waterdens, gravity, anu, 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, DIMENSION(pstart:pend) :: vST, tau, St

       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

       tauk = (anu/epsilon)**0.5       !in s
       eta  = (anu**3.0/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
          vST(i) = winf(i)           !in cm/s
          tau(i) = vST(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.0)
       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.0)
       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 = vST(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
             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.0/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.0/t2    !in cm**2/s**2

              vrms2xy= sqrt(v2xysq)             !in cm/s

             IF(vST(i).ge.vST(j)) THEN
                v1 = vST(i)
                t1 = tau(i)
                v2 = vST(j)
                t2 = tau(j)
             ELSE
                v1 = vST(j)
                t1 = tau(j)
                v2 = vST(i)
                t2 = tau(i)
             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.0/tau(i)/tau(j)    !in cm**2/s**2

             wrtur2xy=vrms1xy**2.0 + vrms2xy**2.0 - 2.*v1v2xy  !in cm**2/s**2

              IF (wrtur2xy.le.0.0) wrtur2xy=0.0

              wrgrav2=pi/8.*(vST(j)-vST(i))**2.0

              wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s


!***** TO CALCULATE gr ********************************

             IF(St(j).gt.St(i)) THEN
                SSt = St(j)
             ELSE
                SSt = St(i)
             ENDIF

             XX = -0.1988*SSt**4.0 + 1.5275*SSt**3.0 &
                  -4.2942*SSt**2.0 + 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.0
             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.)
             IF (grFIN.lt.1.0) grFIN = 1.0

             gck(i,j) = 2. * pi * rrp**2.0 * wrFIN * grFIN       ! in cm**3/s
             gck(j,i) = gck(i,j)

          ENDDO
       ENDDO


       RETURN
    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.0  !in s
       RETURN
    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.0 - (1./tau2 + 1./a)**2.0
        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)   &
                                                           * 1./2./b/aa3             ! in s**2
        RETURN
    END FUNCTION ZHI


!------------------------------------------------------------------------------!
! SUBROUTINE for calculation of terminal velocity winf
!------------------------------------------------------------------------------!

   SUBROUTINE fallg

       USE constants
       USE cloud_parameters
       USE particle_attributes
       USE arrays_3d

      IMPLICIT NONE

      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, &
             -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.0) * (rhoa**2.0)/((eta**4.0) * grav * (rhow-rhoa))
      py   = phy**(1.0/6.0)

!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.0)
            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

            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
            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)
            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
            ENDIF

         ENDIF
      ENDDO
      RETURN
   END SUBROUTINE fallg

!------------------------------------------------------------------------------!
! SUBROUTINE for calculation of collision efficencies
!------------------------------------------------------------------------------!

   SUBROUTINE effic

       USE constants
       USE cloud_parameters
       USE particle_attributes
       USE arrays_3d

!collision efficiencies of hall kernel
      IMPLICIT NONE

      INTEGER :: i, ir, iq, j, k, kk
      REAL    :: ek, rq, pp, qq

      REAL, DIMENSION(1:21) :: rat
      REAL, DIMENSION(1:15) :: r0
      REAL, DIMENSION(1:15,1:21) :: ecoll

      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/)

! two-dimensional linear interpolation of the collision efficiency
! radius has to be in µm

      DO j = pstart, pend
         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
               ENDIF
            ENDDO

            rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06)

            DO kk = 2, 21
               IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
            ENDDO

            IF (ir.lt.16) 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))
                  ec(j,i) = (1.-pp) * (1.-qq) * ecoll(ir-1,iq-1) +  &
                          pp * (1.-qq) * ecoll(ir,iq-1) +         &
                          qq * (1.-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.-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).lt.1.e-20) stop 99
         ENDDO
      ENDDO
      RETURN
   END SUBROUTINE effic

!------------------------------------------------------------------------------!
! SUBROUTINE for calculation of enhancement factor collision efficencies
!------------------------------------------------------------------------------!

   SUBROUTINE turb_enhance_eff

       USE constants
       USE cloud_parameters
       USE particle_attributes
       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/)

! 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 /)

! two-dimensional linear interpolation of the collision efficiency
      DO j =  pstart, pend
         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
               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

! 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

 END MODULE wang_kernel_mod