C
C________________________________________________________________
C
C     SUBROUTINE  C O N T T R
C________________________________________________________________
C
      subroutine conttr(x,y,z,f,c,nrow,ncol,xxscal,yyscal,
     & gflg,yo,dely,iquad,npass,jres,igrid,ixad,iyad,ierr)
c
c     basic input and tier control
c
      dimension x(1),y(1),z(1),c(1)
      integer f(1)

      ierr=0
c
c     setup mode control
c
      if(iquad.eq.3) go to 10
      call coorda(ixad,iyad,ncol,iquad)
      itype=-1
      if(dely.eq.0) itype=0
      go to 20
   10 call coorda(ixad,iyad,ncol,iquad)
      itype=1
c
c     remaining initialization
c
   20 klrow=(nrow-1)*ncol
      last=0
      je=nrow
c
c     tier loop
c
      do 190 ipass=1,npass
      if(ipass.eq.npass) je=jres
      if(ipass.eq.1) go to 70
c
c move down last row
c
      js=2
      jj=ncol
      do 30 i=1,ncol
   30 z(i)=z(klrow+i)
      if(itype) 40,40,50
   40 y(1)=y(nrow)
      go to 80
   50 do 60 i=1,ncol
      x(i)=x(klrow+i)
   60 y(i)=y(klrow+i)
      go to 80
c
c     first row
c
   70 js=1
      jj=0
c
c     get data
c
   80 do 140 j=js,je
      if(itype) 90,100,110
   90 read(igrid,end=130) dum,(z(jj+i),i=1,ncol)
      y(j)=yo
      yo=yo+dely
      go to 140
  100 read(igrid,end=130) y(j),(z(jj+i),i=1,ncol)
      y(j)=y(j)*yyscal
      go to 140
  110 read(igrid,end=130) dum,(x(jj+i),y(jj+i),z(jj+i),i=1,ncol)
      do 120 i=1,ncol
      kk=jj+i
      if(z(kk).gt.gflg) go to 120
      x(kk)=x(kk)*xxscal
      y(kk)=y(kk)*yyscal
  120 continue
      go to 140
  130 je=j-1
      if(je.lt.2) return
      last=1
      go to 150
  140 jj=jj+ncol
c
c     call contr
c
  150 call contr(z,f,c,c,je,ipass,npass,ierr)
      if(ierr.ne.0) return
      CALL grdegg
  180 if(last.ne.0) go to 200
  190 continue
  200 call perout(0,1000,npass+1,npass,'COMPLETED')
      return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  C O N T R
C________________________________________________________________
C
      subroutine contr(grid,flags,work,acval,nrowf,ipass,npass,
     & ierr)
c
c  basic contouring control subroutine.
c
      common /contrc/ cmin,cmax,dcval,ncval,nsec,gradi,
     1                acdel(200),acgrad(200),JCDASH(200),
     2                ncval2,acsize(200),acmin(200),acmax(200)
      common /concom/ ncol,nrow,bmin,bmax,
     & grad,gflg,ijs,ije,ijsi,ijei,prime,ij4,cont,
     & fltmax,lmult(0:3),idashs,linet,sigma,nsig
      common /setxy/ xx(2),yy(2),xp1(4),yp1(4),iplotr,sizel,
     & ncharx,nchary,sizex,sizey,adelx,adely,pllx,plly,
     & lintx,linty,xscale,yscale,xxscal,yyscal,mscale,delx,dely
      common /pltcom/ npens
      common/labcom/ dum1(11),size,dum2(4)
      dimension grid(0:*),flags(0:*),work(1),acval(1)
      integer flags
      logical prime

c
c  set gradient.
      nrow=nrowf
c
C
C IF ACVALS AND ACGRAD ARE BEING USED DO SETUP ONLY WHEN
C ACGRAD CHANGES
C
      IF(acgrad(1).NE.0.AND.NCVAL.NE.0) THEN
        nout=0
        ncont = ncval
        PRIME=.TRUE.
        grdtmp = gradi
        GRADI=-1
        GRAD0=acgrad(1)
        GRAD=GRAD0*GRAD0
        CALL SETUP(GRID,FLAGS,WORK)
        IF(BMIN.GE.BMAX) GOTO 140
        do 1 i=1,ncval2
          if(acgrad(i).ne.grad0) then
            grad0=acgrad(i)
            grad=grad0*grad0
            call setup(grid,flags,work)
          end if
          if(jcdash(1).ne.999) linet=iabs(jcdash(i))
          size=acsize(i)
          cmin=acmin(i)
          cmax=acmax(i)
          dcval=abs(acdel(i))
          if(cmin.ne.0..or.cmax.ne.0.) then
            temin=amax1(bmin,cmin)/dcval
            icont=int(temin)
            if(aint(temin).ne.temin.and.temin.gt.0) icont=icont+1
c           icont=amax1(bmin,cmin)/dcval
            amx=amin1(bmax,cmax)
          else
            temin=bmin/dcval
            icont=int(temin)
            if(aint(temin).ne.temin.and.temin.gt.0) icont=icont+1
c           icont=bmin/dcval
            amx=bmax
          end if
          cont=icont*dcval
    5 if(cont.gt.amx) go to 1
            do 2 j=1,i-1
              if(acmin(j).eq.0..and.acmax(j).eq.0.) then
                if(errbar(cont,acdel(j)).eq.0.) go to 201
              else if(cont.ge.acmin(j).and.cont.le.acmax(j)) then
                if(errbar(cont,acdel(j)).eq.0.) go to 201
              endif
c             temod=abs(amod(cont,acdel(j)))
c             if(temod.lt.1.e-6) goto 201
c             if(acdel(j)-temod.lt.1.e-6) goto 201
    2       continue
            do 3 j=ncval2+1,ncval
              if(acval(j).eq.cont) goto 201
    3       continue
            
            call setlab
            call scan(grid,flags,work,ierr)
c***********OUTPUT PERCENTAGE COMPLETED
            nout=nout+1
            call perout(nout,ncont,ipass,npass,'COMPLETED')
            if(ierr.ne.0) return
  201       icont=icont+1
            cont=icont*dcval
            go to 5
    1 continue
        do 6 i=ncval2+1,NCVAL
          CONT=ACVAL(I)
          IF(CONT.GT.BMIN.AND.CONT.LT.BMAX) THEN
            IF(acgrad(I).NE.GRAD0) THEN
              GRAD0=acgrad(I)
              GRAD=GRAD0*GRAD0
              CALL SETUP(GRID,FLAGS,WORK)
            END IF
            IF(JCDASH(1).NE.999) LINET=IABS(JCDASH(I))
            size=acsize(i)
            CALL SETLAB
            CALL SCAN(GRID,FLAGS,WORK,ierr)
c***********OUTPUT PERCENTAGE COMPLETED
            nout=nout+1
            call perout(nout,ncont,ipass,npass,'COMPLETED')
            if(ierr.ne.0) return
          END IF
    6 continue
        gradi = grdtmp
        GOTO 140
      END IF
C
C IF ACGRAD IS NOT BEING USED THEN DO SETUP ONLY ONCE
C
      if (nsec.le.0.or.ncval.gt.0) then
      grad=abs(gradi)**2
      else
      grad=(abs(gradi)*dcval)**2
      endif
c
c  set flags, etc..
      call setup(grid,flags,work)
10    if (bmin.ge.bmax) go to 140
      if (ncval.eq.0) go to 50
c
c  find lower limit.
      do 20 i=1,ncval
      if (acval(i).le.bmin) go to 20
      ii=i
      go to 30
   20 continue
c
c  contouring loop
c
   30 prime=.true.
c*****FIND UPPER LIMIT
      do 35 i=ii,ncval
       if(acval(i).gt.bmax) go to 37
   35 continue
      i=ncval+1
   37 ncmax=i-1

      nout=0
      ncont = ncmax-ii + 1

      do 40 i=ii,ncmax
      cont=acval(i)
c      if (cont.gt.bmax) go to 140
      IF(JCDASH(1).NE.999) LINET=IABS(JCDASH(I))
      call setlab
      call scan(grid,flags,work,ierr)
      if(ierr.ne.0) return
c*******OUTPUT PERCENTAGE COMPLETED
      nout=nout+1
      call perout(nout,ncont,ipass,npass,'COMPLETED')
   40 continue
      go to 140
c
c  execution for delta contour levels.
   50 if (cmin.eq.-fltmax.and.cmax.eq.fltmax) go to 60
      temin=amax1(bmin,cmin)/dcval
      icont=int(temin)
      if(aint(temin).ne.temin.and.temin.gt.0) icont=icont+1
c      icont=amax1(bmin,cmin)/dcval
      amx=amin1(bmax,cmax)
      go to 65
c   60 icont=bmin/dcval
   60 temin=bmin/dcval
      icont=int(temin)
      if(aint(temin).ne.temin.and.temin.gt.0) icont=icont+1
c     icont=bmin/dcval
      amx=bmax

   65 cont=icont*dcval
      nout=0
      ncont = ifix((amx-cont)/dcval) + 1
c
c  contouring loop.
   70 cont=icont*dcval
      if (cont.gt.amx) go to 140
      prime=nsec.gt.0.and.mod(icont,nsec).eq.0
      if (prime) go to 100
      if (idashs) 80,135,90
c  decode primary contour line thickness
   80 jcol=(-idashs/100)*100
      linet= jcol + mod((-idashs-jcol),8)
c  80 linet=iabs(mod(idashs,8))
      go to 135
   90 linet=0 + (idashs/100) * 100
      go to 135
  100 if (idashs) 110,130,120
  110 jcol=(-idashs/100)*100
      linet=iabs((idashs+jcol)/8)*8
c 110 linet=iabs(idashs/8)*8
      go to 130
  120 linet=idashs
  130 call setlab
C
  135 call scan(grid,flags,work,ierr)
      if(ierr.ne.0) return
c*******OUTPUT PERCENTAGE COMPLETED
      nout=nout+1
      call perout(nout,ncont,ipass,npass,'COMPLETED')
      icont=icont+1
      go to 70
c
c  done with block
  140 return
      end
C________________________________________________________________
C
C     FUNCTION  E R R B A R
C________________________________________________________________
C
C FUNCTION ERROR BAR FINDS THE DISTANCE ALONG AN INTERVAL OF
C LENGTH B THAT A DEPARTS FROM THE NEAREST MULTIPLE OF B.
C
C   A      - REAL*4.  TEST NUMBER
C   B      - REAL*4.  LENGTH OF INTERVAL.
C
C            EXAMPLES:  IF A= 21 AND B= 20 THEN ERRBAR= .05
C                       IF A=101 AND B= 20 THEN ERRBAR= .05
C                       IF A= 49 AND B= 10 THEN ERRBAR=-.1
C
C FUNCTION ERRBAR WRITTEN BY ROB BRACKEN, USGS, 2MAR88.
C
C
      function errbar(a,b)
C
      real*8 aob
C
      aob=dble(a)/dble(b)
      errbar=aob-dnint(aob)
C
      return
      end
C________________________________________________________________
C
C     SUBROUTINE  S E T U P
C________________________________________________________________
C
      subroutine setup(grid,flags,work)
      dimension grid(0:*),flags(0:*),work(1)
      integer flags
      common /concom/ ncol,nrow,bmin,bmax,
     & grad,gflg,ijs,ije,ijsi,ijei,prime,ij4,cont,
     & fltmax,lmult(0:3),idashs,linet,sigma,nsig
      logical prime
c initialization of all parameters (arithmetic and logical)
c  unique to row block being contoured
c
c system of flags
c     bit         definition (for bit on)
c     0       edge-side 4 possible cut
c     1       edge-side 3 possible cut
c     2       edge-side 2 possible cut
c     3       edge-side 1 possible cut
c     4       edge 1 cut not made and not checked
c     5       interior scan (edge 1)
c     6       gradient (drop secondary contours)
c     7       block good (contourable)
c
      integer fg
      logical ccc,collrc,rowlg,bit,kk,collg
        DATA IDVAL/#7FFFFFFF/
c
c start of operations
c
      nrow1=nrow-1
      ncol1=ncol-1
      ij=0
      i1j=1
      ij1=ncol
      i1j1=ij1+1
      ijn1=-ncol
      ijs=(ncol*nrow+3)/4-1
      do 10 i=0,ijs
   10 flags(i)=0
      ijs=idval
      ije=0
c
c preset bmin,bmax
c
      bmin=fltmax
      bmax=-bmin
c
c start of setup scanning
c
      do 130 j=0,nrow1
      rowlg=.false.
      collrc=.false.
      if (grad.eq.0.or.j.eq.nrow1) go to 20
      call coord(work,0,j,xll,yll)
      call coord(work,0,j+1,xul,yul)
   20 continue
      do 120 i=0,ncol1
      if (j.eq.0) go to 30
      collg=bit(flags,ijn1,1)
      go to 40
   30 collg=.false.
   40 if (i.eq.ncol1) go to 70
      ccc=collg
      if (j.eq.nrow1) go to 70
c
c determine if mesh block flagged
c
      if (gflg.ne.0.and.
     & (grid(ij).gt.gflg.or.grid(i1j).gt.gflg.or.
     &  grid(ij1).gt.gflg.or.grid(i1j1).gt.gflg)) go to 70
c
c mesh block contourable
c
      fg=1
      ijs=min0(ijs,ij)
      ije=max0(ije,ij)
c
c check and set left edge (2)
c
      if (rowlg) go to 50
      if (grid(ij1).gt.grid(ij)) fg=fg+32
      rowlg=.true.
c
c check and set lower edge (1)
c
   50 l=0
      if (grid(ij).gt.grid(i1j)) l=4
      if (.not.collg) l=l*4
      fg=fg+l
c
c check and set gradient
c
      if (grad.eq.0.) go to 60
      call coord(work,i+1,j,xlr,ylr)
      call coord(work,i+1,j+1,xur,yur)
      if ((grid(ij)-grid(i1j))**2.gt.
     &      ((xll-xlr)**2+(yll-ylr)**2)*grad
     & .or.(grid(ij)-grid(ij1))**2.gt.
     &      ((xll-xul)**2+(yll-yul)**2)*grad) go to 55
      if((grid(ij1)-grid(i1j1))**2.gt.
     &      ((xul-xur)**2+(yul-yur)**2)*grad
     & .or.(grid(ij1)-grid(i1j1))**2.gt.
     &      ((xur-xlr)**2+(yur-ylr)**2)*grad) go to 55
      if((grid(ij)-grid(i1j1))**2.lt.
     &      ((xll-xur)**2+(yll-yur)**2)*grad
     & .and.(grid(i1j)-grid(ij1))**2.lt.
     &      ((xul-xlr)**2+(yul-ylr)**2)*grad) go to 57
   55 fg=fg+2
   57 xll=xlr
      yll=ylr
      xul=xur
      yul=yur
c
c determine bmin,bmax
c
   60 bmin=amin1(bmin,grid(ij))
      bmax=amax1(bmax,grid(ij))
      call set(flags,ij,fg)
      go to 110
c
c mesh block non-contourable
c
c check and set side 3 and 4 flags
c
   70 kk=.false.
      if (.not.rowlg) go to 80
      if (grid(ij).gt.grid(ij1)) call set(flags,ij-1,128)
      kk=.true.
   80 if (.not.collg) go to 90
      if (grid(ij).lt.grid(i1j)) call set(flags,ijn1,64)
      kk=.true.
   90 if (.not.(collrc.or.kk)) go to 100
      bmin=amin1(bmin,grid(ij))
      bmax=amax1(bmax,grid(ij))
  100 rowlg=.false.
      fg=0
c
c end of mesh
c
  110 ij=i1j
      i1j=i1j+1
      ij1=i1j1
      i1j1=i1j1+1
      collrc=ccc
      ijn1=ijn1+1
  120 continue
c
  130 continue
c
      ijs=ijs/4
      ije=ije/4
      ijsi=(ijs+ncol)/4
      ijei=ije
      return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  S C A N
C________________________________________________________________
C
      subroutine scan(grid,flags,work,ierr)
c
c  scans grid for undrafted contour level.
c
      dimension grid(0:*),flags(0:*),work(1)
      integer flags
      common /concom/ ncol,nrow,bmin,bmax,
     & grad,gflg,ijs,ije,ijsi,ijei,prime,ij4,cont,
     & fltmax,lmult(0:3),idashs,linet,sigma,nsig
      LOGICAL PRIME,EDGE,bit
C
C       CHANGED DATA STATEMENTS TO BE COMPATIBLE WITH THE
C    3   2 BIT VAX.  C. COOPER, OCTOBER 7, 1981.
C       THESE STATEMENTS ARE FOR THE 36 BIT HONEYWELL/MULTICS
C     data nask1/o4004004004/,nask2/o360360360360/,nask3/o360/
C    1  ,NASK4/O20/
C     data nask5/o010010010010/,nask6/o010/,nask7/o4/,
C    & nask8/o200/,nask9/o100/,nask10/o40/,nask11/o004004004004/
      DATA MASK1/#04040404/,MASK2/#F0F0F0F0/,MASK3/#F0/,
     1     MASK4/#10/
      DATA MASK5/#08080808/,MASK6/#08/,MASK7/#4/,
     # MASK8/#80/,MASK9/#40/,MASK10/#20/,MASK11/#04040404/
      data izero/0/,mask12/#18/
c
c  set interior flags.
      if (ijsi.lt.0) go to 20
      do 10 ij=ijsi,ijei
   10 flags(ij)=ior(flags(ij),ishft(iand(flags(ij),mask1),1))
c
c  setup for edge scan
   20 if (ijs.lt.0) go to 180
      assign 40 to iswa
      assign 160 to iswb
      edge=.true.
      ija=ijs
      ijb=ije
      ijsnew=-1
      mask=mask2
      maskb=mask3
      maskd=mask4
      maskf=1
   30 ij4b=ija*4+3
c
c  basic word scan
      do 170 ij=ija,ijb
      mflag=flags(ij)
      if (iand(mflag,mask).eq.izero) go to iswb,(160,150)
      if (ijsnew.lt.0) ijsnew=ij
      ijenew=ij
      ij4=ij4b
      ij4n=ij4+ncol
c
c  sub word scan
      do 140 k=1,4
      if (iand(mflag,maskb).eq.izero) go to 130
      go to iswa,(40,110)
c
c  right edge
   40 if (iand(mflag,mask8).eq.izero) go to 60
      if (cont.le.grid(ij4+1)) go to 50
      CALL CRESET(flags,ij4,mask8)
      go to 80
   50 if (cont.le.grid(ij4n+1)) go to 60
      call trace(grid,flags,work,4,maskf,ierr)
      if(ierr.ne.0) return
      go to 80
c
c  top edge
   60 if (iand(mflag,mask9).eq.izero) go to 80
      if (cont.le.grid(ij4n+1)) go to 70
      CALL CRESET(flags,ij4,mask9)
      go to 100
   70 if (cont.le.grid(ij4n)) go to 80
      call trace(grid,flags,work,3,maskf,ierr)
      if(ierr.ne.0) return
      go to 100
c
c  left edge
   80 if (iand(mflag,mask10).eq.izero) go to 100
      if (cont.le.grid(ij4n)) go to 90
      CALL CRESET(flags,ij4,mask10)
      go to 130
   90 if (cont.le.grid(ij4)) go to 100
      call trace(grid,flags,work,2,maskf,ierr)
      if(ierr.ne.0) return
      go to 130
c
c  bottom edge and interior scan
  100 if (iand(mflag,mask4).eq.izero) go to 130
  110 if (cont.le.grid(ij4)) go to 120
      CALL CRESET(flags,ij4,maskd)
      go to 130
c     mod by r.godson 6/89 after finding problem when
c     trace is reentered from subroutine scan in the
c     do 140 loop(sub word scan).A grid flag within
c     the sub word scan can be reset by the previous call
c     to trace and would not be caught until two
c     points are traced and therefore there is a
c     duplication of those two points.
c  120 if (cont.gt.grid(ij4+1))
c     & call trace(grid,flags,work,1,maskf)
  120 if (cont.gt.grid(ij4+1)) then
      if(bit(flags,ij4,mask12))
     & call trace(grid,flags,work,1,maskf,ierr)
      if(ierr.ne.0) return
      endif
c
C 130 mflag=mflag/512
130   mflag=ishft(mflag,-8)
      ij4=ij4-1
      ij4n=ij4n-1
  140 continue
      go to 160
c
c  interior range set
  150 if (iand(mflag,mask11).eq.izero) go to 160
      if (ijsnew.lt.0) ijsnew=ij
      ijenew=ij
  160 ij4b=ij4b+4
  170 continue
c
c  first time?
      if (.not.edge) go to 190
c
c  update edge range
      ijs=ijsnew
      ije=ijenew
  180 if (ijsi.lt.0) go to 200
c
c  setup for interior scan
      edge=.false.
      assign 110 to iswa
      assign 150 to iswb
      ijsnew=-1
      ija=ijsi
      ijb=ijei
      mask=mask5
      maskb=mask6
      maskd=mask7
      maskf=maskb
      go to 30
c
c  update interior boundary indicies
  190 ijsi=ijsnew
      ijei=ijenew
c
c  end of scan
  200 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  T R A C E
C________________________________________________________________
C
      subroutine trace(grid,flags,work,iside,mask,ierr)
c
c  follows contour through grid until
c  edge or closure found.
c
      dimension grid(0:*),flags(0:*),work(1)
      dimension istart(100),iend(100),ic(100)
      integer flags,ITURN(4,4)
      common/chksiz/conlim
      common /contrc/ cmin,cmax,dcval,ncval,nsec,gradi,
     1                ACDEL(200),ACGRAD(200),JCDASH(200),
     2                ncval2,acsize(200),acmin(200),acmax(200)
      common /concom/ ncol,nrow,bmin,bmax,
     & grad,gflg,ijs,ije,ijsi,ijei,prime,ij4,cont,
     & fltmax,lmult(0:3),idashs,linet,sigma,nsig
c
      common/labcom/ dum1(11),size,dum2(4)
      logical prime,bit,post,ksw,break,PRIME2
c
c trace and plot contour through grid
c
c
c grid indexing
c                  side 3
c          i01 +           + i11
c
c      side 2                 side 4
c
c          i00 +           + i10
c                  side 1
c
c in-- entrance side
c high-equal point(inn) always on left when looking along
c     contour line
c iol-- point opposite entrance side on left
c ior-- point opposite entrance side on right
c
      dimension x(101),y(101)
c     data mask6/o10/
      data mask6/#8/
      DATA ITURN /-2,1,0,-1,-1,-2,1,0,0,-1,2,1,1,0,-1,2/
c
c set up tracing start
c
      kk=1
      ksw=.false.
      break=.false.
      istart(1)=0
      NTURN=0
      post=prime
      i=mod(ij4,ncol)
      j=ij4/ncol
      npts=1
      ic(1)=0
      i00=ij4
      i10=i00+1
      i01=i00+ncol
      i11=i01+1
      go to (10,20,30,40),iside
   10 CALL CRESET(flags,i00,mask6)
      fract=(cont-grid(i00))/(grid(i10)-grid(i00))
      call coord(work,i,j,xa,ya)
      call coord(work,i+1,j,xb,yb)
      iol=i01
      ior=i11
      in=1
      go to 51
   20 fract=(cont-grid(i00))/(grid(i01)-grid(i00))
      call coord(work,i,j,xa,ya)
      call coord(work,i,j+1,xb,yb)
      iol=i11
      ior=i10
      in=2
      go to 51
   30 fract=(cont-grid(i01))/(grid(i11)-grid(i01))
      call coord(work,i,j+1,xa,ya)
      call coord(work,i+1,j+1,xb,yb)
      iol=i10
      ior=i00
      in=3
      go to 51
   40 fract=(cont-grid(i10))/(grid(i11)-grid(i10))
      call coord(work,i+1,j,xa,ya)
      call coord(work,i+1,j+1,xb,yb)
      iol=i00
      ior=i01
      in=4
c
c start trace loop
c
   51 xbeg=xa+(xb-xa)*fract
      ybeg=ya+(yb-ya)*fract
      inbeg=in
C
C INITIALIZE VARIABLES TO MEASURE LENGTH OF TRACE
C
      TLEN=0.
      X0=XBEG
      Y0=YBEG
C
   50 X1=xa+(xb-xa)*fract
      Y1=ya+(yb-ya)*fract
      X(NPTS)=X1
      Y(NPTS)=Y1
C
C CALCULATE TRACE LENGTH
C
      DX=X1-X0
      DY=Y1-Y0
      DS=SQRT(DX*DX+DY*DY)
      X0=X1
      Y0=Y1
      TLEN=TLEN+DS
c
c     if gradi is negative bypass check for primaries
c
      if(gradi.ge.0.) then
        if (.not.prime.and.bit(flags,i00,2)) go to 80
      else
        if(bit(flags,i00,2)) go to 80
      endif
c
c     get starting point for line segment kk
c
      if(.not.ksw) then
        istart(kk)=npts
        ksw=.true.
      endif
      go to 85
c
c     get ending point for line segment kk
c
   80 if(ksw) then
        iend(kk)=npts
        kk=kk+1
        ic(kk)=0
        istart(kk)=0
        ksw=.false.
      endif
      break=.true.
c
c     accumulate points whether lines are broken or not
c
   85 npts=npts+1
      if (npts.le.100) go to 90
      call conplt(x,y,npts,ic,post,1,ksw,
     & break,istart,iend,kk)
      if(ksw) ic(1)=2
c      if(ksw.and.(istart(kk).eq.100)) ic(1)=2
   87 kk=1
      npts=1
      ksw=.false.
      break=.false.
      istart(1)=0
c
c determine exit side
c
   90 k=in+2
      if (grid(ior).lt.cont) k=k-1
      if (grid(iol).lt.cont) k=k-2
      if (k-in) 110,100,120
c
c saddle decision (dayhoff)
c
  100 if ((grid(i00)+grid(i10)+grid(i01)+grid(i11))*0.25.lt.cont)
     & go to 120
      k=in+2
      go to 120
  110 k=in
c
c     compute side branch
c
  120 go to (140,150,160,130,140,150),k
c
c exit bottom -- side 1
c
  130 fract=(cont-grid(i00))/(grid(i10)-grid(i00))
      call coord(work,i,j,xa,ya)
      call coord(work,i+1,j,xb,yb)
      NTURN=NTURN+ITURN(IN,1)
      if (j.eq.0) THEN
        IOUT=1
        go to 170
      END IF
      i01=i00
      i00=i00-ncol
      if (.not.bit(flags,i00,1)) THEN
        IOUT=1
        go to 170
      END IF
      in=3
      i11=i10
      i10=i10-ncol
      iol=i10
      ior=i00
      j=j-1
      go to 50
c
c exit left -- side 2
c
  140 fract=(cont-grid(i00))/(grid(i01)-grid(i00))
      call coord(work,i,j,xa,ya)
      call coord(work,i,j+1,xb,yb)
      NTURN=NTURN+ITURN(IN,2)
      if (i.eq.0) THEN
        IOUT=2
        go to 170
      END IF
      i10=i00
      i00=i00-1
      if (.not.bit(flags,i00,1)) THEN
        IOUT=2
        go to 170
      END IF
      in=4
      i11=i01
      i01=i01-1
      iol=i00
      ior=i01
      i=i-1
      go to 50
c
c exit top -- side 3
c
  150 fract=(cont-grid(i01))/(grid(i11)-grid(i01))
      call coord(work,i,j+1,xa,ya)
      call coord(work,i+1,j+1,xb,yb)
      NTURN=NTURN+ITURN(IN,3)
      if (.not.bit(flags,i01,mask)) THEN
        IOUT=3
        go to 170
      END IF
      CALL CRESET(flags,i01,mask6)
      in=1
      i00=i01
      i10=i11
      i01=i01+ncol
      i11=i11+ncol
      iol=i01
      ior=i11
      j=j+1
      go to 50
c
c exit right -- side 4
c
  160 fract=(cont-grid(i10))/(grid(i11)-grid(i10))
      call coord(work,i+1,j,xa,ya)
      call coord(work,i+1,j+1,xb,yb)
      NTURN=NTURN+ITURN(IN,4)
      if (.not.bit(flags,i10,1)) THEN
        IOUT=4
        go to 170
      END IF
      in=2
      i00=i10
      i01=i11
      i10=i10+1
      i11=i11+1
      iol=i11
      ior=i10
      i=i+1
      go to 50
c
c end of trace loop
c
  170 x(npts)=xa+(xb-xa)*fract
      y(npts)=ya+(yb-ya)*fract
      TLEN=TLEN+SQRT((X(NPTS)-X0)*(X(NPTS)-X0)+
     1     (Y(NPTS)-Y0)*(Y(NPTS)-Y0))
      if (npts.le.1.and.ic(1).eq.0) go to 190
      call conplt(x,y,npts,ic,post,3,ksw,
     & break,istart,iend,kk)
  190 XEND=X(NPTS)
      YEND=Y(NPTS)
      IF(GRADI.GE.0.) THEN
        PRIME2=PRIME
      ELSE
        PRIME2=.TRUE.
      END IF
      CALL LAYEGG(NTURN,XBEG,XEND,YBEG,YEND,TLEN,CONT,linet,
     1 PRIME2,size,ierr)
      RETURN
C  190 return
c
c end of tracing
c
      end
C
C________________________________________________________________
C
C     FUNCTION  B I T
C________________________________________________________________
C
      LOGICAL FUNCTION BIT(LF,I,MASK)
C     The byte is ANDed with MASK.  If any bits of the result
C     are set, the result is .TRUE., otherwise the result is
C     .FALSE.
      DIMENSION LF(0:*)
      J=I/4
      INDEX=IAND(I,3)
      K=ISHFT(MASK,8*(3-INDEX))
      BIT=IAND(LF(J),K).NE.0
      END
C
C________________________________________________________________
C
C     FUNCTION  I T H C H A R
C________________________________________________________________
C
      FUNCTION ITHCHAR(STRING,I)
      CHARACTER*(*) STRING
      ITHCHAR=ICHAR(STRING(I:I))
      END
C
C________________________________________________________________
C
C     SUBROUTINE  S E T
C________________________________________________________________
C
      SUBROUTINE SET(LF,I,MASK)
C     Or's MASK with the Ith byte stored in LF, where I is in
C     the MULTICS HONEYWELL order:
C   MULTICS: (3)(2)(1)(0)(7)(6)(5)(4)(11)(10)( 9) . . .
C       VAX: (0)(1)(2)(3)(4)(5)(6)(7) (8)(9)(10) . . .
      DIMENSION LF(0:*)
      J=I/4
      INDEX=IAND(I,3)
      K=ISHFT(MASK,8*(3-INDEX))
      LF(J)=IOR(LF(J),K)
      END
C
C________________________________________________________________
C
C     SUBROUTINE  C R E S E T
C________________________________________________________________
C
      SUBROUTINE CRESET(LF,I,MASK)
C     Subroutine to reset logical bits -- provides VAX calls
C     compatible with Perkin-Elmer calls.
C     Does appropriate subscript manipulation for
C     character-string
C     ordering:
C       P-E: (3)(2)(1)(0)(7)(6)(5)(4)(11)(10) . . .
C       VAX: (0)(1)(2)(3)(4)(5)(6)(7) (8)(9)(10) . . .
C
      DIMENSION LF(0:*)
      J=I/4
      INDEX=IAND(I,3)
      K=ISHFT(MASK,8*(3-INDEX))
      LF(J)=IAND(LF(J),NOT(K))
      END
C
C________________________________________________________________
C
C     SUBROUTINE  C O O R D A
C________________________________________________________________
C
      subroutine coorda(ixad,iyad,ncol,iquad)
      dimension x(1)
      k=ixad
      l=iyad
      ncoll=ncol
      isw=iquad
      return
C
C________________________________________________________________
C
C     ENTRY  C O O R D
C________________________________________________________________
C
      entry coord(x,i,j,xv,yv)
      go to (10,10,20),isw
   10 xv=x(k+i)
      yv=x(l+j)
      return
   20 ll=j*ncoll+i
      xv=x(k+ll)
      yv=x(l+ll)
      return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  C O N P L T
C________________________________________________________________
C
      subroutine conplt(x,y,npts,ic,post,isw,ksw,
     & break,istart,iend,kk)
C
      common/concom/ncol,nrow,bmin,bmax,grad,gflg,ijs,ije,
     & ijsi,ijei,prime,ij4,cont,fltmax,lmult(0:3),idashs,linet,
     & sigma,nsig
      COMMON /eggchu/ GRADEA(4),NOHEN
      dimension x(101),y(101),xx(1001),yy(1001)
      dimension istart(1),iend(1),ic(1)
      logical post,prime,ksw,break,NOHEN
c
      icc=ic(1)
      if(npts.le.1) go to 90
      if(isw.eq.1) npts=100
      if(nsig.eq.0.and..not.break) go to 110
      if(.not.break) then
c
c       spline but no breaks in segment
c
        call bsplin(x,y,npts,xx,yy,nnt,ninc,sigma)
        go to 56
      else
        if(ksw) iend(kk)=npts
        if(.not.ksw.and.istart(kk).eq.0) kk=kk-1
        if(nsig.ne.0) then
c
c         spline and breaks in segment
c
          call bsplin(x,y,npts,xx,yy,nnt,ninc,sigma)
          call feed(xx,yy,nnt,iccc,-999,100+kk)
          if(kk.eq.0) go to 50
          if(ic(1).ne.2) go to 20
          if(istart(1).ne.1) then
            ic(1)=0
          else
            ic(1)=1
          endif
   20     do 30 i=1,kk
          is=(istart(i)-1)*ninc+1
          ie=(iend(i)-1)*ninc+1
          nnnt=ie-is+1
          icc=ic(i)
C          call line(xx(is),yy(is),nnnt,icc,linet)
          call feed(xx(is),yy(is),nnnt,icc,LINET,isw)
   30     continue
        else
c
c         no spline but breaks in segment
c
          call feed(x,y,npts,iccc,-999,100+kk)
          if(kk.eq.0) go to 50
          if(ic(1).ne.2) go to 35
          if(istart(1).ne.1) then
            ic(1)=0
          else
            ic(1)=1
          endif
   35     do 40 i=1,kk
          is=istart(i)
          ie=iend(i)
          nnnt=ie-is+1
          icc=ic(i)
C          call line(x(is),y(is),nnnt,icc,linet)
                      call feed(x(is),y(is),nnnt,icc,LINET,isw)
   40     continue
        endif
      endif
   50 ic(1)=0
      return
c
c     after splining and no breaks
c
   56 go to (57,80,95),isw
   57 if (post) go to 60
C      call line(xx,yy,nnt,icc,linet)
      call FEED(xx,yy,nnt,icc,LINET,ISW)
      post=prime
      go to 70
   60 IF(NOHEN) THEN
        call label(xx,yy,nnt,icc,post,ISW)
      ELSE
        call feed(xx,yy,nnt,icc,linet,ISW)
      END IF
   70 npts=1
      ic(1)=1
      return
C   80 call line(xx,yy,nnt,icc,linet)
   80 call FEED(xx,yy,nnt,icc,LINET,ISW)
      ic(1)=0
      npts=1
      return
   90 xx(1)=x(1)
      yy(1)=y(1)
      nnt=npts
   95 if(post) go to 100
C      call line(xx,yy,nnt,icc,linet)
      call FEED(xx,yy,nnt,icc,LINET,ISW)
      return
  100 IF(NOHEN) THEN
       call label(xx,yy,nnt,icc,post,ISW)
      ELSE
       call feed(xx,yy,nnt,icc,linet,ISW)
      END IF
      return
c
c     no splining and no breaks
c
  110 go to(120,150,160),isw
  120 if(post) go to 130
C      call line(x,y,npts,icc,linet)
      call FEED(x,y,npts,icc,LINET,ISW)
      post=prime
      go to 140
  130 IF(NOHEN) THEN
       call label(x,y,npts,icc,post,ISW)
      ELSE
       call feed(x,y,npts,icc,linet,ISW)
      END IF
  140 npts=1
      ic(1)=1
      return
C  150 call line(x,y,npts,icc,linet)
  150 call FEED(x,y,npts,icc,LINET,ISW)
      npts=1
      ic(1)=0
      return
  160 if(post) go to 170
C      call line(x,y,npts,icc,linet)
      call FEED(x,y,npts,icc,LINET,ISW)
      return
  170 IF(NOHEN) THEN
       call label(x,y,npts,icc,post,ISW)
      ELSE
       call feed(x,y,npts,icc,linet,ISW)
      END IF
      return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  L A B E L
C________________________________________________________________
C
      subroutine label(x,y,npts,icc,post,ISW)
c
c  label scans contour line array looking for
c  a straight segment.
c  if found,  the contour is labelled.
c
      common /concom/ ncol,nrow,bmin,bmax,
     & grad,gflg,ijs,ije,ijsi,ijei,prime,ij4,cont,
     & fltmax,lmult(0:3),idashs,linet,sigma,nsig
      common /setxy/ xx(2),yy(2),xp1(4),yp1(4),iplotr,sizel,
     & ncharx,nchary,sizex,sizey,adelx,adely,pllx,plly,
     & lintx,linty,xscale,yscale,xxscal,yyscal,mscale,delx,dely
      common/labcom/ lchars(6),fmtc,nchar,size,wdist,wdist2,ccor,
     1 delb
      common/pltcon/npens
c
      dimension x(1),y(1),xxa(1),yya(1)
      character fmtc*16,chars*24
      logical nolab,post,prime
      equivalence (f,sum),(isa,ya),(lchars,chars)
c
c setup and scan
c
      nolab=.not.(nchar.gt.0.and.size.gt.0.)
      if (nolab) go to 50
      is=2
      isa=3
      sum=0.0
      in1=2
      ie=npts-1
      if (ie.lt.3) go to 50
      xb=x(is)
      yb=y(is)
      do 40 i=isa,ie
      sum=sum+sqrt((x(i)-x(in1))**2+(y(i)-y(in1))**2)
   10 dist=(x(i)-xb)**2 + (y(i)-yb)**2
      if (dist.lt.wdist2) go to 30
      dist=sqrt(dist)
      if ((sum/dist).gt.1.02) go to 20
      ie=i
      go to 60
   20 is=is+1
      sum=sum-sqrt((x(is)-xb)**2+(y(is)-yb)**2)
      xb=x(is)
      yb=y(is)
      if (is.lt.i) go to 10
   30 in1=i
   40 continue
c
c  can't fine spot, continue line without labelling.
C   50 call line(x,y,npts,icc,linet)
   50 call FEED(x,y,npts,icc,LINET,ISW)
      go to 70
c
c  plottable location.
c
C   60 call line(x,y,is,icc,linet)
   60 call FEED(x,y,is,icc,LINET,ISW+10)
      call FEED(x(is+1),y(is+1),ie-is,1,-999,ISW+40)
      xa=x(ie)-x(in1)
      ya=y(ie)-y(in1)
      f=1.-(dist-wdist)/sqrt(xa*xa+ya*ya)
      xa=f*xa + x(in1)
      ya=f*ya + y(in1)
      phi=atan2(ya-yb,xa-xb)
      if (abs(phi).gt.1.5707963) phi=phi-sign(3.1415927,phi)
C
C     WRITE LABELS OUT TO A FILE SO THE MAY BE DRAWN ALL
C     AT ONCE IF THE HP7580A IS USED.  THEN CALL VCHAR IN
C     THE MAIN PROGRAM TO READ THE FILE.
      IF (IPLOTR .EQ. 5) THEN
               ipen=linet/100
               if(ipen.lt.1) ipen=1
               if(ipen.gt.npens) ipen=npens
               call newpen(ipen)
               call VCHAR(0.5*(xa+xb),0.5*(ya+yb),lchars,nchar,
     1                     2,size,phi,ccor,0.)
      ELSE
                ilinet = sign(2 +(linet-mod(linet,100)),linet)
                call VCHAR(0.5*(xa+xb),0.5*(ya+yb),lchars,nchar,
     1                     ilinet,size,phi,ccor,0.)
      ENDIF
      xxa(1)=xa
      yya(1)=ya
C      call line(xxa(1),yya(1),1,0,linet)
      call FEED(xxa(1),yya(1),1,0,LINET,ISW+30)
C      call line(x(ie),y(ie),npts-ie+1,1,linet)
      call FEED(x(ie),y(ie),npts-ie+1,1,LINET,ISW+20)
      post=.false.
c
   70 return
C
C________________________________________________________________
C
C     ENTRY  S E T L A B
C________________________________________________________________
C
c
c  entry to establish labeling character string.
c  must be called before label.
c
      entry setlab
      nolab=.not.(nchar.gt.0.and.size.gt.0.)
      if (nolab) go to 100
C
C ADDITION MADE 20FEB87 TO ALLOW VARIABLE CONTOUR LINE LABEL
C FORMATS
C
      CALL FMTMOD(CONT,NCHAR,FMTC)
C
      write(chars,fmtc) cont
      wdist=(nchar+1)*size-.4*SIZE
      wdist2=wdist*wdist
      ccor=-.5*wdist+size
  100 return
      end
C
C________________________________________________________________
C
C       SUBROUTINE  F M T M O D
C________________________________________________________________
C
C SUBROUTINE FORMAT MODIFY SETS UP THE NUMBER OF CHARACTERS
C (NCHAR) AND THE MINIMUM SPACE FORMAT (FMT) FOR THE VALUE PASSED
C (CONT).  THE FIRST TIME FMTMOD IS CALLED IT UTILIZES THE
C INFORMATION IN NCHAR AND FMT TO DETERMINE THE ALLOWABLE NUMBER
C OF DECIMAL PLACES TO THE RIGHT OF THE DECIMAL POINT AND THE
C NUMBER OF TO CUT FROM THE RIGHT.  IN SUBSEQUENT CALLS NCHAR
C AND FMT ARE DETERMINED BASED ON THE CONSTANTS PICKED BY THE
C FIRST CALL.  THEREFORE IT IS ADVISABLE THAT THIS SUBROUTINE
C BE USED FOR ONLY ONE STREAM OF INPUT VALUES.
C
C SUBROUTINE FMTMOD WRITTEN BY ROB BRACKEN, USGS, 19FEB87
C
C
      subroutine fmtmod(cont,nchar,fmt)
      save
C
C VARIABLE DECLARATION
C
      character*(*) fmt
      character*16 fmt2
      real*8 cont2,tendec,fcont2
      data nchar2 /-999/
C
C DETERMINE IF LENDEC AND LENCUT HAVE BEEN FOUND
C
      if(nchar2.ne.-999) goto 30
C
C CONVERT THE FORMAT LENGTH SPECIFIER TO AN INTEGER
C
      fmt2=fmt
      if(fmt2(2:2).ne.'f'.and.fmt2(2:2).ne.'F') return
      do 5 i=3,16
        if(fmt2(i:i).eq.'.') goto 10
    5 continue
      return
   10 if(i.eq.3) return
      read(fmt2(3:i-1),1000) lenfmt
 1000 format(i2)
C
C CONVERT THE FORMAT DECIMAL SPECIFIER TO AN INTEGER
C
      do 15 j=i+1,16
        if(fmt2(j:j).eq.')'.or.fmt2(j:j).eq.' ') goto 20
   15 continue
      return
   20 if(j.eq.i+1) return
      read(fmt2(i+1:j-1),1000) lendec
C
C SET NCHAR2 TO INDICATE THAT LENDEC AND LENCUT HAVE BEEN FOUND
C
      nchar2=nchar
      lencut=lenfmt-nchar2
C
C ROUND CONT TO THE NUMBER OF DECIMAL PLACES STATED IN THE FORMAT
C
   30 cont2=dble(cont)
      tendec=10.d0**lendec
      fcont2=dnint(cont2*tendec)/tendec
C
C DETERMINE THE NUMBER OF DIGITS TO THE LEFT OF THE DECIMAL
C
      if(dabs(fcont2).ge.1.d0) then
        nleft=1+idnint(dlog10(dint(dabs(fcont2))))
      else
        nleft=1
      end if
      if(fcont2.lt.0.d0) nleft=nleft+1
C
C SET UP THE NEW NCHAR FOR THE PRESENT CONT
C
      lenfmt=nleft+1+lendec
      nchar=lenfmt-lencut
      if(dabs(fcont2).lt.1.d0.and.nchar.lt.1) nchar=1
C
C SET UP THE NEW FMT FOR THE PRESENT CONT
C
      llf=1
      lld=1
      if(lenfmt.ge.10) llf=2
      if(lendec.ge.10) lld=2
      if(llf.eq.1.and.lld.eq.1)write(fmt,801) lenfmt,lendec
  801 format('(f',i1,'.',i1,')')
      if(llf.eq.1.and.lld.eq.2)write(fmt,802) lenfmt,lendec
  802 format('(f',i1,'.',i2,')')
      if(llf.eq.2.and.lld.eq.1)write(fmt,803) lenfmt,lendec
  803 format('(f',i2,'.',i1,')')
      if(llf.eq.2.and.lld.eq.2)write(fmt,804) lenfmt,lendec
  804 format('(f',i2,'.',i2,')')
      return
      end
C
C_______________________________________________________________
C
C     SUBROUTINE  B S P L I N E
C_______________________________________________________________
C
      subroutine bsplin(x,y,n,xx,yy,nn,ninc,sigmap)
      save
c
c     b-spline subroutine
c     r.godson,usgs,denver,co.
c
      dimension x(1),y(1),xx(1),yy(1),c1(10),c2(10),c3(10),c4(10)
      logical close
      data sa/.1666667/,sb/.6666667/,isw/0/
      if(isw.eq.1) go to 10
      call setbsp(sigmap,step)
      k=0
      s1=1.-step + .000001
      do 5 t=0.,s1,step
      k=k+1
      t1=t*.5
      t2=t*t
      t2a=t2*.5
      t3=t2*t
      t3a=t3*.5
      c4(k)=sa*t3
      c1(k)=-c4(k)+t2a-t1+sa
      c2(k)=t3a-t2+sb
      c3(k)=-t3a+t2a+t1+sa
    5 continue
      isw=1
   10 if(n.lt.3) then
      do 20 i=1,n
      xx(i)=x(i)
      yy(i)=y(i)
   20 continue
      nn=n
      ninc=1
      return
      endif
      ninc=k
      if(x(1).eq.x(n).and.y(1).eq.y(n)) then
        close=.true.
        x1=x(n-1)
        y1=y(n-1)
        x(n+1)=x(2)
        y(n+1)=y(2)
      else
        close=.false.
        x11=x(2)-x(1)
        x1=x(1)-x11
        y11=y(2)-y(1)
        y1=y(1)-y11
        xnn=x(n)-x(n-1)
        x(n+1)=x(n)+xnn
        ynn=y(n)-y(n-1)
        y(n+1)=y(n)+ynn
      end if
      nn=1
      xx1=x(1)
      xx2=x(2)
      xx3=x(3)
      yy1=y(1)
      yy2=y(2)
      yy3=y(3)
      do 25 j=1,k
      xx(nn)=c1(j)*x1+c2(j)*xx1+c3(j)*xx2+c4(j)*xx3
      yy(nn)=c1(j)*y1+c2(j)*yy1+c3(j)*yy2+c4(j)*yy3
      nn=nn+1
   25 continue
      do 30 i=2,n-1
      i1=i-1
      i2=i+1
      i3=i+2
      xxi1=x(i1)
      xxi=x(i)
      xxi2=x(i2)
      xxi3=x(i3)
      yyi1=y(i1)
      yyi=y(i)
      yyi2=y(i2)
      yyi3=y(i3)
      do 27 j=1,k
      xx(nn)=c1(j)*xxi1+c2(j)*xxi+c3(j)*xxi2+c4(j)*xxi3
      yy(nn)=c1(j)*yyi1+c2(j)*yyi+c3(j)*yyi2+c4(j)*yyi3
      nn=nn+1
   27 continue
   30 continue
      if(close) then
        xx(nn)=xx(1)
        yy(nn)=yy(1)
      else
        xx(nn)=x(n)
        yy(nn)=y(n)
      end if
      return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  S E T B S P
C________________________________________________________________
C
C SUBROUTINE SET B-SPLINE RETURNS THE STEP SIZE FOR THE B-SPLINE
C WHEN GIVEN THE GRID CELL SIZE AND SCALE.
C
C PASSED VARIABLES:
C
C   DELX   - REAL*4.  X SIZE OF GRID CELL IN DATA UNITS.
C   DELY   - REAL*4.  Y SIZE
C   XXSCAL - REAL*4.  X DIRECTION MAP_INCHES/DATA_UNITS.
C   YYSCAL - REAL*4.  Y DIRECTION
C   SIGMAB - REAL*4.  MAXIMUM ALLOWABLE AVERAGE DISTANCE BETWEEN
C            SPLINED POLYGON CURVE AND PERFECT CURVE. THE MINIMUM
C            GRID CELL AREA (IN MAP INCHES) THAT WILL BE SPLINED
C            IS:  8 * [ PI*SIGMAB/( PI-2 ) ]**2
C            (RECOMMEND SIGMAB = .0073 INCHES FOR VERSATEC PLOTS)
C   STEP   - REAL*4.  NORMALIZED STEP DISTANCE FOR INPUT TO
C            B-SPLINE SUBROUTINE.
C
C SUBROUTINE SET B-SPLINE WRITTEN BY ROB BRACKEN, USGS, 6APR87
C
C
      subroutine setbsp(sigmab,step)
C
      parameter(pi=3.141592654,pio2=pi/2.)
      common/setxy/dum(26),xxscal,yyscal,mscale,delx,dely
C
C --- DISTANCE FROM GRID CELL CORNER TO CENTER IN MAP INCHES ---
      r=sqrt(delx*delx*xxscal*xxscal+dely*dely*yyscal*yyscal)/2.
C
      do 10 n=1,9
C ---   DIFFERENCE IN AREA BETWEEN CIRCLE AND POLYGON ---
        area=pi*r*r*(1.-(sin(pio2/n)*n)/pio2)
C ---   DISTANCE AROUND CIRCLE
        s=2.*pi*r
C ---   AVERAGE DISTANCE BETWEEN CIRCLE AND POLYGON BOUNDARY ---
        aos=area/s
C ---   IF DISTANCE IS WITHIN TOLERANCE THEN FINISHED ---
        if(aos.le.sigmab) goto 101
   10 continue
C
  101 step=1./n
      return
      end
