c     file: /export/software/puddle0/dems/timearea.f
c     Last Modified: 11/13/96
c
c     Written by: Glenn E. Moglen
c
c     Program calculates the time-area curve for a DEM basin given the 
c     area and drainage direction files.  Program reads the input file, timearea.in.
c
      program timearea
c
      parameter (igrid = 2404)
      integer*4 area(igrid, igrid), nx, ny
      integer*2 dir(igrid, igrid),  iroot, jroot
      real*4 dist(igrid, igrid), maxdist, dx, dy
      character*80 areafile, dirfile, file2
      open (15, file = 'timearea.in', form = 'formatted', status = 'old')
c      print *, 'Enter name of input area file: '
      read (15, 100) areafile
c      print *, 'Enter name of input direction file: '
      read (15, 100) dirfile
c      print *, 'Enter name of output file:'
      read (15, 100) file2
      read (15, *) ithresh   ! threshold area for channel formation (in pixels)
      read (15, *) npts      ! number of ordinates on time-area curve
      read (15, *) factor    ! velocity ratio hillslope / channel
      read (15, *) vout      ! velocity at outlet (m/s)
      close (15, status = 'keep')
      CALL demread (area, areafile, nx, ny, igrid, dx, dy)
      CALL demread (dir, dirfile, nx, ny, igrid, dx, dy)
      call getroot (area, nx, ny, iroot, jroot, igrid, vout)
      dir(jroot, iroot) = 9
      call getdist (area, dir, dist, nx, ny, dx, dy, iroot, jroot,
     +              maxdist, igrid, ithresh, factor, vout)
      call printout (area, dist, maxdist, nx, ny, file2, igrid, npts, dx, dy)
 100  format (a80)
      end
*
*
*   
      subroutine getroot(area, nx, ny, iroot, jroot, igrid, vout)
      integer*4 area(igrid, igrid), maxarea, nx, ny
      integer*2 iroot, jroot
      maxarea = 0
      do j = 1, ny
         do i = 1, nx
            if (area(j, i) .gt. maxarea) then
               maxarea = area(j, i)
               iroot = i
               jroot = j
            end if
         end do
      end do
      print *, 'Root at: ', iroot, jroot, maxarea
      vfact = vout / float(maxarea) ** 0.1
      print *, 'Velocity factor: ', vfact
      vout = vfact
      return
      end 
*
*
*
      subroutine getdist (area, dir, dist, nx, ny, dx, dy, iroot, jroot,
     +                    maxdist, igrid, ithresh, factor, vfact)
      integer*4 area(igrid, igrid), nx, ny
      integer*2 dir(igrid, igrid),  iroot, jroot
      real*4 dist(igrid, igrid), length, maxdist
      logical done
      print *, nx, ny, dx, dy, iroot, jroot, igrid
      do j = 1, ny
         do i = 1, nx
            dist(j, i) = -1.0
         end do
      end do
      dist(jroot, iroot) = 0.0
      maxdist = 0.0
      do j = 1, ny
         print *, j, ny, maxdist
         do i = 1, nx
            if (area(j, i) .gt. 0) then
               length = 0.0
               done = .false.
               iold = i
               jold = j
 1             continue
                  length = length + deltal(dir(jold, iold), dx, dy)
                  call nextpixel (iold, jold, dir, inew, jnew, igrid)
c
c  Vfactor = 1.0    for hillslope
c          = factor for channel
c
                  if (area(jold, iold) .lt. ithresh) then
                     vfactor = 1.0
                  else
                     vfactor = factor
                  end if
c
c  From Leopold & Maddock (1953) v = k Q ^ 0.1 
c  Assuming that Q ~= A
c
                  vfactor = vfactor * vfact * area(jold, iold) ** 0.1
c                 print *, vfactor, area(jold, iold), vfact
                  if (dist(jnew, inew) .ge. 0.0) then
                     dist(j, i) = dist(jnew, inew) + length * vfactor
                     done = .true.
                  else if (dir(jnew, inew) .eq. 9) then
                     print *, 'ERROR:', inew, jnew, area(jnew, inew)
                     do i4 = iold - 1, iold + 1
                        print *, (area(j4, i4), j4 = jold - 1, jold + 1)
                     end do
                     do i4 = iold - 1, iold + 1
                        print *, (dir(j4, i4), j4 = jold - 1, jold + 1)
                     end do
                     dist(j, i) = -1.0
                     done = .true.
c                     stop 
                  end if
               if (.not. done) then
                  iold = inew
                  jold = jnew
                  go to 1
               end if
               maxdist = max(dist(j, i), maxdist)
            end if
         end do
      end do
      do j = 1, ny
         do i = 1, nx
c           if (area(j, i) .gt. 0) print *, dist(j, i)
         end do
      end do
      return
      end
*
*
*
      real function deltal (dir, dx, dy)
      integer*2 dir
      if (dir .eq. 1 .or. dir .eq. 5) then
         deltal = dx
      elseif(dir .eq. 3 .or. dir .eq. 7) then
         deltal = dy
      else if (dir .eq. 2 .or. dir .eq. 4 .or. 
     +         dir .eq. 6 .or. dir .eq. 8) then
         deltal = (dx ** 2 + dy ** 2) ** 0.5
      else  
         deltal = 0.0
      end if
      return
      end
*
*
*
      subroutine nextpixel (i1, j1, dir, i2, j2, igrid)
      integer*2 dir(igrid, igrid), changex(9), changey(9)
      data changey/0, -1, -1, -1, 0, 1, 1, 1, 0/
      data changex/1, 1, 0, -1, -1 ,-1, 0, 1, 0/
      i2 = i1 + changex(dir(j1, i1)) 
      j2 = j1 + changey(dir(j1, i1))
      return
      end
*
*
*
      subroutine printout (area, dist, maxdist, nx, ny, file2, igrid,
     +                     npts, dx, dy)
      character*80 file2
      integer*2 group(4096), index
      integer*4 area(igrid, igrid), nx, ny
      real*4 dist(igrid, igrid), maxdist, distinc
      open (11, file = file2, form = 'formatted', status = 'unknown')
      distinc = maxdist / float(npts - 1)
      if (npts .gt. 4096) then
         print *, 'Need to redimension group(4096)', npts
         stop
      end if
      do i = 1, npts
         group(i) = 0
      end do
      do j = 1, ny
         do i = 1, nx
            if (area(j, i) .ge. 1) then
               index = int (dist(j, i) / distinc) + 1
               group(index) = group(index) + 1
            end if
         end do
      end do
      do i = 1, npts
         write (11, *) float(i) * distinc / 3600.0, group(i) * dx * dy / 100.0
      end do
      close (11, status = 'keep')
      return
      end