c  terrace.for - dos version
c  changes from vms version are indicated by cv
c  vms version by L.Cordell and A.McCafferty, USGS, Denver, Colorado
c
c   To steepen gradients on basis of local curvature.
c
      character file1*50, file2*50, file3*50
      character pgm*8
      dimension a1(10000), a2(10000), a3(10000), id(14)
      real it, ib, il, ir, ic2, ic3
      data ddval / 1.e30 /
  111 write(6, 1) 
    1 format(1x,22hEnter input file name:/1x1h*$)
      read(5, 101) file1
  101 format(a50)
cv    open(10,file=file1,form='unformatted',status='old',iostat=ios)
cv    if(ios.ne.0) then
cv    open(10,file=file1,form='unformatted',status='old',readonly)
cv    endif
      open(10, file=file1, form='unformatted', status='old', err=111)
      write(6, 2) 
    2 format(1x,23hEnter output file name:/1x1h*$)
      read(5, 101) file2
cv    open (11,file='terrace.tmp',status='new',form='unformatted')
cv    open (12,file=file2,form='unformatted',status='new')
      open(11, file='terrace.tmp', status='unknown', 
     &form='unformatted') 
      open(12, file=file2, form='unformatted', status='unknown')

      file1='terrace.tmp'

      file2='terrace.tm1' 
      write(6, 7) 
    7 format(1x,27hEnter number of iterations;/1x1h*$)
      read(5, *) iterations
      write(6, 4) 
    4 format(1x,21hEnter output-file id:/1x1h*$)
      read(5, 5) (id(i),i = 1, 14)
    5 format(14a4)
      read(10) id, pgm, ncol, nrow, nz, x0, dx, y0, dy
      pgm = 'terrace '
      if ((nz .ne. 1) .or. (ncol .gt. 10000)) goto 99
      write(12) id, pgm, ncol, nrow, nz, x0, dx, y0, dy
      do 30 iter = 1, iterations
      icount = 0
      call rowio(a1, y, nz, ncol, 10, 1)
      call rowio(a2, y, nz, ncol, 10, 1)
      call rowio(a1, y, nz, ncol, 11, 2)
      do 10 j = 3, nrow
      call rowio(a3, y, nz, ncol, 10, 1)
      do 21 i = 2, ncol - 1
      if (a2(i) .ge. ddval) goto 13
      ic2 = 0.
      nc2 = 0
      if (a1(i) .lt. ddval) then
        it = a1(i)
        ic2 = ic2 + it
        nc2 = nc2 + 1
      else
        it = a2(i)
      end if
      if (a3(i) .lt. ddval) then
        ib = a3(i)
        ic2 = ic2 + ib
        nc2 = nc2 + 1
      else
        ib = a2(i)
      end if
      if (a2(i - 1) .lt. ddval) then
        il = a2(i - 1)
        ic2 = ic2 + il
        nc2 = nc2 + 1
      else
        il = a2(i)
      end if
      if (a2(i + 1) .lt. ddval) then
        ir = a2(i + 1)
        ic2 = ic2 + ir
        nc2 = nc2 + 1
      else
        ir = a2(i)
      end if
      ic3 = ic2 - (nc2 * a2(i))
      if (ic3) 12, 13, 14
   12 a1(i) = max(it,ib,il,ir,a2(i))
      goto 20
   13 a1(i) = a2(i)
      goto 20
   14 a1(i) = min(it,ib,il,ir,a2(i))
   20 if (a1(i) .eq. a2(i)) icount = icount + 1
   21 continue
      a1(1) = a2(1)
      a1(ncol) = a2(ncol)
      call rowio(a1, y, nz, ncol, 11, 2)
      do 11 i = 1, ncol
      a1(i) = a2(i)
   11 a2(i) = a3(i)
   10 continue
      call rowio(a2, y, nz, ncol, 11, 2)
      count = (icount * 100.) / float(ncol * nrow)
      write(6, 31) iter, count
   31 format(1x,11hIteration =,i5,27h   Percent of flat slopes =,f9.5)
      close(10) 
      close(11) 
cv    open (10,file='terrace.tmp',status='old',disp='delete',
cv   1	form='unformatted')
      open(10, file=file1, status='old', form='unformatted'
     &) 
      if (iter .eq. iterations) goto 30
cv    open (11,file='terrace.tmp',status='new',form='unformatted')
      open(11, file=file2, status='unknown',form='unformatted') 

      file3=file1

      file1=file2

      file2=file3
   30 continue
   33 continue
      do 32 j = 1, nrow
      call rowio(a1, y, nz, ncol, 10, 1)
   32 call rowio(a1, y, nz, ncol, 12, 2)
      goto 999
   99 write(6, 98) 
   98 format(1x,16hCan`t handle it.)
  999 continue
      stop 
      end
      subroutine rowio(a[HUGE], y, nz, ncol, ld, key)
      dimension a(nz, ncol)
      goto (1, 2), key
    1 read(ld) y, a
      goto 90
    2 write(ld) y, a
   90 return 
      end
