Listing 2: FORTRAN code to calculate the integral form of Kh in 
      [7] using single-segment 10th-order Gaussian quadrature 
      
c     The variables p1 and p1 are the matric suctions, y1 and y2, at the 
c     centers of vertically adjacent grid cells.  The function, ptc(p), 
c     returns the relative conductivity, kr, for a given matric suction, 
c     p > 0.  The calculated value of Kh is the returned function value, 
c     dink10.  There is not apparent division by Dy in the routine 
c     because it cancels the dy term in single-segment numerical 
c     integration. 

      double precision function dink10 (p1, p2)

      implicit real*8 (a-h, k, o-z)
      dimension g(5), w(5)

      data (g(i), i=1,5) / 0.1488743389d0, 0.4333953941d0,
     + 0.6794095682d0, 0.8650633666d0, 0.9739065285d0 /
      data (w(i), i=1,5) / 0.2955242247d0, 0.2692667193d0,
     + 0.2190863625d0, 0.1494513491d0, 0.0666713443d0 /

      if (p1 .eq. p2) then
         dink10 = ptc(p1)
         return
      end if
      dink10 = 0.d0
      pm = .5d0*(p1+p2)
      pr = .5d0*(p2-p1)
      do i = 1, 5
         pg = pr*g(i)
         dink10 = dink10 + w(i)*( ptc(pm+pg) + ptc(pm-pg) )
      end do
      dink10 = 0.5d0*dink10
      
      return
      end