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

MODULE wang_kernel_mod

!------------------------------------------------------------------------------!
! Current revisions:
! -----------------
! 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:
! -----------------
! $Id: wang_kernel.f90 799 2011-12-21 17:48:03Z franke $
!
! 790 2011-11-29 03:11:20Z raasch
! initial revision
!
! 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  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 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

    CONTAINS

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

       USE arrays_3d
       USE cloud_parameters
       USE constants
       USE cpulog
       USE indices
       USE interfaces
       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

!       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 ( turbulence_effects_on_collision )  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
             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 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
!------------------------------------------------------------------------------!
    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 :: 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

       LOGICAL, SAVE ::  first = .TRUE.

       INTEGER :: i, j

!
!--   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/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
       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
                v1 = winf(i)
                t1 = tau(i)
                v2 = winf(j)
                t2 = tau(j)
             ELSE
                v1 = winf(j)
                t1 = tau(j)
                v2 = winf(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/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)
             ELSE
                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
             gck(j,i) = gck(i,j)

          ENDDO
       ENDDO

    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)

       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

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

         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

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

 END MODULE wang_kernel_mod