Listing 2: tbnewtc3.for, Mass-Conservative Infiltration model, Arithmetic 
Mean, Adaptive DIRK2 Time Step 

Copies of this program can be made, used and distributed freely as long as 
they are whole and accompanied by a citation of this paper.  Modifications can 
be made and used likewise as long as they are accompanied in addition by a 
disclaimer that removes any liability from these authors. 

      program tbnewtc3

c     Celia's method with arithmetic mean & adaptive time step
c     & DIRK2 time step

      implicit real*8 (a-h, k, o-z)
      implicit integer (i-j, l-n)
      character fil*8
      parameter (idm=500)
      dimension k(0:idm), x(0:idm), h(0:idm), ht(0:idm),
     + th(0:idm), c(0:idm), thj(0:idm), hj(0:idm), thp(0:idm)
      common / r1 / tin, dx, dt, dx2, rxx, rmbtol, prk, a1, a2, flpt, 
     + flpb, flxt, flxb, sumth, dtt, rmb
      common / r2 / k, x, h, ht, th, c, thj, hj, thp
      common / i1 / np, maxit

      fil = 'au0598e0'
      xl = 120.d0
      np = 1000
      tl = 1000.d0
      dt = 1.d-3
      rmbtol = 1.d-8
      maxit = 4
      ptop = 20.7d0
      pbot = 929.8d0
      pini = pbot
      dtt = 1000.d0
      prk = 1.d0/3.d0
      a1 = 0.5d0/(1.d0-prk)
      a2 = (1.d0-2.d0*prk)*a1
      t = 0.d0 
      nt = 0
      tin = 0.d0
      dx = xl/dble(np)
      am = 0.d0
      ftt = 0.d0
      ftb = 0.d0
      dx2 = dx*dx
      
c      call timer (istr)
      open (1,file=fil//'.prn',status='unknown')
      write (1, *) '"'//fil//'"'
      write (1, '(f9.2/i6/f9.1)') xl, np, tl
      write (1, '(f6.2)') dtt
c      write (1, *)

      p = pbot
      h(0) = -p
      x(0) = 0.d0
      ht(0) = x(0) - pbot
      call wprop (p, k(0), th(0), c(0))
      do i = 1, np-1
         p = pini
         h(i) = -p
         x(i) = dble(i)*dx
         ht(i) = x(i) - p
         call wprop (p, k(i), th(i), c(i))
      end do
      p = ptop
      h(np) = -p
      x(np) = xl
      ht(np) = x(np) - ptop
      call wprop (p, k(np), th(np), c(np))
      do i = 0, np
         hj(i) = h(i)
         thj(i) = th(i)
         thp(i) = th(i)
      end do

      do while (t .lt. tl)
         nt = nt + 1
         if ((t+dt) .gt. tl) dt = tl-t

   10    continue
         t = t + dt
      
         do i = 1, np-1
            hj(i) = h(i)
         end do

         rxx = dt/dx2
         call pnewt (it)
         
         if (it .ge. (maxit+1)) then
            write (*, *) 
            write (*, *) 'restart iteration'
            write (*, *) 
            t = t - dt
            dt = 0.619d0*dt
            do i = 1, np-1
               h(i) = hj(i)
               p = -h(i)
               ht(i) = x(i) - p
               call wprop (p, k(i), th(i), c(i))
               thp(i) = th(i)
            end do
            if (dt .le. 1.d-100) then
               write (*, *) 'endless loop', t
               write (1, *) 'endless loop', t
               pause
               close (1)
               stop
            end if
            go to 10
         end if

         do i = 1, np-1
            th(i) = thp(i)
         end do

         call newt (it)
         dtl = dt
         if (it .le. (maxit-2)) then
            dt = 1.062d0*dt
            if (dt .gt. dtt) dt = dtt
         else if (it .eq. (maxit-1)) then 
            dt = 1.001d0*dt
            if (dt .gt. dtt) dt = dtt
         else if (it. eq. (maxit)) then
            dt = 0.619d0*dt
         else if (it .ge. (maxit+1)) then
            write (*, *) 
            write (*, *) 'restart iteration'
            write (*, *) 
            t = t - dt
            dt = 0.619d0*dt
            do i = 1, np-1
               h(i) = hj(i)
               p = -h(i)
               ht(i) = x(i) - p
               call wprop (p, k(i), th(i), c(i))
            end do
            if (dt .le. 1.d-100) then
               write (*, *) 'endless loop', t
               write (1, *) 'endless loop', t
               pause
               close (1)
               stop
            end if
            go to 10
         end if

         flx = (flxb - flxt)*dtl
         am = am + sumth
         ftt = ftt + flxt*dtl
         ftb = ftb + flxb*dtl
         ft = ftb - ftt
         cmb = (am-ft)/ft
         write (*, *)
         write (*, '(f9.2,i3,8g15.6)')
     +               t, it, dt, am, ftt, ftb, cmb, sumth, flx, rmb
         write (1, '(f9.2,i3,8g15.6)')
     +               t, it, dt, am, ftt, ftb, cmb, sumth, flx, rmb
c         write (*, *)
         do i = 1, np-1
            thj(i) = th(i)
            thp(i) = th(i)
         end do

      end do
      
c      write (1, *)
      write (1, *)
      write (1, '(g15.6)') am, ftt, ftb, cmb
      write (1, '(i7)') nt
      write (1, '(g16.7)') tin
      write (1, *)
      do i = 0, np
c      do i = 0, np, 2
         write (1, '(3(g15.6,1h,),g15.6)') (xl-x(i)), h(i), th(i), k(i)
      end do
      close (1)

      print *, fil, maxit, nt, tin
      print *, 'all done'
      pause 'hit return to finish'
            
      end


      subroutine tridia (nn, a, b, c, d)
      
      implicit integer (I-J), integer (L-N)
      implicit double precision (A-H), double precision (K)
      implicit double precision (O-Z)
      parameter (idm=500)
      dimension a(0:idm), b(0:idm), c(0:idm), d(0:idm)
      
      do 10 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)
   10 continue
      d(nn) = d(nn) / b(nn)
      do 20 n = nn - 1 , 1 , -1
          d(n) = (d(n) - c(n) * d(n + 1)) / b(n)
   20 continue
      
      return
      end


      subroutine wprop (p, k, th, c)

      implicit real*8 (a-h, k, o-z)
      implicit integer (i-j, l-n)

      if (p .gt. 0.d0) then
         k = 9.44d-3*(1.175d+6)/(1.175d+6+p**4.74d0)
         th = 1.611d+6*(2.12d-1)/(1.611d+6 + p**3.96d0)
         c = 3.96d0*th*(p**2.96d0)/
     +           (1.611d+6 + p**3.96d0)
         th = th + 0.075d0
      else
         k = 0.9999999999d0*9.44d-3
         th = 0.287d0
         c = 1.d-15
      end if

      return
      end


      subroutine pnewt (it)

      implicit real*8 (a-h, k, o-z)
      implicit integer (i-j, l-n)
      parameter (idm=500)
      dimension k(0:idm), x(0:idm), h(0:idm), ht(0:idm),
     + th(0:idm), c(0:idm), thj(0:idm), hj(0:idm), thp(0:idm),
     + aa(0:idm), bb(0:idm), cc(0:idm), dd(0:idm)
      common / r1 / tin, dx, dt, dx2, rxx, rmbtol, prk, a1, a2, flpt,
     + flpb, flxt, flxb, sumth, dtt, rmb
      common / r2 / k, x, h, ht, th, c, thj, hj, thp
      common / i1 / np, maxit
      
      rxp = prk*rxx
      do it = 1, maxit+1

         tin = tin + 1.d0
         do i = 1, np-1
            aa(i) = -0.5d0*(k(i)+k(i-1))*rxp
            cc(i) = -0.5d0*(k(i+1)+k(i))*rxp
            bb(i) = c(i) -aa(i) -cc(i)
            dd(i) = (thp(i)-thj(i))
     +              -0.5d0*( (k(i+1)+k(i))*(ht(i+1)-ht(i))
     +              -(k(i)+k(i-1))*(ht(i)-ht(i-1)) )*rxp
            dd(i) = -dd(i)
         end do
            
         call tridia (np-1, aa, bb, cc, dd)

         sumth = 0.d0
         do i = 1, np-1
            h(i) = h(i) + dd(i)
            p = -h(i)
            ht(i) = x(i) - p
            call wprop (p, k(i), thp(i), c(i))
            sumth = sumth + thp(i) - thj(i)
         end do

         sumth = sumth*dx
         flpt = -0.5d0*(k(np)+k(np-1))*(ht(np)-ht(np-1))/dx
         flpb = -0.5d0*(k(0)+k(1))*(ht(1)-ht(0))/dx
         flx = (flpb-flpt)*dt*prk
         rmb = (sumth-flx)/flx
         write (*, '(1h+,i5,2g16.6)') it, rmb, dt
         if (dabs(rmb) .le. rmbtol) go to 100

      end do
      it = it - 1

  100 continue
      write (*, *)

      return
      end


      subroutine newt (it)

      implicit real*8 (a-h, k, o-z)
      implicit integer (i-j, l-n)
      parameter (idm=500)
      dimension k(0:idm), x(0:idm), h(0:idm), ht(0:idm),
     + th(0:idm), c(0:idm), thj(0:idm), hj(0:idm), thp(0:idm),
     + aa(0:idm), bb(0:idm), cc(0:idm), dd(0:idm)
      common / r1 / tin, dx, dt, dx2, rxx, rmbtol, prk, a1, a2, flpt,
     + flpb, flxt, flxb, sumth, dtt, rmb
      common / r2 / k, x, h, ht, th, c, thj, hj, thp
      common / i1 / np, maxit
      
      rx2 = rxx*a2
      ap = a1/prk
      do it = 1, maxit+1

         tin = tin + 1.d0
         do i = 1, np-1
            aa(i) = -0.5d0*(k(i)+k(i-1))*rx2
            cc(i) = -0.5d0*(k(i+1)+k(i))*rx2
            bb(i) = c(i) -aa(i) -cc(i)
            dd(i) = th(i)-thj(i) -ap*(thp(i)-thj(i))
     +              -0.5d0*( (k(i+1)+k(i))*(ht(i+1)-ht(i))
     +              -(k(i)+k(i-1))*(ht(i)-ht(i-1)) )*rx2
            dd(i) = -dd(i)
         end do
            
         call tridia (np-1, aa, bb, cc, dd)

         sumth = 0.d0
         do i = 1, np-1
            h(i) = h(i) + dd(i)
            p = -h(i)
            ht(i) = x(i) - p
            call wprop (p, k(i), th(i), c(i))
            sumth = sumth + th(i) - thj(i)
         end do

         sumth = sumth*dx
         flxt = -0.5d0*(k(np)+k(np-1))*(ht(np)-ht(np-1))/dx
         flxt = a1*flpt + a2*flxt
         flxb = -0.5d0*(k(0)+k(1))*(ht(1)-ht(0))/dx
         flxb = a1*flpb + a2*flxb
         flx = (flxb-flxt)*dt
         rmb = (sumth-flx)/flx
         write (*, '(1h+,i5,2g16.6)') it, rmb, dt
         if (dabs(rmb) .le. rmbtol) go to 100

      end do
      it = it - 1

  100 continue

      return
      end