      character id*8, ifile*80, ans*1, ans1*1, prj*1
      dimension post(6)
      pi=3.1415927
      rad=pi/180.
      print*,'Input post file:'
      read(5,100) ifile
  100 format(a)
      open(10,file=ifile,form='unformatted',status='old')
      print*,'Output .raw file:'
      read(5,100) ifile
      open(11,file=ifile,form='formatted',status='unknown')
      print*,'Is post file projected?:'
      read(5,100) prj
      if(prj.eq.'Y') prj='y'
      rfac=360./(4.*pi*6371.)
      print*,'Do you want to plot strike symbols?:'
      read(5,100) ans
      if(ans.eq.'y'.or.ans.eq.'Y') then
        print*,'Length of strike symbols in km (or'
        print*,'  units of dx for None projection):'
        read(5,*) bkm
        b=bkm*rfac
      else
        bkm=-1.
        b=-1.
      endif
      print*,'Do you want to plot depth circles?:'
      read(5,100) ans
      if(ans.eq.'y'.or.ans.eq.'Y') then
        print*
        print*,'The circle diameter will be set equal to'
        print*,'  factor * (depth + reference_elevation).'
        print*
        print*,'Circles will not be plotted for depths above'
        print*,'  the reference elevation.'
        print*
        print*,'Enter the reference elevation in same units as dx:'
        read(5,*) elev
        print*,'Enter the factor:'
        read(5,*) factor
        elev=-elev
        print*,'Enter maximum percent error in depth:'
        read(5,*) pct
        dmax=-1.e38
        dmin=1.e38
        dave=0.
        nn=0
        mm=0
   10   read(10,end=15) id,alon,alat,post
        mm=mm+1
        if(post(2).le.pct) then
          if(post(1).gt.dmax) dmax=post(1)
          if(post(1).lt.dmin) dmin=post(1)
          dave=dave+post(1)
          nn=nn+1
        endif
        go to 10
   15   rewind(10)
        dave=dave/float(nn)
        print*
        print*,'reference elevation  = ',elev
        print*,'   number of strikes = ',mm
        print*,'   number of depths  = ',nn
        print*,'      minimum depth  = ',dmin
        print*,'      average depth  = ',dave
        print*,'      maximum depth  = ',dmax
        print*
c        print*,'set new maximum depth?:'
c        read(5,100) ans1
c        if(ans1.eq.'y'.or.ans1.eq.'Y') then
c          print*,'enter new maximum depth:'
c          read(5,*) dmax
c        endif
c        range=dmax-elev
c         range=2.*bkm
      endif
        print*,'Enter min,max susc contrast:'
        read(5,*) smin,smax
      nn=0
      mm=0
      itype=1
   20 read(10,end=99) id,alon,alat,post
      if(post(6).lt.smin.or.post(6).gt.smax) go to 20
      mm=mm+2
c      if(mm.gt.10987) stop 'too many vectors for reorder'
      if(mm.gt.10987) then
        print*,'Output file full'
        close(11)
        print*,'New .raw file:'
        read(5,100) ifile
        open(11,file=ifile,form='formatted',status='unknown')
        mm=1
      endif
      if(prj.ne.'y') then
        car=cos(alat*rad)
        if(car.eq.0.0) car=1.0
        if(b.ge.0.0) then
          theta=rad*post(3)
          stheta=b*sin(theta)
          stheta=stheta/car
          ctheta=b*cos(theta)
          write(11,*) itype
          write(11,*) alat-ctheta,alon-stheta
          write(11,105) alat+ctheta,alon+stheta
          theta=rad*post(4)
          stheta=b*sin(theta)
          stheta=stheta/car
          ctheta=b*cos(theta)
          write(11,*)itype
          write(11,*) alat,alon
          write(11,105) alat+ctheta,alon+stheta
        endif
        if(ans.eq.'y'.or.ans.eq.'Y') then
          if(post(2).le.pct) then
            radius=rfac*factor*(post(1)-elev)
            if(radius.le.0.0) then
              nn=nn+1
              go to 20
            endif
            write(11,*) itype
            write(11,*) alat, alon+radius/car
            do 30 i=30,330,30
            ang=rad*i
            write(11,*) alat+radius*sin(ang), alon+radius*cos(ang)/car
   30       continue
            write(11,105) alat, alon+radius/car
          endif
        endif
      else
        if(bkm.ge.0.0) then
          theta=rad*post(3)
          stheta=bkm*sin(theta)/2.
          ctheta=bkm*cos(theta)/2.
          write(11,*) itype
          write(11,*) (alat-ctheta)*1000.,(alon-stheta)*1000.
          write(11,105) (alat+ctheta)*1000.,(alon+stheta)*1000.
          theta=rad*post(4)
          stheta=bkm*sin(theta)/2.
          ctheta=bkm*cos(theta)/2.
          write(11,*) itype
          write(11,*) alat*1000.,alon*1000.
          write(11,105) (alat+ctheta)*1000.,(alon+stheta)*1000.
        endif
        if(ans.eq.'y'.or.ans.eq.'Y') then
          if(post(2).le.pct) then
            radius=factor*(post(1)-elev)/2.
            if(radius.le.0.0) then
              nn=nn+1
              go to 20
            endif
            write(11,*) itype
            write(11,*) alat*1000., (alon+radius)*1000.
            do 40 i=30,330,30
            ang=rad*i
            write(11,*) (alat+radius*sin(ang))*1000.,
     1                  (alon+radius*cos(ang))*1000.
   40     continue
          write(11,105) alat*1000., (alon+radius)*1000.
          endif
        endif
      endif
      go to 20
   99 if(nn.gt.0) then
        print*,'Warning: ',nn,' depths above reference elevation.'
      endif
      stop
105   format(2f16.6,'//')
      end
