c  spike.for
c
c  Spike filter on a grid.  Output is the difference between the value
c  of the array at a point and the mean value on the border of a circular
c  window centered on the point. Dvals are ok.
c
c  Lin Cordell, 11 May 91.
c
      character*50 ifile,ofile
      dimension a(250,250),ar(250),id(14),pgm(2)
      data dval/0.170412e+39/
      dval2=1.e20
      print*,' Enter input file name:'
      print 1
1     format(' *'$)
      read(5,2)ifile
2     format(a)
      open(10,file=ifile,form='unformatted',status='old')
      read(10)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
      if(nz.ne.1.or.dx.eq.0.e0) go to 900
      if(ncol.gt.250.or.nrow.gt.250)go to 900
      if(dx.ne.dy)go to 900
      print*,'Enter radius of the window:'
      print 1
      read*,r
      if(dx.gt.r.or.dy.gt.r)then
        print*,'The window must be larger than one grid interval.'
        print*,'Run aborted'
        go to 998
      endif
      print*,'Enter output file name:'
      print 1
      read(5,2)ofile
      open(11,file=ofile,form='unformatted',status='new')
      print*,' Enter new id:'
      print 1
      read(5,3)(id(i),i=1,14)
3     format(14a4)
      write(11)id,pgm,ncol,nrow,nz,x0,dx,y0,dy
c  Read the data grid
      do 5 j=1,nrow
5     read(10) dum,(a(i,j),i=1,ncol)
c  Computations
      r=r/dx
      ir=r+0.5
      rsq=ir*ir
      yline=y0-dy
      do 200 j=1,nrow
      yline=yline+dy
      do 100 i=1,ncol
      if(a(i,j).gt.dval2)then
        ar(i)=dval
        go to 100
      endif
c  Work border of the circle with radius r.
      n=0
      sum=0
      do 20 ii=0,ir
      x=ii
      y=sqrt(rsq-x*x)
      jj=y+0.5
      ib=i+ii
      if(ib.gt.ncol)go to 10
      jb=j-jj
      if(jb.gt.1)then
        if(a(ib,jb).gt.dval2)then
          ar(i)=dval
        else
          n=n+1
          sum=sum+a(ib,jb)
        endif
      endif
      jb=j+jj
      if(jb.lt.nrow)then
        if(a(ib,jb).gt.dval2)then
          ar(i)=dval
        else
          n=n+1
          sum=sum+a(ib,jb)
        endif
      endif
10    continue
      ib=i-ii
      if(ib.lt.1)go to 11
      jb=j-jj
      if(jb.gt.1)then
        if(a(ib,jb).gt.dval2)then
          ar(i)=dval
        else
          n=n+1
          sum=sum+a(ib,jb)
        endif
      endif
      jb=j+jj
      if(jb.lt.nrow)then
        if(a(ib,jb).gt.dval2)then
          ar(i)=dval
        else
          n=n+1
          sum=sum+a(ib,jb)
        endif
      endif
11    continue
20    continue
      if(n.ne.0)then
        rn=n
        ar(i)=a(i,j)-(sum/rn)
      else
        ar(i)=dval
      endif
100   continue
      call rowio(ar,yline,nz,ncol,11,2)
200   continue
998   close(10)
      close(11)
      go to 999
900   print*,'Cant handle it.'
      go to 999
999   stop
      end
      subroutine rowio(a,y,nz,ncol,ld,key)
      dimension a(nz,ncol)
      go to (1,2),key
1     read(ld) y,a
      go to 90
2     write(ld)y,a
90    return
      end
