Listing 2: FORTRAN Program fdflow.for, implicit head-based finite difference method to solve [13] c fdflow implicit double precision (a-h, k, o-z) parameter (iz=500) dimension aa(0:iz), bb(0:iz), cc(0:iz), dd(0:iz) dimension th(0:iz), hj(0:iz), h(0:iz), k(0:iz), c(0:iz), z(0:iz) open (1, file='outflw.dat', status='unknown') zl = 0.2d0 tin = 3600.d0*2.d0 qin = -1.389d-6 tl = 3600.d0*12.d0 np = 100 dt = 10.0d0 hini = -100.d0 t = 0.d0 dz = zl/dble(np) rz = dt/dz rzz = dt/(dz*dz) do i = 0, np-1 h(i) = hini hj(i) = hini p = -h(i) call wprop (p, th(i), c(i), k(i)) j = np - i z(j) = (0.5d0-dble(i))*dz end do z(0) = -zl+0.5d0*dz write (*, *) ip = 0 do while ((t+dt) .le. tl) t = t + dt write (*, '(1h+, f10.2)') t do i = 1, np-2 km = 0.5d0*(k(i-1)+k(i)) kp = 0.5d0*(k(i)+k(i+1)) aa(i) = -rzz*km bb(i) = c(i) + rzz*(km+kp) cc(i) = -rzz*kp dd(i) = c(i)*hj(i) + rz*(kp-km) end do dd(1) = dd(1) + rzz*0.5d0*(k(0)+k(1))*hini km = 0.5d0*(k(np-2)+k(np-1)) aa(np-1) = -rzz*km bb(np-1) = c(i) + rzz*km cc(np-1) = 0.d0 dd(np-1) = c(np-1)*hj(np-1) -rz*(qin+km) call tridia (np-1, aa, bb, cc, dd) do i = 1, np-1 h(i) = dd(i) hj(i) = h(i) p = -h(i) call wprop (p, th(i), c(i), k(i)) end do if ((t. ge. tin) .and. (ip .eq. 0)) then c write out profile data after infiltration do i = 0, np-1 write (1, '(f7.4,f9.3,f7.4)') z(i), h(i), th(i) end do write (1, *) c set infiltration to zero ip = 1 qin = 0.05556d-6 end if end do write (*, *) c write out profile data after redistribution do i = 0, np-1 write (1, '(f7.4,f9.3,f7.4)') z(i), h(i), th(i) end do close (1) end subroutine wprop (p, th, c, k) implicit double precision (a-h, k, o-z) if (p .gt. 0.d0) then te = 1.d0/(1.d0 + (1.04d0*p)**1.3954d0) se = te**0.28335961d0 th = 0.106d0 + 0.3626d0*se c = 0.149106921d0*((1.04d0*p)**0.3954d0)*se*te + 1.d-20 sk = 1.d0-(1.d0-se**3.529084471d0)**0.28335961d0 k = 1.5162d-6*dsqrt(se)*sk*sk else th = .4686 k = 1.5162d-6 c = 1.d-20 end if return end subroutine tridia (nn, a, b, c, d) implicit double precision (a-h, o-z) parameter (iz=500) dimension a(0:iz), b(0:iz), c(0:iz), d(0:iz) do n = 2, nn r = a(n)/b(n-1) b(n) = b(n) - r*c(n-1) d(n) = d(n) - r*d(n-1) end do d(nn) = d(nn)/b(nn) do n = nn-1, 1, -1 d(n) = (d(n) - c(n)*d(n+1))/b(n) end do return end