c g_screen.for
c  To remove or average close-together gravity stations.
c  Lin Cordell, 24 July 90
      character*50 ifile,ofile
      character*1 ans
      dimension id(14)
      dimension x(2000),y(2000),f(2000,6),l(2000)
c      data dval/0.170412e+39/
c  Read data.
      print*,'Enter input gravity post file name:'
      print 1
1     format(' *',$)
      read(5,800)ifile
800   format(a)
      open(10,file=ifile,form='unformatted',status='old')
      k=0
      xmin=1.e20
      xmax=-1.e20
      ymin=1.e20
      ymax=-1.e20
c      hmin=1.e20
c      hmax=-1.e20
c      fmin=1.e0
c      fmax=-1.e0
40    k=k+1
      if(k.gt.2000) then
        print*,'WARNING - processing only the first 2000 points'
        go to 41
      endif
      read(10,end=41)id(1),id(2),x(k),y(k),(f(k,j),j=1,6)
      if(x(k).gt.xmax)xmax=x(k)
      if(x(k).lt.xmin)xmin=x(k)
      if(y(k).gt.ymax)ymax=y(k)
      if(y(k).lt.ymin)ymin=y(k)
c      if(h(k).gt.hmax)hmax=h(k)
c      if(h(k).lt.hmin)hmin=h(k)
c      if(f(k).gt.fmax)fmax=f(k)
c      if(f(k).lt.fmin)fmin=f(k)
      goto 40
41    kmax=k-1
      close(10)
      do 71 k=1,kmax
71    l(k)=1
c      write(6,42)kmax,xmin,xmax,ymin,ymax,hmin,hmax,fmin,fmax
      write(6,42)kmax,xmin,xmax,ymin,ymax
42    format(i5,' Data points.'/' xmin,max = ',2e12.3,/' ymin,max'
c     &'= ',2e12.3,/' hmin,max = ',2e12.3,/' fmin,max = ',2e12.3)
     &'= ',2e12.3)
c  Read parameters.
      print*,'Enter proximity threshold (distance):'
      print 1
      read*,p
      psq=p*p
43    print*,'Enter test channel number (1-6):'
      print 1
      read*, itst
      if(itst.lt.1.or.itst.gt.6) go to 43
      print*,' Want manual or automatic mode?'
      print 1
      read(5,800)ans
      if(ans.eq.'a')then
       mode=2
       write(6,73)
73    format(' Enter difference limit (stations within the proximity'/
     &' limit having functional values in the test channel that are'/
     &' within the difference limit are averaged, otherwise both'/
     &' stations are deleted):')
       print 1
       read*,difflim
      else
       mode=1
      endif
c  Test for data proximity.
      rsqmin=1.e20
      do 51 k=1,kmax
      if(l(k).eq.0)go to 51
      do 55 j=1,kmax
      if(l(j).eq.0)go to 55
      if(j.eq.k)go to 55
      xsq=(x(k)-x(j))**2
      ysq=(y(k)-y(j))**2
      rsq=xsq+ysq
      if(rsq.gt.psq)go to 56
      r=sqrt(rsq)
      diff=abs(f(k,itst)-f(j,itst))
      if(mode.eq.1)go to 58
      ans='a'
      if(diff.gt.difflim)ans='b'
      go to 83
58    continue
      print*,k,j,itst
      print*,' f(i) =',f(k,itst),', f(j) =',f(j,itst)
      print*,' Distance =',r,', Difference =',diff
      write(6,57)
57    format(/' Average, remove both, remove i, remove j,'/
     &' or keep both (a,b,i,j,or k)?')
      print 1
      read(5,800)ans
83    if(ans.eq.'a')go to 59
      if(ans.eq.'b')go to 60
      if(ans.eq.'i')go to 61
      if(ans.eq.'j')go to 62
      if(ans.eq.'k')go to 56
      print*,'Answer error, run aborted.'
      go to 999
59    x(k)=0.5*(x(k)+x(j))
      y(k)=0.5*(y(k)+y(j))
      do i=1,6
c      h(k)=0.5*(h(k)+h(j))
      f(k,i)=0.5*(f(k,i)+f(j,i))
      end do
      l(j)=0
      go to 55
60    l(k)=0
      l(j)=0
      go to 55
61    l(k)=0
      go to 55
62    l(j)=0
      go to 55
56    if(rsq.lt.rsqmin)rsqmin=rsq
55    continue
51    continue
      kcount=0
      do 72 k=1,kmax
72    if(l(k).eq.1)kcount=kcount+1
      print*,'Number of points is',kcount
      rmin=sqrt(rsqmin)
      print*,'Minimum distance is',rmin
      print*,'Enter output file name:'
      print 1
      read(5,800)ofile
      open(11,file=ofile,form='unformatted',status='unknown')
c      open(10,file=ifile,form='unformatted',status='old')
      do 70 k=1,kmax
c      read(10)id(1),id(2),xk,yk,v1,fk,hk,v2,v3,v4
      if(l(k).eq.0)go to 70
c      write(11)id(1),id(2),xk,yk,v1,fk,hk,v2,v3,v4
      write(11)id(1),id(2),x(k),y(k),(f(k,i),i=1,6)
70    continue
999   continue
c      close(10)
      stop
      end
   
