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