c     plug1.for
c     J.Phillips, USGS, Denver, Colorado
c
c     fills holes using the median of surrounding values
c
      dimension a(1000), b(1000, 3), c(1000), id(14), pgm(2)
      character ifile*50
      data ddval / 1.0e+36 /, dval / 1.0e+38 /
      print*,'PLUG1 can either:'
      print*,'  (1) add single-cell rinds to a grid or fill single-'
      print*,'      cell holes in a grid using the median of'
      print*,'      surrounding values, or'
      print*,'  (2) subtract single-cell rinds or delete isolated'
      print*,'      single-cell data values in a grid.'
      print*
    1 write(*, 200) 
  200 format(14h Enter ifile: ,$)
      read(5, 100) ifile
  100 format(a50)
      open(10, file=ifile, status='old', form='unformatted', err=1) 
      write(*, 201) 
  201 format(14h Enter ofile: ,$)
      read(5, 100) ifile
      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
      read(10) dlt, (b(i,3),i = 1, nc)
      write(11) dlt, (b(i,3),i = 1, nc)
      read(10) dlt, (b(i,2),i = 1, nc)
      print*
      print*,'Holes to be filled or enlarged are defined by the'
      print*,'minimum number of surrounding values required:'
      print*,'Specifying -8 will erase data in grid cells that are'
      print*,'              completely surrounded by dvals.'
      print*,'Specifying -7 will erase data in grid cells that are'
      print*,'              surrounded by at least 7 dvals.'
      print*,'Specifying -1 will erase a rind one grid cell thick from'
      print*,'              around the data area of the grid.'
      print*,'Specifying +1 will add a rind one grid cell thick'
      print*,'              around the data area of the grid.'
      print*,'Specifying +7 will fill single cell holes that are'
      print*,'              surround by at least 7 data values.'
      print*,'Specifying +8 will fill single cell holes that are'
      print*,'              completely surround by data.'
      print*
      write(*, 202)
  202 format(' What is the minimum number of surrounding values',/,
     1' required? (-8 to -1, 1 to 8): ',$)
      read(5,*) nv
      if(nv.lt.-8) nv=-8
      if(nv.gt.8) nv=8
      if(nv.eq.0) stop 'invalid response'
      if(nv.lt.0) then
        print*,'Requiring ',-nv,' surrounding dvals before erasing'
        print*,'  a data value.'
      else
        print*,'Requiring ',nv,' surrounding data values before'
        print*,'  interpolating a new data value.'
      endif
      nd = 0
      do 16 ii = 3, nr
        read(10) dlt, (b(i,1),i = 1, nc)
        c(1) = b(1,2)
        c(nc) = b(nc,2)
        do 14 i = 2, nc - 1
          c(i) = b(i,2)
          if (c(i) .ge. ddval .and. nv .gt. 0) then
            j = 0
            do 11 k = 1, 3
              do 10 l = i - 1, i + 1
                if (b(l,k) .lt. ddval) then
                  j = j + 1
                  a(j) = b(l,k)
                end if
10            continue
11          continue
            if (j .ge. nv) then
              call median(j, a)
                c(i) = a(1)
                nd=nd+1
            end if
          else if(c(i) .lt. ddval .and. nv .lt. 0) then
            j = 0
            do 13 k = 1, 3
              do 12 l = i - 1, i + 1
                if (b(l,k) .ge. ddval) then
                  j = j + 1
                end if
12            continue
13          continue
            if (j .ge. -nv) then
                c(i) = dval
                nd=nd+1
            end if
          end if
14      continue
        write(11) dlt, (c(i),i = 1, nc)
        do 15 i = 1, nc
          b(i,3) = b(i,2)
          b(i,2) = b(i,1)
15      continue
16    continue
      write(11) dlt, (b(i,1),i = 1, nc)
      if(nv.gt.0) write(*, 101) nd
  101 format(i10,' dvals plugged.')
      if(nv.lt.0) write(*, 102) nd
  102 format(i10,' dvals added.')
      stop 
      end
c
      subroutine median(n, a)
      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
