c
c     file: /export/software/puddle0/dems/slopearea.f
c     Last Modified: 11/19/96
c
c     Written by: Glenn E. Moglen
c
c     Program calculates the slope-area distribution
c
      program slopearea
      integer*2 d1(8), d2(8)
      real*4 dx, dy, vec(400000, 2)
      integer*2 dir(1450, 1450), elev(1450, 1450)
      integer*4 area(1450, 1450), nx, ny
      character*80 areafile, elevfile, dirfile, trash
      data d1/0, -1, -1, -1, 0, 1, 1, 1/
      data d2/1, 1, 0, -1, -1, -1, 0, 1/
      igridy = 1450
      open (15, file = 'branch.in', form = 'formatted', status = 'old')
      read (15, 100) trash
      read (15, 100) dirfile
      read (15, 100) areafile
      read (15, 100) elevfile
      close (15, status = 'keep')
      call demread (area, areafile, nx, ny, igridy, dx, dy)
      call demread (dir, dirfile, nx, ny, igridy, dx, dy)
      call demread (elev, elevfile, nx, ny, igridy, dx, dy)
      k = 1
      do i = 1, ny
         do j = 1, nx
            if (area(i, j) .gt. 0) then
               ii = i + d1(dir(i, j))
               jj = j + d2(dir(i, j))
               run = ((dy*d1(dir(i, j))) ** 2 + (dx*d2(dir(i,j))) ** 2) ** 0.5
               slope = float(elev(i, j) - elev(ii, jj)) / run
               if (slope .lt. 0.0) then
                  print *, i, j, dir(i, j), elev(i, j), elev(ii, jj)
c                 stop
               end if 
               vec(k, 1) = float(area(i, j))
               vec(k, 2) = slope
               k = k + 1
               if (k .gt. 400000) go to 20
            end if
         end do
      end do
 20   n = k - 1
c     print *, 'Number of points read: ', n
      call sort (vec, n)
      call group (vec, n)
 100  format(a80)
      end
*
*
*
      subroutine group (a, n)
      real*4 a(400000, 2), g(1000, 2)
      integer*2 binsize 
      integer*4 kk
      binsize = 100
      ig = 0
      kk = 1
 20   if ((n - kk) .gt. binsize) then 
          ig = ig + 1
          s1 = 0.0
          s2 = 0.0
          nobin = 0
 30       if (nobin .lt. binsize) then
             s1 = s1 + a(kk, 1)
             s2 = s2 + a(kk, 2)
             nobin = nobin + 1
             kk = kk + 1
             go to 30
          end if
          g(ig, 1) = s1 / float(binsize)
          g(ig, 2) = s2 / float(binsize)
          go to 20  
      end if   
      nobin = ig - 1
      do i = 1, ig
         print *, g(i, 1), g(i, 2)
      end do
      return
      end 
*
*
*
      subroutine sort (a, n)
      real*4 a(400000, 2)
c     print *, 'Sorting: '
      do i = 1, n - 1
c        if (mod(i, 1000) .eq. 0) print *, i, n
         do j = i + 1, n
            if (a(i, 1) .gt. a(j, 1)) then
               temp = a(i, 1)
               a(i, 1) = a(j, 1)
               a(j, 1) = temp
               temp = a(i, 2)
               a(i, 2) = a(j, 2)
               a(j, 2) = temp
            end if
         end do
      end do
      return
      end