c     medifilt.for - dos version
c     J.Phillips, USGS, Denver, Colorado
c     changes from vms version are indicated by cv
c
c     median filtering using a 3x3, 5x5, 7x7, 9x9, or 11x11 window
c
      dimension a(10000), b(10000, 11), c(10000), id(14), pgm(2)
      character ifile*50
      data ddval / 1.0e+36 /
    1 write(*, 200) 
  200 format(14h enter ifile: ,$)
      read(5, 100) ifile
  100 format(a50)
cv    open(10,file=ifile,status='old',form='unformatted')
      open(10, file=ifile, status='old', form='unformatted', err=1) 
      write(*, 201) 
  201 format(14h enter ofile: ,$)
      read(5, 100) ifile
cv    open(11,file=ifile,status='new',form='unformatted')
      open(11, file=ifile, status='unknown', form='unformatted') 
      read(10) id, pgm, nc, nr, nz, xo, dx, yo, dy
      write(11) id, pgm, nc, nr, nz, xo, dx, yo, dy
  203 write(*, 202) 
  202 format(43h enter window index (1=3x3, 2=5x5, 3=7x7): ,$)
      read(5, *) iw
      if ((iw .lt. 1) .or. (iw .gt. 5)) goto 203
      do 10 j = (2 * iw) + 1, iw + 2, -1
        read(10) dlt, (b(i,j),i = 1, nc)
        write(11) dlt, (b(i,j),i = 1, nc)
10    continue

      do 11 j = iw + 1, 2, -1
        read(10) dlt, (b(i,j),i = 1, nc)
11    continue
      nd = 0
      do 18 ii = (2 * iw) + 1, nr
        read(10) dlt, (b(i,1),i = 1, nc)
        do 12 i = 1, nc
          c(i) = b(i,iw + 1)
12      continue
        do 15 i = iw + 1, nc - iw
          j = 0
          do 14 k = 1, (2 * iw) + 1
            do 13 l = i - iw, i + iw
              if (b(l,k) .lt. ddval) then
                j = j + 1
                a(j) = b(l,k)
              end if
13          continue
14        continue
          if (j .gt. 0) call median(j, a)
          if (j .lt. 1) then
            nd = nd + 1
          else
            c(i) = a(1)
          end if
15      continue
        write(11) dlt, (c(i),i = 1, nc)
        do 17 i = 1, nc
          do 16 j = (2 * iw) + 1, 2, -1
            b(i,j) = b(i,j - 1)
16        continue
17      continue
18    continue
      do 19 j = iw, 1, -1
        write(11) dlt, (b(i,j),i = 1, nc)
19    continue
      write(*, 101) nd
  101 format(i10,13h dvals remain)
      stop 
      end
c
      subroutine median(n, a[HUGE])
      dimension a(n)
   10 if (n .le. 2) return 
      ax = a(1)
      an = a(1)
      iax = 1
      ian = 1
      do 11 i = 2, n
        if (a(i) .gt. ax) then
          ax = a(i)
          iax = i
        end if
        if (a(i) .lt. an) then
          an = a(i)
          ian = i
        end if
11    continue
      if (an .eq. ax) then
        a(1) = an
        return 
      end if
      j = 0
      do 20 i = 1, n
        temp = a(i)
        if ((i .eq. iax) .or. (i .eq. ian)) goto 20
        j = j + 1
        a(j) = temp
   20 continue
      n = j
      goto 10
      end
