C
C   SUBROUTINE  H A C H E R 
C________________________________________________________________
C
C SUBROUTINE HACHER IS DESIGNED TO PUT HACHERS ON LOW CLOSURES
C WHEN CALLED BY PROGRAM CONTOUR.  THIS SUBROUTINE WILL BRIDGE
C TIER BOUNDARIES.  AN INITIAL CALL TO THE SUBROUTINE ITSELF
C MUST BE MADE BEFORE ANY CONTOURING BEGINS IN ORDER TO
C INITIALIZE VARIABLES.  THEN A CALL TO SUBROUTINE FEED MUST
C BE MADE EACH TIME A CONTOUR LINE SEGMENT IS DRAWN.  A CALL
C TO ENTRY LAYEGG MUST BE MADE AT THE END OF EACH CONTOUR TRACE.
C ENTRY GRADEGG MUST BE CALLED AFTER EACH TIER IS COMPLETED AND
C ENTRY GATHEREGGS MUST BE CALLED AT THE END OF THE CONTOURING.
C FOLLOWING IS A LIST OF THE PASSED VARIABLES:
C
C   NPASS  - INTEGER*4.  THE NUMBER OF TIERS TO BE CONTOURED.
C   NROWT  - INTEGER*4.  THE NUMBER OF ROWS IN EACH TIER.
C   TY0    - REAL*4.  THE Y POSITION IN MAP INCHES OF THE FIRST
C            ROW IN THE GRID.
C   TDELY  - REAL*4.  THE DISTANCE BETWEEN ROWS IN MAP INCHES.
C   HLEN   - REAL*4.  THE LENGTH IN INCHES OF THE HACHER TICKS.
C   HSPACE - REAL*4.  THE DISTANCE IN INCHES BETWEEN THE HACHER
C            TICKS.  HACHERS WILL NOT BE MADE IF HSPACE = 0.
C   HGAP   - REAL*4.  THE RATIO (DIRECT DISTANCE BETWEEN CONTOUR
C            END POINTS) / (THE TOTAL LENGTH OF THE CONTOUR).
C            BY APROPRIATELY SPECIFYING THIS VARIABLE LOW
C            CLOSURES WHICH INTERSECT MAP EDGES AND DVALED AREAS
C            WILL BE HACHERED.  SUGGESTED VALUE IS .1
C   HMAX   - REAL*4.  THE MAXIMUM NUMBER OF HACHURED CONTOUR
C            LINES THAT MAY BE STACKED AROUND ANY GIVEN LOW.
C            THIS NUMBER WILL BE CONVERTED TO AN INTEGER.
C   HLIM   - REAL*4.  THE MINIMUM CIRCUMFERENCE OF A CONTOUR
C            LINE THAT WILL BE CONSIDERED PART OF THE STACK
C            AND THEREFORE SUBJECT TO THE HMAX LIMIT.
C   HVB    - REAL*4  HIGHS WHICH ARE SHORTER THAN HVB WILL BE
C            IGNORED IF THEY ARE INSIDE A LOW CLOSURE.
C   CLIM   - REAL*4.  CONTOUR LIMIT.  THE SHORTEST LENGTH (IN
C            MAP INCHES) CONTOUR LINE WHICH WILL BE PLOTTED.
C   DELB   - REAL*4.  THE DESIRED DISTANCE IN MAP INCHES BETWEEN
C            CONTOUR LINE LABELS.
C
C SUBROUTINE HACHERY WRITTEN BY ROB BRACKEN, USGS, 17MAR87
C
C
      subroutine hacher(npass,nrowt,ty0,tdely,hlen,hspace,hgap,
     1 hmax,hlim,hvb,clim,delb)
      save
C
C VARIABLE DECLARATION
C
      integer*4 lenxy,lental,lenclo,nparm
      parameter(lenxy=1001,lental=5000,lenptr=lental+40)
      parameter(lenclo=5000,nparm=5)
      parameter(tol=.000001)
      real*4 tal(2,lental),tal2(2,lental),xx(lenxy),yy(lenxy)
      real*4 xa(lenclo),ya(lenclo)
      integer*4 ital(2,lental),lpoint(lenptr)
      integer*2 istack(lental)
      integer*4 jzero(13),ikey(10),iparm(nparm)
      integer*4 dsk1,dsk2,dsk3,dsk4,dsk5,south,north
      parameter(dsk1=14,dsk2=15,dsk3=16,dsk4=17,dsk5=18)
      equivalence (tal(1,1),ital(1,1))
      logical*4 nohach,skip,noclim,nodelb,nohen
      common /eggchu/ jsrec,nohach,noclim,nodelb,nohen
      common/captur/icapt
      data jzero /0,0,0,0,0,0,0,0,0,0,0,0,0/
c
c     open file for capturing contours if requested
c
      if(icapt.ne.0) then
      open(22,file='contour.cor',status='unknown',
     & form='formatted')
      open(23,file='contour.val',status='unknown',
     & form='formatted')
      endif
C
C DETERMINE IF HACHERS ARE GOING TO BE MADE
C
  100 nohach=.false.
      if(hlen.lt.1.e-30) nohach=.true.
      if(hspace.lt.1.e-30) nohach=.true.
      if(nint(hmax).eq.0) nohach=.true.
C
C DETERMINE IF THE CONTOUR LINE LENGTH IS TO HAVE A MINIMUM
C
      noclim=.false.
      if(clim.lt.1.e-30) noclim=.true.
C
C DETERMINE IF THE DISTANCE BETWEEN LABELS HAS BEEN SPECIFIED
C
      nodelb=.false.
      if(delb.lt.1.e-30) nodelb=.true.
C
C DETERMINE IF HACHERY IS TO BE USED AT ALL
C
      nohen=.false.
      if(nohach.and.noclim.and.nodelb.and.(icapt.eq.0)) nohen=.true.
      if(nohen) goto 199
C
C OPEN SCRATCH FILES
C
C --- BARN1 HOLDS TRACE END POINTS AND TRACE SEGMENT POINTERS ---
      open(unit=dsk1,form='unformatted',
     1 status='scratch',recl=24,
     2 access='direct')
C --- BARN2 HOLDS CONTOUR LINE POINTERS ---
      open(unit=dsk2,form='unformatted',
     1 status='scratch',recl=20,
     2 access='direct')
C --- PATH HOLDS TRACE LINKING POINTERS ---
      open(unit=dsk3,form='unformatted',
     1 status='scratch',recl=52,
     2 access='direct')
C --- THISTLES HOLDS NON-LOW CLOSURE SAMPLE POINTS ---
      open(unit=dsk5,form='unformatted',
     1 status='scratch',recl=8,
     2 access='direct')
C --- MEGARRAY HOLDS CONTOUR TRACE SEGMENTS ---
      call megaop(dsk4)
C
C INITIALIZE VARIABLES
C
C --- INITIALIZE HACHER DIMENSIONS ---
      glen=hlen
      gspace=hspace
C --- INITIALIZE CLOSURE RATIO ---
      gratio=hgap
C --- INITIALIZE MAX HACHURE IN A STACK ---
      mxhach=nint(hmax)
      if(mxhach.gt.999.or.mxhach.lt.0) mxhach=999
C --- INITIALIZE STACK QUALIFICATION MINIMUM CONTOUR LENGTH ---
      glim=hlim
C --- INITIALIZE MAXIMUM ENCLOSED HIGH LIMIT
      hachvb=hvb
C --- INITIALIZE MINIMUM CONTOUR LINE LENGTH ---
      conlim=clim
      if(glim.lt.conlim) glim=conlim
C --- INITIALIZE SPACE BETWEEN LABELS ---
      bspace=delb
      if(nodelb.or.delb.lt.1.e-30) bspace=7.0
C --- INITIALIZE FILE POINTERS ---
      jsrec=1
      nsrec=1
      ndrecs=0
      ndrecn=0
      ndrec=0
C --- INITIALIZE TIER AND STORAGE POINTERS ---
      south=dsk1
      north=dsk5
      itier=1
      ntier=npass
      nrowt1=nrowt-1
      ty1=ty0
      ty4=ty0
      dely=tdely
      do 10 i=1,nrowt1
        ty4=ty4+dely
   10 continue
C
C RETURN TO CALLING ROUTINE
C
  199 return
C
C________________________________________________________________
C
C  ENTRY  L A Y E G G
C________________________________________________________________
C
C ENTRY LAYEGG DETERMINES IF A COMPLETE TRACE MEETS THE INITIAL
C CONDITIONS FOR A LOW CLOSURE.  IF IT DOES NOT THE TRACE IS
C DISCARDED.  IF IT IS NOT A CLOSURE BUT IT ENTERS AND EXITS THE
C TIER AT THE SOUTH OR NORTH SIDE IT IS STORED FOR LATER
C EVALUATION.  IF IT IS A LOW CLOSURE IT IS HACHERED AND
C DISCARDED.
C
C  NTURN  - INTEGER*4.  THE SUM OF THE DIRECTIONS TAKEN BY THE
C           TRACE.  NEG=PREDOMINANT LEFT TURN.  ZERO=PREDOMINANT
C           STRAIGHT (NO TURN).  POS=PREDOMINANT RIGHT TURN.  A
C           CLOSURE WILL NECESSARILY BE NEG OR POS FOUR.  SINCE
C           THE CONTOUR IS DRAWN IN SUCH A WAY THAT THE SLOPE
C           ALWAYS DECREASES TO THE RIGHT A LOW CLOSURE WILL BE
C           SIGNIFIED BY NTURN=+4.
C  XB,XE  - REAL*4.  BEGINNING AND ENDING X-POSITIONS OF TRACE.
C  YB,YE  - REAL*4.                       Y-POSITIONS
C  TLEN   - REAL*4.  TOTAL LENGTH OF CONTOUR TRACE.
C  CONT   - REAL*4.  VALUE OF THE TRACE'S CONTOUR LEVEL.
C
      entry layegg(nturn,xb,xe,yb,ye,tlen,cont,linetp,ipost,size)
  300 continue
C
C DETERMINE IF HACHERS ARE TO BE MADE
C
      if(nohen) goto 399
C
C CHECK FOR COMPLETE CLOSURE OF THE TRACE WITHIN A SINGLE TIER.
C
      if(xe.eq.xb.and.ye.eq.yb) then
C ---   COMPLETE CLOSURES WILL BE PRELINKED ---
        lnkon=1
      else
C ---   NON-CLOSURES WILL BE LINKED AFTER MAP IS COMPLETED ---
        lnkon=-1
      end if
C
C STORE TRACE FOR LATER EVALUATION
C
      ndrec=ndrec+2
      ndrec1=ndrec-1
      ndrecb=ndrec1
      ndrece=ndrec
      if(lnkon.eq.1) then
C ---   STORE CLOSURES PRE-LINKED ---
        write(dsk3,rec=ndrecb) lnkon,ndrece,ndrece,
     1   nturn,nsrec,jsrec,tlen,cont,xb,yb,linetp,ipost,size
        write(dsk3,rec=ndrece) lnkon,ndrecb,ndrecb,
     1   nturn,nsrec,jsrec,tlen,cont,xe,ye,linetp,ipost,size
      else
C ---   STORE NON-CLOSURES UNLINKED ---
        write(dsk3,rec=ndrecb) lnkon,lnkon,ndrece,
     1   nturn,nsrec,jsrec,tlen,cont,xb,yb,linetp,ipost,size
        write(dsk3,rec=ndrece) lnkon,lnkon,ndrecb,
     1   nturn,nsrec,jsrec,tlen,cont,xe,ye,linetp,ipost,size
      end if
C
C DIFFERENTIATE BETWEEN NORTH AND SOUTH ENTRY AND EXIT POINTS
C
      if((ndrecs+1).ge.lental.or.(ndrecn+1).ge.lental) then
      write(*,*) 'tier boundary array size(lental value) exceeded'
      stop
      endif
      if(abs(yb-ty1).le.tol) then
C ---   TRACE BEGINNING AT SOUTH TIER BOUNDARY ---
        ndrecs=ndrecs+1
        tal(1,ndrecs)=xb
        ital(2,ndrecs)=ndrecb
      else if(abs(yb-ty4).le.tol) then
C ---   TRACE BEGINNING AT NORTH TIER BOUNDARY ---
        ndrecn=ndrecn+1
        write(north,rec=ndrecn) xb,ndrecb
      end if
      if(abs(ye-ty1).le.tol) then
C ---   TRACE ENDING AT SOUTH TIER BOUNDARY ---
        ndrecs=ndrecs+1
        tal(1,ndrecs)=xe
        ital(2,ndrecs)=ndrece
      else if(abs(ye-ty4).le.tol) then
C ---   TRACE ENDING AT NORTH TIER BOUNDARY ---
        ndrecn=ndrecn+1
        write(north,rec=ndrecn) xe,ndrece
      end if
C
C ADVANCE POINTER
C
  301 nsrec=jsrec
C
C RETURN TO CALLING ROUTINE
C
  399 return
C
C________________________________________________________________
C
C     ENTRY  G R A D E G G
C________________________________________________________________
C
C ENTRY GRADEGG LINKS ALL TRACE END POINTS AT THE SOUTH SIDE OF
C THE TIER WITH THEIR CONTINUATIONS AT THE NORTH SIDE OF THE 
C PREVIOUS TIER.  IT IS ASSUMED THAT ALL Y VALUES AT THE TIER
C BOUNDARY ARE EQUAL.  THEREFORE ONLY THE X VALUES NEED TO BE
C MATCHED IN ORDER TO FIND THE CONTINUATION.
C
C
      entry grdegg
  400 continue
C
C DETERMINE IF HACHERS ARE TO BE MADE
C
      if(nohen) goto 499
C
C IF NOT ENOUGH POINTS AT SOUTH BOUNDARY OF TIER TO LINK THEN
C PROCEDE TO NEXT TIER.
C
      if(ndrecs.le.1) goto 401
C
C SORT TRACE END POINTS AT SOUTH BOUNDARY OF TIER.
C
      nkey=1
      ikey(1)=1
      call sortil(tal,2,lental,ndrecs,ikey,nkey,lpoint,lenptr,
     1 ier)
C
C LINK ALL END POINTS AT SOUTH BOUNDARY WITH MATCHING END POINTS
C AT NORTH BOUNDARY OF PREVIOUS TIER.
C
      ll=lpoint(ndrecs+1)
      x1=tal(1,ll)
      ix1=ital(2,ll)
      do 11 i=2,ndrecs
        ll=lpoint(ll)
        x2=tal(1,ll)
        ix2=ital(2,ll)
        if(x2.eq.x1) then
          read(dsk3,rec=ix1) lnkon1,lnk1,il1,
     1     nturn1,nsrec1,jsrec1,tlen1,cont1,x1,y1,linet1,ipost1,
     2     size1
          read(dsk3,rec=ix2) lnkon2,lnk2,il2,
     1     nturn2,nsrec2,jsrec2,tlen2,cont2,x2,y2,linet2,ipost2,
     2     size2
          if(lnkon1.eq.-1.and.lnkon2.eq.-1) then
            lnkon1=1
            lnkon2=1
            write(dsk3,rec=ix1) lnkon1,ix2,il1,
     1       nturn1,nsrec1,jsrec1,tlen1,cont1,x1,y1,linet1,
     2       ipost1,size1
            write(dsk3,rec=ix2) lnkon2,ix1,il2,
     1       nturn2,nsrec2,jsrec2,tlen2,cont2,x2,y2,linet2,
     2       ipost2,size2
          end if
        end if
        x1=x2
        ix1=ix2
   11 continue
C
C INCREMENT TIER NUMBER AND SWITCH TRACE END POINT STORAGE FILES.
C
  401 if(itier.lt.ntier) then
        itier=itier+1
        ty1=ty4
        do 20 i=1,nrowt1
          ty4=ty4+dely
   20 continue
        do 30 i=1,ndrecn
          read(north,rec=i) tal(1,i),ital(2,i)
   30 continue
        ndrecs=ndrecn
        ndrecn=0
      end if
C
C RETURN TO CALLING ROUTINE
C
  499 return
C
C________________________________________________________________
C
C     ENTRY  G A T H E R E G G S
C________________________________________________________________
C
C ENTRY GATHEREGGS DETERMINES WHERE THE LOW CLOSURES ARE
C AFTER THE ENTIRE PLOT IS COMPLETE AND THEN RECORDS THEIR
C LOCATION ON DSK1 AND DSK2.
C
      entry gather
  500 continue
C
C CHECK IF HACHERS ARE TO BE MADE OR MINIMUM LENGTH TO BE CHECKED
C
      if(nohen) goto 699
C
C FOLLOW EACH OPEN LINKED LIST AND CHECK IF IT QUALIFIES AS A LOW
C CLOSURE BASED ON HACH GAP AND THE DIRECTION OF TURN.
C THEN FOLLOW ALL CLOSURES AND CHECK IF THEY TURN CLOCKWISE.
C SAVE QUALIFYING LOWS IN FILES WHICH REFERENCE THEM TO THE 
C MEGARRAY SEGMENT STORAGE FILE.
C
C --- INITIALIZE FILE POINTERS ---
      ist=0
      ist2=0
      ist3=0
      level=1
C
C CHECK FOR OPEN TRACES AND OPEN TRACE CHAINS
C
      do 31 i=1,ndrec1,2
C ---   READ IN BEGINNING OF A TRACE ---
        read(dsk3,rec=i) lnkon1,lnk1,il1,
     1   nturn1,nsrec1,jsrec1,tlen1,cont1,x1,y1,linet1,ipost1,
     2   size1
C ---   SKIP DELETED (0) OR LINKED (1) TRACES ---
        if(lnkon1.ne.-1) goto 513
        read(dsk3,rec=il1) lnkon2,lnk2,il2,
     1   nturn2,nsrec2,jsrec2,tlen2,cont2,x2,y2,linet2,ipost2,
     2   size2
        if(lnkon2.eq.1) goto 515
        if(lnkon2.eq.0) goto 513
C ---   DELETE BOTH ENDS OF TRACE ---
        write(dsk3,rec=i) jzero
        write(dsk3,rec=il1) jzero
C ---   CHECK GAP FOR TRACES WHICH ARE UNLINKED AT BOTH ENDS ---
        gap=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))
C ---   STORE POINTERS FOR CONTOUR LINE ---
        ist=ist+1
        write(dsk1,rec=ist) nsrec1,jsrec1,cont1,linet1,ipost1,
     1   size1
        if(nturn1.gt.0.and.gap.le.gratio*tlen1) then
C ---     STORE LOW CLOSURE ---
          ist2=ist2+1
          write(dsk2,rec=ist2) cont1,tlen1,ist,ist,level
        else if(tlen1.gt.conlim.and.tlen1.gt.hachvb) then
C ---     STORE SAMPLE POINT FOR NON-LOW CLOSURE ---
          ist3=ist3+1
          write(dsk5,rec=ist3) x1,y1
        end if
        goto 513
C ---   FOLLOW CHAIN OF TRACES UNTIL UNLINKED END ---
  515   jst=0
        nturnt=0
        tlent=0.
        x0=x1
        y0=y1
        j1=i
  512   jst=jst+1
C ---   KEEP RECORD OF EACH LINKED TRACE IN CHAIN ---
        write(dsk1,rec=ist+jst) nsrec1,jsrec1,cont1,linet1,
     1   ipost1,size1
        jrev=il1-j1
        nturnt=nturnt+nturn1*jrev
        tlent=tlent+tlen1
        j2=il1
C ---   READ IN END OF TRACE ---
        read(dsk3,rec=j2) lnkon2,lnk2,il2,
     1   nturn2,nsrec2,jsrec2,tlen2,cont2,x2,y2,linet2,ipost2,
     2   size2
C ---   DELETE BOTH ENDS OF TRACE ---
        write(dsk3,rec=j1) jzero
        write(dsk3,rec=j2) jzero
C ---   IF UNLINKED END FIND GAP AND STOP FOLLOWING CHAIN ---
        if(lnkon2.ne.1) then
          gap=sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))
          goto 514
        end if
C ---   IF LINKED BACK ON ITSELF (SHOULD NEVER HAPPEN) STOP ---
        if(lnk2.eq.i) goto 514
        j1=lnk2
C ---   READ IN BEGINNING OF NEXT TRACE ---
        read(dsk3,rec=j1) lnkon1,lnk1,il1,
     1   nturn1,nsrec1,jsrec1,tlen1,cont1,x1,y1,linet1,ipost1,
     2   size1
C ---   IF UNLINKED BEGINNING FIND GAP AND STOP FOLLOWING ---
        if(lnkon1.ne.1) then
          gap=sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))
          goto 514
        end if
        goto 512
  514   if(nturnt.gt.0.and.gap.le.gratio*tlent) then
C ---     STORE POINTERS FOR LOW CLOSURE CHAIN ---
          ist2=ist2+1
          write(dsk2,rec=ist2) cont1,tlent,ist+1,ist+jst,level
        else if(tlent.gt.conlim.and.tlent.gt.hachvb) then
C ---     STORE SAMPLE POINT FOR NON-LOW CLOSURE CHAIN ---
          ist3=ist3+1
          write(dsk5,rec=ist3) x0,y0
        end if
C ---   STORE CONTOUR CHAIN ---
        ist=ist+jst
  513   continue
   31 continue
C
C CHECK FOR CLOSED TRACES AND CLOSED TRACE CHAINS
C
  501 do 32 i=1,ndrec1,2
C ---   READ IN BEGINNING OF TRACE ---
        read(dsk3,rec=i) lnkon1,lnk1,il1,
     1   nturn1,nsrec1,jsrec1,tlen1,cont1,x1,y1,linet1,ipost1,
     2   size1
C ---   BREAK INTO LINKED SECTION ---
        if(lnkon1.eq.1) then
          jst=0
          nturnt=0
          tlent=0.
          x0=x1
          y0=y1
          j1=i
  502     jst=jst+1
C ---     KEEP RECORD OF EACH LINKED TRACE ---
          write(dsk1,rec=ist+jst) nsrec1,jsrec1,cont1,linet1,
     1     ipost1,size1
          jrev=il1-j1
          nturnt=nturnt+nturn1*jrev
          tlent=tlent+tlen1
          j2=il1
C ---     READ IN END OF TRACE ---
          read(dsk3,rec=j2) lnkon2,lnk2,il2,
     1     nturn2,nsrec2,jsrec2,tlen2,cont2,x2,y2,linet2,ipost2,
     2     size2
C ---     DELETE TRACE ---
          write(dsk3,rec=j1) jzero
          write(dsk3,rec=j2) jzero
C ---     IF UNLINKED END STOP FOLLOWING CHAIN ---
          if(lnkon2.ne.1) goto 503
C ---     CHAIN CLOSES ON ITSELF ---
          if(lnk2.eq.i) goto 504
          j1=lnk2
C ---     READ BEGINNING OF NEXT TRACE ---
          read(dsk3,rec=j1) lnkon1,lnk1,il1,
     1     nturn1,nsrec1,jsrec1,tlen1,cont1,x1,y1,linet1,ipost1,
     2     size1
C ---     IF BEGINNING IS UNLINKED STOP FOLLOWING CHAIN ---
          if(lnkon1.ne.1) goto 503
          goto 502
  504     if(ist2.ge.lental.or.ist3.ge.lental) then
          write(*,*) 'closure array size(lenclo value) exceeded'
          stop
          endif
          if(nturnt.gt.0) then
C ---       STORE POINTERS FOR LOW CLOSURES ---
            ist2=ist2+1
            write(dsk2,rec=ist2) cont1,tlent,ist+1,ist+jst,level
          else if(tlent.gt.conlim.and.tlent.gt.hachvb) then
C ---       STORE SAMPLE POINT FOR NON-LOW CLOSURES ---
            ist3=ist3+1
C           if(gratio.ge.1.e-30) write(dsk5,rec=ist3) x0,y0
            write(dsk5,rec=ist3) x0,y0
          end if
C ---     STORE CONTOUR CHAIN ---
          ist=ist+jst
        end if
  503   continue
   32 continue
      close(unit=dsk3)
C
C________________________________________________________________
C
C     ENTRY  E G G Y O L K
C________________________________________________________________
C
C ENTRY EGGYOLK SORTS EACH QUALIFYING LOW CLOSURE BY TOTAL LENGTH
C (SECONDARY SORT) AND CONTOUR LEVEL (PRIMARY SORT) AND THEN
C STARTING AT THE LOWEST CONTOUR LEVEL AND SHORTEST LENGTH
C CONSTRUCTS A HEIRARCHY OF CONTOUR LEVELS.  THE LEVEL OF EACH
C CONTOUR CLOSURE IS WRITTEN BACK INTO THE FILE WHICH KEEPS TALLY
C OF THEM.
C
C
      entry eggylk
  700 continue
      if(nohach) go to 600
C
C READ IN NON-LOW CLOSURE LIST AND SORT ON X AND Y
C
C      if(gratio.lt.1.e-30) goto 701
      do 33 i=1,ist3
        read(dsk5,rec=i) tal(1,i),tal(2,i)
   33 continue
      nkey=2
      ikey(1)=1
      ikey(2)=2
      call sortil(tal,2,lental,ist3,ikey,nkey,lpoint,lenptr,
     1 ier)
      ll=lpoint(ist3+1)
      do 34 i=1,ist3
        tal2(1,i)=tal(1,ll)
        tal2(2,i)=tal(2,ll)
        ll=lpoint(ll)
   34 continue
  701 continue
C
C READ IN LOW CLOSURE LIST AND SORT ON CONTOUR VALUE AND LENGTH
C
      do 35 i=1,ist2
        read(dsk2,rec=i) tal(1,i),tal(2,i),idum3,idum4,idum5
   35 continue
      if(mxhach.lt.999) then
        nkey=2
        ikey(1)=1
        ikey(2)=2
        call sortil(tal,2,lental,ist2,ikey,nkey,lpoint,lenptr,
     1   ier)
      else
        lpoint(ist2+1)=1
        do 36 i=1,ist2
          lpoint(i)=i+1
   36   continue
      end if
C
C LOAD ALL SEGMENTS OF EACH CLOSURE INTO AN ARRAY
C
C --- INITIALIZE NUMBER OF INDEPENDANT CLOSURE STACKS ---
      nstack=0
C --- FIND POINTER OF THE LOWEST AND SHORTEST CLOSURE ---
      ll=lpoint(ist2+1)
C --- LOOP THROUGH CHECKING ONE COMPLETE CLOSURE AT A TIME ---
      do 105 jll=1,ist2
C ---   READ IN POINTERS LOCATIONS FOR THE NEXT LOW CLOSURE ---
        read(dsk2,rec=ll) conval,conlen,ibeg,iend
C ---   INITIALIZE OFFSET POINTER IN ARRAYS FOR HOLDING THE LOW
C       CLOSURE POINTS ---
        joff=0
C ---   LOOP THROUGH LOW CLOSURE POINTERS ---
        do 60 i=ibeg,iend
          read(dsk1,rec=i) nsrec,jsrec
C ---     INITIALIZE MEGARRAY POINTER ---
          iptr=nsrec
          skip=.false.
   45     if(iptr.ge.jsrec) go to 60
            call megard(iptr,nparm,iparm)
            if(skip) then
C ---         SKIP SEGMENT BECAUSE IT'S ALREADY BEEN LOADED ---
              iskip=iskip+1
              if(iskip.gt.nskip) skip=.false.
            else
C ---         LOAD SEGMENT INTO ARRAY ---
              iptr=iptr+nparm
              npts=iparm(2)
              call rmegar(iptr,npts,xx)
              call rmegar(iptr+npts,npts,yy)
              do 50 j=1,npts
                xa(joff+j)=xx(j)
                ya(joff+j)=yy(j)
   50         continue
C ---         INCREMENT OFFSET POINTER ---
              joff=joff+npts
C ---         CHECK IF SEGMENT CONTAINS SOME OTHER SEGMENTS ---
              isegtp=iparm(4)
              if(isegtp.gt.100) then
                skip=.true.
                iskip=1
                nskip=isegtp-100
              end if
            end if
C ---       POSITION MEGARRAY POINTER TO BEG OF NEXT SEGMENT ---
            iptr=iparm(1)
         go to 45
   60   continue
C
C CHECK FOR GAPS AT CORNERS AND MAKE THEM SQUARE
C
        if(xa(joff)-xa(1)) 711,715,712
  711   if(ya(joff)-ya(1)) 714,715,713
  712   if(ya(joff)-ya(1)) 713,715,714
C ---   ADD POINT TO MAKE SQUARE GAPS AT SW AND NE CORNERS ---
  713   joff=joff+1
        xa(joff)=xa(1)
        ya(joff)=ya(joff-1)
        goto 715
C ---   ADD POINT TO MAKE SQUARE GAPS AT NW AND SE CORNERS ---
  714   joff=joff+1
        xa(joff)=xa(joff-1)
        ya(joff)=ya(1)
  715   continue
C
C FIND MIN AND MAX X AND Y IN CURRENT LOW CLOSURE
C
        xamin=xa(1)
        xamax=xamin
        yamin=ya(1)
        yamax=yamin
        do 70 i=2,joff
          if(xa(i).lt.xamin) then
            xamin=xa(i)
          else if(xa(i).gt.xamax) then
            xamax=xa(i)
          end if
          if(ya(i).lt.yamin) then
            yamin=ya(i)
          else if(ya(i).gt.yamax) then
            yamax=ya(i)
          end if
   70 continue
c
C FIND IF A NON-LOW CLOSURE IS CONTAINED IN THE LOW CLOSURE
C
C        if(gratio.lt.1.e-30) goto 721
        call binsrc(tal2,ist3,1,xamin,xamax,minx,maxx)
        do 80 i=minx+1,maxx-1
          if(tal2(2,i).ge.yamin.and.tal2(2,i).le.yamax) then
      if(inside(xa,ya,joff,tal2(1,i),tal2(2,i)).ne.0) then
              level=1000
              goto 731
            end if
          end if
   80   continue
  721     continue
C
C FIND LEVEL IN STACK OF PRESENT LOW CLOSURE
C
      if(mxhach.ge.999) goto 732
C ---   INITIALIZE STACK POINTER ---
        jstack=1
C ---   INITIALIZE STACK POSITION HOLDING POINTER ---
        kstack=0
   85   if(jstack.gt.nstack) go to 95
          if(tal(1,jstack).lt.xamin) goto 741
          if(tal(1,jstack).gt.xamax) goto 741
          if(tal(2,jstack).lt.yamin) goto 741
          if(tal(2,jstack).gt.yamax) goto 741
      if(inside(xa,ya,joff,tal(1,jstack),tal(2,jstack)).ne.0) then
C ---       PRESENT STACK IS CONTAINED IN PRESENT LOW CLOSURE ---
            if(kstack.eq.0) then
C ---         SET STACK POSITION HOLDING POINTER ---
              kstack=jstack
            else
C ---         MERGE STACKS BY SETTING MAXIMUM STACK LEVEL IN
C             FIRST STACK POSITION ---
              if(istack(jstack).gt.istack(kstack))
     1         istack(kstack)=istack(jstack)
C ---         ELIMINATE SECOND STACK AFTER MERGE ---
              nstack=nstack-1
              do 90 j=jstack,nstack
                j1=j+1
                istack(j)=istack(j1)
                tal(1,j)=tal(1,j1)
                tal(2,j)=tal(2,j1)
   90         continue
              jstack=jstack-1
            end if
          end if
C ---     INCREMENT STACK POINTER ---
  741     jstack=jstack+1
        go to 85
   95   continue
        if(kstack.eq.0) then
C ---     BEGIN A NEW CLOSURE STACK ---
          nstack=nstack+1
          kstack=nstack
          istack(kstack)=1
          tal(1,kstack)=xa(2)
          tal(2,kstack)=ya(2)
        else
C ---     ADD NEW LEVEL TO EXISTING CLOSURE STACK ---
          istack(kstack)=istack(kstack)+1
        end if
C ---   IF CONTOUR DOES NOT QUALIFY THEN DO NOT COUNT IT ---
        if(conlen.lt.glim) istack(kstack)=istack(kstack)-1
C ---   WRITE OUT CLOSURE LEVEL AND POINTER LOCATIONS ---
        level=istack(kstack)
  731   write(dsk2,rec=ll) conval,conlen,ibeg,iend,level
C ---   FIND POINTER OF NEXT LOWEST AND SHORTEST CLOSURE ---
  732   ll=lpoint(ll)
  105 continue
C
C________________________________________________________________
C
C     ENTRY  H C O N P L
C________________________________________________________________
C
C ENTRY HCONPL CALLS AND CONTROLS BCONPL, THE SUBROUTINE
C WHICH DRAWS CONTOUR LINES, LABELS, AND HACHER TICKS.
C
      entry hconpl
  600 close(unit=dsk5)
C
C DETERMINE IF HACHERS ARE TO BE MADE
C
      if(nohach) goto 601
C
C FROM THE FILE WHICH CONTAINS ONE ENTIRE CLOSURE PER RECORD
C READ THE BEGINNING AND ENDING POINTERS FOR THE FILE WHICH
C CONTAINS ONE TRACE PER RECORD
C
      do 38 k=1,ist2
        read(dsk2,rec=k) dum1,dum2,ibeg,iend,level
C ---   DETERMINE IF THE LEVEL IS WITHIN ACCEPTABLE LIMITS ---
        if(level.le.mxhach) then
C
C FROM THE FILE WHICH CONTAINS ONE TRACE PER RECORD READ THE
C BEGINNING AND ENDING POINTERS FOR THE ARRAY WHICH CONTAINS THE
C TRACE SEGMENTS
C
          do 37 j=ibeg,iend
            read(dsk1,rec=j) nsrec,jsrec,contlb,linet1,ipost1,
     1       size1
            write(dsk1,rec=j) -1*nsrec,jsrec,contlb,linet1,
     1       ipost1,size1
   37     continue
        end if
   38 continue
  601 close(unit=dsk2)
      call bconpl(dsk1,ist,nparm,iparm,xa,ya,conlim,
     1 bspace,glen,gspace)
C
C RETURN TO CALLING ROUTINE
C
      close(unit=dsk1)
      call megacl(dsk4)
  699 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  H A C H
C________________________________________________________________
C
C SUBROUTINE HACH DRAWS HACHER MARKS ON LINES WHEN CALLED BY A
C PROGRAM USING THE RON WAHL PLOT SYSTEM.
C FOLLOWING IS A DESCRIPTION OF THE PASSED VARIABLES.
C
C   X,Y    - REAL*4. ARRAYS CONTAINING X AND Y LOCATIONS FOR THE
C            POINTS IN A SEGMENT OF THE CONTOUR TO BE HACHURED.
C            THE LOW SIDE IS ASSUMED TO BE ON THE RIGHT (HACHURE
C            TICS WILL BE PLACED ON THE RIGHT SIDE OF THE LINE).
C            IF THE SEGMENT IS SUPPOSED TO BE CLOSED THE FIRST
C            POINT MUST BE THE SAME AS THE LAST.
C   NPTS   - INTEGER*4.  THE NUMBER OF POINTS IN THE SEGMENT.
C   SDIST  - REAL*4.  THE CURRENT LENGTH OF THE CONTOUR.
C            (I.E. PREVIOUS SEGMENTS MAY HAVE ALREADY BEEN DRAWN)
C   HDIST  - REAL*4.  THE CURRENT LENGTH ALONG THE CONTOUR FROM
C            THE BEGINNING TO THE FIRST HACHURE TIC ON THIS
C            SEGMENT.  HDIST MUST BE GREATER THAN OR EQUAL TO
C            SDIST.
C   HLEN   - REAL*4.  LENGTH (IN MAP INCHES) OF THE HACHERS TICS.
C   HSPACE - REAL*4.  DISTANCE (IN MAP INCHES) BETWEEN HACHER
C            TICS.
C   LINETP - INTEGER*4.  A NUMBER DESIGNATING THE LINE TYPE OF
C            THE HACHURE TICS.
C
C SUBROUTINE HACH WRITTEN BY ROB BRACKEN, USGS, 24MAR87
C
C
      subroutine hach(x,y,npts,sdist,hdist,hlen,hspace,linetp)
      save
C
C VARIABLE DECLARATION
C
      real*4 x(1),y(1)
      real*4 xh(2),yh(2)
C
C CHECK, ADJUST, AND INITIALIZE VARIABLES
C
C --- CHECK HACHURE TIC SPACING ---
      if(hspace.lt.1.e-30) goto 199
C --- ADVANCE HACHURE TIC DISTANCE TO WITHIN ONE TIC SPACING
C     BEYOND CURRENT DISTANCE ---
  105 if(sdist.gt.hdist) then
        hdist=hdist+hspace
        goto 105
      end if
C --- INITIALIZE INDEX OF CURRENT POINT ALONG SEGMENT ---
      j=1
C
C FOLLOW SEGMENT TO ONE POINT BEYOND NEXT HACHER TIC
C
  106 if(sdist.gt.hdist) go to 107
        j=j+1
C ---   EXIT IF NO MORE POINTS IN SEGMENT ---
  104   if(j.gt.npts) goto 199
C ---   FIND DISTANCE BETWEEN CURRENT AND PREVIOUS POINT ---
        dx=x(j-1)-x(j)
        dy=y(j-1)-y(j)
        ds=sqrt(dx*dx+dy*dy)
        if(ds.lt.1.e-30) then
C ---     ELIMINATE DUPLICATE POINT ---
          npts=npts-1
          do 10 i=j,npts
            x(i)=x(i+1)
            y(i)=y(i+1)
   10     continue
          goto 104
        end if
C ---   INCREMENT CURRENT CONTOUR DISTANCE ---
        sdist=sdist+ds
      go to 106
  107 continue
C
C PUT ONE HACHURE TIC ON SEGMENT
C
C --- FIND POSITION ON CONTOUR LINE ---
      t=(sdist-hdist)/ds
      xh(1)=dx*t+x(j)
      yh(1)=dy*t+y(j)
C --- FIND DIRECTION OF PERPENDICULAR TO CONTOUR ---
      cosi=dx/ds
      sine=dy/ds
      cosip=-1*sine
      sinep=cosi
C --- FIND POSITION HLEN AWAY FROM CONTOUR.  ON RIGHT SIDE
C     LOOKING ALONG CONTOUR LINE (OR ON LEFT LOOKING BACK) ---
      xh(2)=cosip*hlen+xh(1)
      yh(2)=sinep*hlen+yh(1)
C
C................................................................
C
C     CALL TO PLOT SYSTEM -- LINE PLOTTING ROUTINE
C
      call line(xh,yh,2,0,linetp)
C
C ...............................................................
C
C --- INCREMENT HACHURE TIC DISTANCE TO NEXT TIC POSITION ---
      hdist=hdist+hspace
      goto 106
C
C RETURN TO CALLING ROUTINE
C
  199 return
      end
C
C________________________________________________________________
C
C   SUBROUTINE  F E E D
C________________________________________________________________
C
C SUBROUTINE FEED STORES CONTOUR LINE SEGMENTS IN A SEQUENTIAL
C FILE UNTIL THEY ARE DISCARDED OR USED FOR HACHER POSTITIONS.
C
C   X      - REAL*4 ARRAY.  ARRAY OF CONTOUR LINE SEGMENT
C                           X POSITIONS.
C   Y      - REAL*4 ARRAY.  Y POSITIONS.
C   NPTS   - INTEGER*4.  NUMBER OF POINTS IN SEGMENT.
C   IC     - INTEGER*4.  0=BEGIN NEW SEGMENT.  1=CONTINUE FROM
C            PREVIOUS SEGMENT.
C   LINETP - INTEGER*4.  TYPE OF LINE:  PLOT SYSTEM DEPENDANT
C            -999 - SEGMENT IS NOT TO BE PLOTTED
C   ISEGTP - INTEGER*4.  TYPE OF SEGMENT:  1=MID LINE,
C            2=BEFORE GRADIENT FLAG,  3=END OF TRACE,
C            11-13 = BEFORE A LABEL,  21-23 = AFTER A LABEL.
C            31-33 = IMMEDIATELY AFTER A LABEL.
C            41-43 = POINTS REMOVED FOR THE LABEL.
C            100+KK = SEGMENT INCLUDES KK FOLLOWING SEGMENTS
C
C
      subroutine feed(x,y,npts,ic,linetp,isegtp)
      save
C
C VARIABLE DECLARATION
C
      integer*4 dsk4
      parameter(dsk4=17,nparm=5)
      real*4 x(1),y(1)
      integer*4 iparm(nparm)
      common /eggchu/ jemeg,nohach,noclim,nodelb,nohen
      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
      logical*4 nohach,noclim,nodelb,nohen
C
C DETERMINE IF HACHERS ARE TO BE MADE OR MINIMUM CONTOUR LENGTH
C HAS BEEN SPECIFIED
C
      if(nohen) then
C
C................................................................
C
C       CALL TO PLOT SYSTEM -- LINE PLOTTING ROUTINE
C
        if(iplotr.eq.5) then
        ipen=linetp/100
        if(ipen.lt.1) ipen=1
        if(ipen.gt.npens) ipen=npens
        call newpen(ipen)
        endif
        if(linetp.gt.-999) call line(x,y,npts,ic,linetp)
C
C................................................................
C
        goto 299
      end if
C
C WRITE TO FILE AND INCREMENT COUNTER
C
      j=jemeg
      jemeg=jemeg+nparm+npts+npts
      iparm(1)=jemeg
      iparm(2)=npts
      iparm(3)=ic
      iparm(4)=isegtp
      iparm(5)=linetp
      call megawr(j,nparm,iparm)
      j=j+nparm
      call rmegaw(j,npts,x)
      j=j+npts
      call rmegaw(j,npts,y)
  299 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE   S O R T I L
C________________________________________________________________
C
c  Hart modification to Wooodrum sort from Creative Computing
c  1978,vol 4,no. 1,p. 96-101.
c  Translated to FORTRAN by R.R. Wahl UTR Branch
c  US Geol. Survey, Denver, CO
c
c     modified by r.godson,usgs
c     to handle multi-dimensioned input arrays
c     and return sorted data in another array
c     equivalent in size to the input array    
c
c     return array option is deleted in this version
c
C
C  DESCRIPTION OF PASSED VARIABLES:
C
C     AIN    - REAL*4.  2 DIMENSIONAL INPUT UNSORTED ARRAY.
C     IDIM   - INTEGER*4.  FIRST DIMENSION OF AIN.
C     JDIM   - INTEGER*4.  SECOND DIMENSION OF AIN.
C     NSORT  - INTEGER*4.  NUMBER OF COLUMNS TO SORT.
C     IKEY   - INTEGER*4.  1 DIMENSIONAL ARRAY CONTAINING A LIST
C              OF KEY FIELDS (ROW NUMBERS, FIRST INDEX) TO SORT
C              FROM MAJOR TO MINOR.
C     NKEY   - INTEGER*4.  THE NUMBER OF KEY FIELDS.
C     LNK    - INTEGER*4.  1 DIMENSIONAL OUTPUT LINK ARRAY.
C              AIN(ISORT,LNK(NSORT+1)) = LEAST VALUE.
C              AIN(ISORT,LNK(LNK(NSORT+1))) = SECOND LEAST VALUE.
C     LDIM   - INTEGER*4.  DIMENSION OF LINK ARRAY.
C              MINIMUM LDIM = NSORT + LN( NSORT-2 ) + 2
C     IER    - INTEGER*4.  ERROR INDICATOR.  0 = SUCCESSFUL RUN.
C              1 = LESS THAN 2 ELEMENTS TO SORT.  2 = NOT ENOUGH
C              WORK SPACE IN LINK ARRAY.
C
C
      subroutine sortil(ain,idim,jdim,
     1                  nsort,ikey,nkey,lnk,ldim,ier)
C
      real ain
      dimension ain(idim,jdim),lnk(ldim),ikey(nkey)
c
      ier=0
      if (nsort.le.1)  go to 1000
      mt=nsort+ifix(alog(float(nsort))/.69314718+.5)+2
      if (ldim.lt.mt)  go to 1010
c
c  Initialize variables
      k1=0
      i=0
      m1=0
      it2=0
      it4=0
c     Head of sequence 1
      j=nsort+1
      lnk(1)=1
      lnk(j)=1
      k2=1
c     Number of leaves
      nleav=nsort
c
c  Climb the tree
  250 continue
c     Low-order twig value
      if (nleav.lt.4)  go to 320
      ib2=nleav/2
      it4=it4+k2*(nleav-2*ib2)
      nleav=ib2
c     Total number of twigs at an altitude
      k2=k2*2
      go to 250
c
c  Initial calculations
  320 continue
c     Number of low-order twigs
      it4=k2-it4
c     High bit value of binary counter
      ib2=k2/2
c
c  Next twig
  350 continue
c     Exit -- sort complete
      if (k1.eq.k2)  go to 940
c     Twig number
      k1=k1+1
      it1=k1
c     High bit value
      ib1=ib2
c     Previous reflected twig number
      it3=it2
c
c  Add 1 to reflected binary counter and carry
  400 continue
c     Any more carries
      it1s=it1/2
      if (2*it1s.lt.it1)  go to 470
      it1=it1s
c     Number of merges
      m1=m1+1
      it2=it2-ib1
c     Next bit value
      ib1=ib1/2
c     Carry one
      go to 400
c
c  Twig calculations
  470 continue
c     Reflected twig number
      it2=it2+ib1
c     2-twigs and 3-twigs
      if (nleav.eq.2)  go to 550
c
c  3-twigs and 4-twigs
c     low-order twig (3-twig)
      if (it3.lt.it4)  go to 560
c
c  4-twig
c     Dis-engage number of merges
      m1=-m1
      go to 630
  550 continue
c     Low-order twig (2-twig)
      if (it3.lt.it4)  go to 610
c
c  3-twig
  560 continue
c     Number of merges
      m1=m1+1
c     Next leaf
      i=i+1
c     Generate a leaf
      lnk(i)=i
      lnk(j)=i
c     Generate next sequence head
      j=j+1
c
c  2-twig
  610 continue
c     Number of merges
      m1=m1+1
  630 continue
c     Next leaf
      i=i+1
c     Generate a new leaf
      l1=i
      lnk(i)=i
      lnk(j)=i
c     Head of older leaf (last line)
      l0=j
c     Head of latest leaf (next two lines)
      j=j+1
c     Next leaf
      i=i+1
c     Generate a leaf
      l2=i
      lnk(i)=i
      lnk(j)=i
c     Merge leaves
      go to 750
c
c  Merge twigs and branches
  700 continue
c     Head of latest branch or twig
      j=j-1
c     Head of older branch or twig
      l0=j-1
c     Head of sequence 1
      l1=lnk(l0)
c     Head of sequence 2
      l2=lnk(j)
  750 continue
c     Stay in sequence 1
C
      kkey=1
  763 if(ain(ikey(kkey),l1)-ain(ikey(kkey),l2)) 820,762,761
  762 kkey=kkey+1
      if(kkey.le.nkey) goto 763
      goto 820
  761 continue
C
C      if (ain(isort,l1).le.ain(isort,l2))  go to 820
c     Switch to sequence 2
      lnk(l0)=l2
  770 continue
c     Top leaf in sequence 2
      l0=l2
c     Next leaf in sequence 2
      l2=lnk(l0)
c     End of sequence 2
      if (l2.eq.l0)  go to 870
c     Stay in sequence 2
C
      kkey=1
  783 if(ain(ikey(kkey),l1)-ain(ikey(kkey),l2)) 781,782,770
  782 kkey=kkey+1
      if(kkey.le.nkey) goto 783
  781 continue
C
C      if (ain(isort,l1).gt.ain(isort,l2))  go to 770
c     Switch to sequence 1
      lnk(l0)=l1
  820 continue
c     Top leaf in sequence 1
      l0=l1
c     Next leaf in sequence 2
      l1=lnk(l0)
c     Not end of sequence 1
      if (l1.ne.l0)  go to 750
c     Switch to sequence 2
      lnk(l0)=l2
      go to 880
  870 continue
c     Switch to sequence 1
      lnk(l0)=l1
  880 continue
c     Number of merges
      m1=m1-1
      if (m1.gt.0)  go to 700
      if (m1.eq.0)  go to 350
c
c  Generate second half of 4-twig
c     Re-engage number of merges
      m1=1-m1
      go to 630
c
c  All done
c
c     sort output array
c
  940 ll=lnk(nsort+1)
c      do 960 i=1,nsort
c      do 950 j=1,jdim
c      aout(i,j)=ain(ll,j)
c  950 continue
c      ll=lnk(ll)
c  960 continue
  970 return
c     Less than two elements to sort
 1000 ier=1
      lnk(nsort+1)=nsort
      lnk(nsort)=nsort+1
      go to 970
c     Not enough work space in link array
 1010 ier=2
      go to 970
      end
C
C________________________________________________________________
C
C     SUBROUTINE  B I N S R C
C________________________________________________________________
C
C SUBROUTINE BINARY SEARCH USES A BINARY WALK TO SEARCH FOR
C MINIMUM AND MAXIMUM LIMITS IN A PRESORTED ARRAY.
C
C    ARRAY  - REAL*4.  A PRE-SORTED ARRAY CONTAINING VALUES OF
C             WHICH TO FIND THE POSITIONS.
C    NPTS   - INTEGER*4.  NUMBER OF POINTS IN ARRAY.
C    KEY    - INTEGER*4.  THE KEY FIELD TO SEARCH.
C    AMIN   - REAL*4.  THE VALUE OF THE LOWER LIMIT.
C    AMAX   - REAL*4.  THE VALUE OF THE UPPER LIMIT.
C    IMIN   - INTEGER*4.  THE POSITION IN ARRAY IMMEDIATELY
C             BEFORE AMIN.  IF IMIN IS GREATER THAN NPTS THEN ALL
C             VALUES IN ARRAY ARE LESS THAN AMIN.  IF IMIN IS
C             LESS THAN 1 THEN ALL VALUES IN ARRAY ARE GREATER
C             THAN AMIN.
C    IMAX   - INTEGER*4.  THE POSITION IN ARRAY IMMEDIATELY
C             AFTER AMAX.  IF IMAX IS GREATER THAN NPTS THEN ALL
C             VALUES IN ARRAY ARE LESS THAN AMAX.  IF IMAX IS
C             LESS THAN 1 THEN ALL VALUES IN ARRAY ARE GREATER
C             THAN AMAX.
C
C SUBROUTINE BINSRC WRITTEN BY ROB BRACKEN, USGS, 5MAY87.
C
C
      subroutine binsrc(array,npts,key,amin,amax,imin,imax)
      save
C
      real*4 array(2,*)
C
C FIND IMIN
C
C --- FIND POWER OF 2 LIMIT FOR IMIN ---
      limit=1
  104 if(limit.gt.npts) goto 101
      if(array(key,limit).ge.amin) goto 101
      limit=limit*2
      goto 104
C --- SET UP BEGINNING IMIN IN STRATIGIC LOCATION ---
  101 ndiff=limit/4
      imin=limit-ndiff
C --- SET DISTANCE TO WALK ON THE NEXT STEP ---
  102 ndiff=ndiff/2
      if(ndiff.lt.1) goto 103
C --- DETERMINE DIRECTION OF WALK AND MAKE A STEP ---
      if(imin.gt.npts) then
        imin=imin-ndiff
      else if(array(key,imin).ge.amin) then
        imin=imin-ndiff
      else
        imin=imin+ndiff
      end if
      goto 102
C --- CLEAN UP EDGE EFFECTS ---
  103 if(imin.gt.npts) then
        imin=npts
      else if(imin.lt.1) then
        imin=1
      end if
      if(array(key,imin).ge.amin) imin=imin-1
      if(imin.eq.npts) then
        imax=npts+1
        goto 199
      end if
C
C FIND IMAX
C
C --- FIND POWER OF 2 LIMIT FOR IMAX ---
      if(limit.ge.2) limit=limit/2
  114 if(limit.gt.npts) goto 111
      if(array(key,limit).gt.amax) goto 111
      limit=limit*2
      goto 114
C --- SET UP BEGINNING IMAX IN STRATIGIC LOCATION ---
  111 ndiff=limit/4
      imax=limit-ndiff
C --- SET DISTANCE TO WALK ON THE NEXT STEP ---
  112 ndiff=ndiff/2
      if(ndiff.lt.1) goto 113
C --- DETERMINE DIRECTION OF WALK AND MAKE A STEP ---
      if(imax.gt.npts) then
        imax=imax-ndiff
      else if(array(key,imax).gt.amax) then
        imax=imax-ndiff
      else
        imax=imax+ndiff
      end if
      goto 112
C --- CLEAN UP EDGE EFFECTS ---
  113 if(imax.gt.npts) then
        imax=npts
      else if(imax.lt.1) then
        imax=1
      end if
      if(array(key,imax).le.amax) imax=imax+1
      if(imax.eq.1) imin=0
C
C EXIT PROCEDURE
C
  199 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  B C O N P L 
C________________________________________________________________
C
C SUBROUTINE B CONTOUR PLOT IS DESIGNED TO FIND WHETHER OR NOT
C THE CONTOUR LINES ARE LONGER THAN THE MINIMUM CONTOUR LENGTH
C (CONLIM) AND PLOT THEM.  SHORTER CONTOUR LINES WILL NOT BE
C PLOTTED.  LABELS AND HACHURE TICS WILL BE ADDED TO THE CONTOUR
C LINES WHEN APPROPRIATE.
C
C   IUNIT  - INTEGER*4.  THE UNIT NUMBER OF THE DEVICE TO READ
C            ARRAY POINTERS FROM.
C   IST    - INTEGER*4.  THE TOTAL NUMBER OF POINTERS TO READ
C            FROM IUNIT.
C   NPARM  - INTEGER*4.  THE NUMBER OF PARAMETERS TO READ BEFORE
C            EACH LINE SEGMENT.
C   IPARM  - INTEGER*4.  THE ARRAY CONTAINING THE PARAMETERS.
C   X,Y    - REAL*4.  ARRAYS CONTAINING THE X AND Y COORDINATES
C            FOR THE POINTS IN THE SEGMENTS.
C   CONLIM - REAL*4.  THE MINIMUM ALLOWABLE LENGTH FOR A
C            CONTINUOUS CONTOUR SEGMENT.
C   BSPACE - REAL*4.  THE DESIRED DISTANCE IN MAP INCHES BETWEEN
C            CONTOUR LABELS.
C   GLEN   - REAL*4.  THE LENGTH OF HACHURE TICS.
C   GSPACE - REAL*4.  THE DISTANCE BETWEEN HACHURE TICS.  IF
C            GSPACE=0. HACHURE TICS WILL NOT BE MADE AT ALL.
C
C SUBROUTINE BCONPL WRITTEN BY ROB BRACKEN, USGS, 13MAY87
C
C
      subroutine bconpl(iunit,ist,nparm,iparm,x,y,conlim,
     1 bspace,glen,gspace)
      save
C
C VARIABLE DECLARATION
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 /labcom/ dum1(11),size,dum2(4)
      integer*4 iparm(1)
      real*4 x(1),y(1)
      parameter(ddval=1.1e+38,dval=1.0e+37,tol=1.e-3)
      parameter (lental=5000)
      logical*4 post,post2
      common/captur/icapt
C
C INITIALIZE VARIABLES
C
C --- LENGTH OF PREVIOUS SEGMENT ---
      tlen=0.
C --- BEGINNING POINT FOR CURRENT SEGMENT ---
      x0=ddval
      y0=ddval
C --- BEGINNING INDEX OF NEXT SEGMENT ---
      ib=1
C --- ENDING INDEX OF PREVIOUS SEGMENT ---
      ie=ib-1
C
C READ POINTERS FROM IUNIT
C
      do 30 i=1,ist
        read(iunit,rec=i) nsrec,jsrec,contlb,linetp,post,size2
        if(nsrec.lt.0) then
C ---     HACHURES ARE TO BE MADE ---
          nsrec=-1*nsrec
          hspac=gspace
        else
C ---     HACHURES ARE NOT TO BE MADE ---
          hspac=0.
        end if
C ---   INITIALIZE ARRAY POINTER FOR SECTION OF SEGMENTS ---
        iptr=nsrec
C
C EXAMINE CONTOUR SEGMENTS
C
    5   if(iptr.ge.jsrec) go to 30
          call megard(iptr,nparm,iparm)
C ---     SET POINTER FOR NEXT SEGMENT ---
          jptr=iparm(1)
          isegtp=iparm(4)
C ---     SKIP UNACCEPTABLE SEGMENTS ---
          if(isegtp.ge.30) goto 101
C ---     READ IN ACCEPTABLE SEGMENTS ---
          iptr2=iptr+nparm
          npts=iparm(2)
          if((ib+npts-1).gt.lental) then
c         check for number of points in contour segment
c
          write(*,900) lental
  900     format(' number of points in contour line greater',
     1    ' than maximum(',i5,') allowed')
          write(*,950)
  950     format(' you can increase grid interval or turn off splining',
     1    /,' or eliminate hachers by setting mxhach,conlim,delb',
     2    ' & hachsp = 0')
          return
          endif
          call rmegar(iptr2,npts,x(ib))
          call rmegar(iptr2+npts,npts,y(ib))
          ic=iparm(3)
          if(ic.eq.0.and.
     1     (abs(x(ib)-x0).gt.tol.or.abs(y(ib)-y0).gt.tol)) then
c          if(ic.eq.0.and.(x(ib).ne.x0.or.y(ib).ne.y0)) then
C ---       SEGMENT DOES NOT JOIN WITH PREVIOUS SEGMENT ---
            if(tlen.gt.conlim) then
              if(icapt.ne.0) then
              write(23,1200) cont
 1200         format('contour value = ',f15.7)
              write(22,1000) cont,ie
 1000         format('contour value = ',f15.7,
     &        ' number of points = ',i10)
              write(*,1050) cont,ie
 1050         format(' writing contour value = ',f15.7,
     &        ' number of points = ',i10)
              write(22,1100) (x(k),y(k),k=1,ie)
 1100         format(8f10.3)
C ---         PLOT PREVIOUS SEGMENT ---
              else
              call dconlb(x,y,ie,post2,bspace,hspac)
              endif
            end if
C ---       RE-INITIALIZE TOTAL LENGTH ---
            tlen=0.
C ---       RE-INITIALZE BEGINNING POINT ---
            x0=x(ib)
            y0=y(ib)
C ---       SLIDE CURRENT SEGMENT TO BEGINNING OF X,Y ARRAY ---
            if(ie.gt.0) then
              do 10 j=1,npts
                x(j)=x(j+ie)
                y(j)=y(j+ie)
   10         continue
            end if
C ---       RE-INITIALIZE BEGINNING AND ENDING INDICES ---
            ib=1
            ie=ib-1
C ---       SET CONTOUR LEVEL, LINE TYPE, AND HACHURE PARMS ---
            size=size2
            call setcon(contlb,linetp,glen,hspac)
            post2=post
          end if
          if(x0.ge.dval.and.y0.ge.dval) then
C ---       RE-INITIALIZE BEGINNING POINT IF FIRST SEGMENT ---
            x0=x(1)
            y0=y(1)
          end if
C ---     FIND LENGTH OF CURRENT SEGMENT ---
          do 20 k=ib,ie+npts
            if(tlen.gt.conlim) goto 102
            dx=x(k)-x0
            dy=y(k)-y0
            ds=sqrt(dx*dx+dy*dy)
            tlen=tlen+ds
            x0=x(k)
            y0=y(k)
   20     continue
C ---     INCREMENT BEGINNING AND ENDING INDICES ---
  102     ib=ib+npts
          ie=ib-1
C ---     SET X0 AND Y0 TO END OF CURRENT SEGMENT ---
          x0=x(ie)
          y0=y(ie)
C ---     MOVE POINTER TO NEXT SEGMENT ---
  101     iptr=jptr
        go to 5
   30  continue
C --- PLOT LAST SEGMENT ---
      if(tlen.gt.conlim) then
              if(icapt.ne.0) then
              write(23,1200) cont
              write(22,1000) cont,ie
              write(*,1050) cont,ie
              write(22,1100) (x(k),y(k),k=1,ie)
              else
              call dconlb(x,y,ie,post2,bspace,hspac)
              endif
      end if
C
C RETURN TO CALLING ROUTINE
C
      return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  D C O N L B
C________________________________________________________________
C
C SUBROUTINE DCONLB IS A DRIVER FOR CONLAB.  IT USES CONLAB TO
C PUT LABELS ON THE CONTOURS AT INTERVALS ACCORDING TO PARAMETER
C BSPACE.
C
C   X,Y    - REAL*4.  ARRAY CONTAINING THE X & Y LOCATIONS.
C   NPTS   - INTEGER*4.  TOTAL NUMBER OF POINTS IN X & Y.
C   POST   - LOGICAL*4.  IF POST IS TRUE CONTOURS WILL BE LABELED
C            OTHERWISE THEY WILL NOT.
C   BSPACE - REAL*4.  DESIRED DISTANCE BETWEEN CONTOUR LABELS.
C   HSPACE - REAL*4.  EXACT DISTANCE BETWEEN HACHURE TICS.
C
C *** NOTE - ENTRY SETCONLAB MUST BE CALLED BEFORE DCONLAB ***
C
C SUBROUTINE DCONLB WRITTEN BY ROB BRACKEN, USGS, 19MAY87.
C
C
      subroutine dconlb(x,y,npts,post,bspace,hspace)
      save
c
      common /labcom/ dum1(12),wdist,dum2(3)
      real*4 x(1),y(1)
      logical*4 post,post2
      data distl2 /0./
c
      ibeg=1
      icc=0
      post2=post
      distl=distl2*distl2
      sdist=0.
      hdist=hspace*.5
c
  101 call conlab(x,y,ibeg,npts,icc,post2,distl,dista,ie,
     1 sdist,hdist)
      if(post2) then
        ibeg=ie
        post2=post
        distl=bspace*bspace
        goto 101
      else if(post) then
        if(distl2.lt.bspace) then
          distl2=distl2+1.5*wdist
        else
          distl2=0.
        end if
      end if
      return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  C O N L A B
C________________________________________________________________
C
C SUBROUTINE CONLAB IS A MODIFIED VERSION OF SUBROUTINE
C LABEL (IN CONTOUR).  THE MODIFICATIONS INCLUDE:  A) PLOTTING
C ONLY THE FIRST CONTOUR LINE BEFORE THE LABEL, B) CALLING A
C HACHURING ROUTINE TO PUT HACHURE TICKS ON THE PLOTTED LINE,  C)
C RETURNING THE LENGTH ALONG THE FIRST LINE FROM THE BEGINNING
C TO THE RIGHT SIDE OF THE LABEL, AND  D) ADDING THE EXTRA POINT
C TO THE ARRAY BEFORE THE AFTER LABEL LINE AND RETURNING ITS
C LOCATION.
C
C   X,Y    - REAL*4.  ARRAYS OF CONTOUR LINE X & Y LOCATIONS.
C   IBEG   - INTEGER*4.  THE POSITION (INDEX) IN X & Y OF THE
C            FIRST POINT TO BE CONSIDERED IN THE LABELING
C            PROCESS.  ALL PRECEDING POINTS WILL BE IGNORED.
C   NPTS   - INTEGER*4.  TOTAL NUMBER OF POINTS IN X & Y.
C   ICC    - INTEGER*4.  LINE CONTINUATION.  ICC=0 IS NO
C            CONTINUATION FROM LAST LINE.  ICC=1 IS CONTINUATION.
C            ALWAYS SET ICC TO 0 AND SUPPLY THE BEGINNING POINT
C            IN THE X & Y ARRAYS IF CONTINUATION IS DESIRED.
C   POST   - LOGICAL*4.  WHEN CALLING CONLAB IF POST IS TRUE
C            LABELS WILL BE PLOTTED.  IF POST IS FALSE LABELS
C            WILL NOT BE PLOTTED.  CONLAB SETS POST TO TRUE IF A
C            LABEL WAS ACTUALLY PLOTTED AND FALSE IF IT WAS NOT.
C   DISTL  - REAL*4.  DISTANCE TO FOLLOW CONTOUR LINE BEFORE
C            STARTING TO LOOK FOR A LABEL LOCATION.
C   DISTA  - REAL*4.  ACTUAL DISTANCE FROM THE BEGINNING OF THE
C            CONTOUR LINE TO THE FARTHEST RIGHT SIDE OF THE LABEL
C   IE     - INTEGER*4.  POSITION (INDEX) IN X & Y OF THE FIRST
C            POINT AFTER THE LABEL.  WARNING -- CONPLOT CHANGES
C            THE X & Y ARRAYS TO ADD AN EXTRA POINT IMMEDIATLY
C            AFTER THE LABEL.  THE NEW POINT IS IN POSITION IE.
C   SDIST  - REAL*4.  THIS PERTAINS ONLY TO THE HACHURING OPTION.
C            IT IS THE CURRENT DISTANCE FROM AN ARBITRARY
C            BEGINNING POSITION ALONG THE CONTOUR LINE.  IT WILL
C            BE ADDED TO WHEN HACHURE TICS ARE MADE.
C   HDIST  - REAL*4.  THIS PERTAINS ONLY TO THE HACHURING OPTION.
C            IT IS THE CURRENT DISTANCE (FROM THE SAME ARBITRARY
C            BEGINNING POSITION AS SDIST) TO THE NEXT HACHURE
C            TIC.  HDIST MUST BE GREATER THAN OR EQUAL TO SDIST.
C            IT WILL BE ADDED TO WHEN HACHURE TICS ARE MADE.
C
C *** NOTE -- CALL ENTRY SETLAB BEFORE CALLING CONLAB. ***
C
C USGS, 19MAY87
C
C
      subroutine conlab(x,y,ibeg,npts,icc,post,distl,dista,
     1 ie,sdist,hdist)
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/pltcom/npens
c
      dimension x(1),y(1)
      character fmtc*16,chars*24
      logical nolab,post,prime
      equivalence (f,sum),(isa,ya),(lchars,chars)
c
c setup for H-P plotters
c
      if(iplotr.eq.5) then
      ipen=linet/100
      if(ipen.lt.1) ipen=1
      if(ipen.gt.npens) ipen=npens
      call newpen(ipen)
      endif
c
c setup and scan
c
  100 is0=ibeg
      dista=0.
      nolab=.not.(post.and.nchar.gt.0.and.size.gt.0.)
      if (nolab) go to 50
c
c      do  5 is0=is0,npts-1
c        if(dista.ge.distl) goto 101
c        dx=x(is0+1)-x(is0)
c        dy=y(is0+1)-y(is0)
c        dista=dista+sqrt(dx*dx+dy*dy)
c    5 continue
c      goto 50
c
  101 is=is0+1
      isa=is+1
      sum=0.0
      in1=is
      ie=npts-1
      if (ie.lt.isa) 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
      if((xb-x(ibeg))**2+(yb-y(ibeg))**2.lt.distl) goto 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
      if(ibeg.eq.1.and.distl.gt.0) then
        distl=0.
        goto 100
      end if
c
c  can't fine spot, continue line without labelling.
   50 post=.false.
      npts2=npts-ibeg+1
C
C................................................................
C
C     CALL TO PLOT SYSTEM -- LINE PLOTTING ROUTINE
C
      ibeg2=ibeg
      ic=icc
   51 if(ibeg2+100.lt.npts) then
        call line(x(ibeg2),y(ibeg2),101,ic,linet)
        ibeg2=ibeg2+100
        ic=0
        goto 51
      else
        call line(x(ibeg2),y(ibeg2),npts-ibeg2+1,ic,linet)
      end if
C
C................................................................
C
      if(hspace.gt.1.e-30)
     1 call hach(x(ibeg),y(ibeg),npts2,sdist,hdist,hlen,hspace,
     2 linet)
c      do i=is0,npts-1
c        dx=x(i+1)-x(i)
c        dy=y(i+1)-y(i)
c        dista=dista+sqrt(dx*dx+dy*dy)
c      end do
      ie=npts
      go to 70
c
c  plottable location.
c
   60 post=.true.
      npts2=is-ibeg+1
C
C................................................................
C
C     CALL TO PLOT SYSTEM -- LINE PLOTTING ROUTINE
C
      ibeg2=ibeg
      ic=icc
   61 if(ibeg2+100.lt.is) then
        call line(x(ibeg2),y(ibeg2),101,ic,linet)
        ibeg2=ibeg2+100
        ic=0
        goto 61
      else
        call line(x(ibeg2),y(ibeg2),is-ibeg2+1,ic,linet)
      end if
C
C................................................................
C
      if(hspace.gt.1.e-30)
     1 call hach(x(ibeg),y(ibeg),npts2,sdist,hdist,hlen,
     2 hspace,linet)
      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)
      IF (IPLOTR .EQ. 5) THEN
                call VCHAR(0.5*(xa+xb),0.5*(ya+yb),lchars,nchar,
     1                     2,size,phi,ccor,0.)
      ELSE
                ilinet = 2 +(linet-mod(linet,100))
                call VCHAR(0.5*(xa+xb),0.5*(ya+yb),lchars,nchar,
     1                     ilinet,size,phi,ccor,0.)
      ENDIF
c
c      do i=is0,is
c        dx=x(i+1)-x(i)
c        dy=y(i+1)-y(i)
c        dista=dista+sqrt(dx*dx+dy*dy)
c      end do
c
c      dista=dista+sqrt((xa-x(is))**2+(ya-y(is))**2)
      ie=ie-1
      x(ie)=xa
      y(ie)=ya
c
   70 return
C
C________________________________________________________________
C
C     ENTRY  S E T C O N
C________________________________________________________________
C
c
c  entry to establish labeling character string.
c  must be called before CONLAB.
C
C   CONTLB - REAL*4.  VALUE OF CURRENT CONTOUR LEVEL.
C   LINETP - INTEGER*4.  NUMBER OF CURRENT LINE TYPE.
C   HLENG  - REAL*4.  LENGTH OF HACHURE TICS.
C   HSPAC  - REAL*4.  SPACING BETWEEN HACHURE TICS.  SETTING
C            HSPACE TO 0.E0 ELIMINATES HACHURE TIC PLOTTING.
C
      entry setcon(contlb,linetp,hleng,hspac)
C
      cont=contlb
      linet=linetp
      hlen=hleng
      hspace=hspac
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
  199 return
      end
C****************************************************************
C
C________________________________________________________________
C
C     FUNCTION  I N S I D E
C________________________________________________________________
C
C FUNCTION INSIDE IS A LOGICAL FUNCTION WHICH IS TRUE WHEN A TEST
C POINT IS WITHIN THE BOUNDARIES OF A POLYGON AND FALSE WHEN THE
C TEST POINT IS OUTSIDE THE POLYGON. A POINT ON THE BOUNDARY IS
C CONSIDERED INSIDE.
C
C   X      - REAL*4 ARRAY.  X COORDINATES OF A POLYGON.
C   Y      - REAL*4 ARRAY.  Y COORDINATES
C   NPTS   - INTEGER*4.  NUMBER OF POINTS IN THE POLYGON.
C            (NOTE: THE FIRST AND LAST POINT OF THE POLYGON
C                   DO NOT NEED TO HAVE THE SAME COORDINATES.)
C   XTEST  - REAL*4.  X COORDINATE OF THE TEST POINT.
C   YTEST  - REAL*4.  Y COORDINATE
C
C FUNCTION INSIDE WRITTEN BY ROB BRACKEN, USGS, 31MAR87
C
C
       function inside(x,y,npts,xtest,ytest)
C
C VARIABLE DECLARATION
C
      logical*4  inedge
      real*4 x(1),y(1)
C
C INITIALIZE TEST POINT OUTSIDE POLYGON
C
      inside=0
C
C INITIALIZE EDGE MARKER OFF OF POLYGON BORDER LINE
C
      inedge=.false.
C
C INITIALIZE X-AXIS CROSSING COUNTER
C
      ncross=0
C
C INITIALIZE I (CURRENT POINTER) AND J (PREVIOUS POINTER)
C
      i=1
      j=npts
C
C POLYGON ENTRY PROCEDURE
C
      if(y(j)-ytest) 101,102,103
C --- POINT J BELOW X-AXIS ---
  101 if(y(i)-ytest) 111,130,120
C --- POINT J ON X-AXIS ---
  102 j=j-1
      if(j.lt.1) goto 130
      if(y(j)-ytest) 104,102,105
C --- X-AXIS APPROACHED FROM BELOW ---
  104 i=0
      if(x(npts)-xtest) 131,132,133
C --- X-AXIS APPROACHED FROM ABOVE ---
  105 i=0
      if(x(npts)-xtest) 141,142,143
C --- POINT J ABOVE X-AXIS ---
  103 if(y(i)-ytest) 110,140,121
C
C I BELOW X-AXIS. J ON OR ABOVE X-AXIS.
C
  110 if(x(i).gt.xtest) then
        if(x(j).ge.xtest) then
C ---     SEGMENT CROSSES POSITIVE X-AXIS ---
          ncross=ncross+1
          goto 111
        end if
      else
        if(x(j).lt.xtest) then
C ---     SEGMENT CROSSES NEGATIVE X-AXIS ---
          goto 111
        end if
      end if
C --- SEGMENT CROSSES X-AXIS VERY NEAR TEST POINT ---
C --- VECTOR PRODUCT:  POINT BELOW X POINT ABOVE ---
      vprod=(x(i)-xtest)*(y(j)-ytest)-(y(i)-ytest)*(x(j)-xtest)
      if(vprod.gt.0.) then
C ---   SEGMENT CROSSES POSITIVE X-AXIS ---
        ncross=ncross+1
      else if(vprod.eq.0.) then
C ---   SEGMENT CROSSES TEST POINT ---
        inedge=.true.
      end if
C --- FOLLOW POLYGON EDGE ---
  111 do 115 i=i+1,npts
        if(y(i).ge.ytest) then
          if(y(i).eq.ytest) goto 130
C ---     POINT I IS ABOVE X-AXIS ---
          j=i-1
          goto 120
        end if
  115 continue
      goto 190
C
C I ABOVE X-AXIS. J ON OR BELOW X-AXIS.
C
  120 if(x(j).gt.xtest) then
        if(x(i).ge.xtest) then
C ---     SEGMENT CROSSES POSITIVE X-AXIS ---
          ncross=ncross+1
          goto 121
        end if
      else
        if(x(i).lt.xtest) then
C ---     SEGMENT CROSSES NEGATIVE X-AXIS ---
          goto 121
        end if
      end if
C --- SEGMENT CROSSES X-AXIS VERY NEAR TEST POINT ---
C --- VECTOR PRODUCT:  POINT BELOW X POINT ABOVE ---
      vprod=(x(j)-xtest)*(y(i)-ytest)-(y(j)-ytest)*(x(i)-xtest)
      if(vprod.gt.0.) then
C ---   SEGMENT CROSSES POSITIVE X-AXIS ---
        ncross=ncross+1
      else if(vprod.eq.0.) then
C ---   SEGMENT CROSSES TEST POINT ---
        inedge=.true.
      end if
C --- FOLLOW POLYGON EDGE ---
  121 do 125 i=i+1,npts
        if(y(i).le.ytest) then
          if(y(i).eq.ytest) goto 140
C ---     POINT I IS BELOW X-AXIS ---
          j=i-1
          goto 110
        end if
  125 continue
      goto 190
C
C I ON X-AXIS.  J BELOW X-AXIS.
C
  130 if(x(i)-xtest) 131,132,133
C --- SEGMENT CROSSES TESTPOINT ---
  132 inedge=.true.
C --- POINT I ON NEGATIVE X-AXIS ---
  131 i=i+1
      if(i.gt.npts) goto 190
      if(y(i)-ytest) 111,135,121
C --- NEW POINT I ON X-AXIS ---
  135 if(x(i)-xtest) 131,132,134
C --- SEGMENT CROSSES TEST POINT ---
  134 inedge=.true.
C --- POINT I ON POSITIVE X-AXIS ---
  133 i=i+1
      if(i.gt.npts) goto 190
      if(y(i)-ytest) 111,136,137
C --- NEW POINT I ON X-AXIS ---
  136 if(x(i)-xtest) 132,132,133
C --- NEW POINT I CROSSES POSITIVE X-AXIS ---
  137 ncross=ncross+1
      goto 121
C
C I ON X-AXIS.  J ABOVE X-AXIS.
C
  140 if(x(i)-xtest) 141,142,143
C --- SEGMENT CROSSES TESTPOINT ---
  142 inedge=.true.
C --- POINT I ON NEGATIVE X-AXIS ---
  141 i=i+1
      if(i.gt.npts) goto 190
      if(y(i)-ytest) 111,145,121
C --- NEW POINT I ON X-AXIS ---
  145 if(x(i)-xtest) 141,142,144
C --- SEGMENT CROSSES TEST POINT ---
  144 inedge=.true.
C --- POINT I ON POSITIVE X-AXIS ---
  143 i=i+1
      if(i.gt.npts) goto 190
      if(y(i)-ytest) 147,146,121
C --- NEW POINT I ON X-AXIS ---
  146 if(x(i)-xtest) 142,142,143
C --- NEW POINT I CROSSES POSITIVE X-AXIS ---
  147 ncross=ncross+1
      goto 111
C
C EVALUATE INEDGE AND NCROSS
C
  190 if(inedge) then
        inside=1
C --- INSIDE IF NCROSS IS ODD ---
      else if(iand(ncross,1).eq.1) then
        inside=1
      end if
C
C RETURN TO CALLING ROUTINE
C
      return
      end

C****************************************************************
C
C     SUBROUTINE SERIES  M E G A R R A Y
C________________________________________________________________
C
C SUBROUTINE SERIES MEGARRAY IS A BLOCK OF SUBROUTINES WHICH WORK
C TOGETHER TO MAKE A DIRECT ACCESS FILE WHICH "LOOKS" LIKE AN
C ENDLESS SINGLY DIMENSIONED ARRAY.
C
C SUBROUTINES:
C
C   MEGAOPEN  - TO BE CALLED WITH THE UNIT NUMBER BEFORE DOING
C               ANY MEGARRAY OPERATIONS.
C   MEGACLOSE - TO BE CALLED AFTER MEGARRAY OPERATIONS ARE
C               COMPLETED.
C   MEGAWRITE - CALL WHEN WRITING integers TO THE ARRAY.
C   MEGAREAD  - CALL WHEN READING integers FROM THE ARRAY.
C   rMEGAWRITE - CALL WHEN WRITING reals TO THE ARRAY.
C   rMEGAREAD  - CALL WHEN READING reals FROM THE ARRAY.
C
C PARAMETERS:
C
C   IUNIT  - INTEGER*4.  UNIT NUMBER TO BE USED FOR ARRAY STORAGE
C   NBEGIN - INTEGER*4.  BEGINNING LOCATION (IN WORDS) IN THE
C            MEGARRAY FOR THE SMALL ARRAY THAT IS BEING STORED OR
C            RETRIEVED.  THIS MAY BE ANY NUMBER WHICH WILL FIT ON
C            THE UNIT BEING ACCESSED.
C   LENGTH - INTEGER*4.  LENGTH (IN WORDS) OF THE SMALL ARRAY.
C   ARRAY  - INTEGER*4.  A SMALL ARRAY WHICH HAS A BEGINNING
C            LOCATION AT NBEGIN WITHIN THE MEGARRAY AND A SIZE OF
C            LENGTH.  THE SMALL ARRAY CAN BE ANY SIZE.
C
C SUBROUTINE SERIES MEGARRAY WRITTEN BY ROB BRACKEN, USGS, 9APR87
C
C________________________________________________________________
C
C     SUBROUTINE  M E G A O P
C________________________________________________________________
C
C
      subroutine megaop(iunit)
C
C VARIABLE DECLARATION
C
      parameter(lenbuf=450,nbuf=8)
      integer*4 junit,nummax,jboff(nbuf),jbnum(nbuf),jbmod(nbuf)
      INTEGER*4 buffer(lenbuf,nbuf)
      common /bufcom/ junit,nummax,jboff,jbnum,jbmod,buffer
C
C OPEN SCRATCH FILE
C
      open(unit=iunit,form='unformatted',
     1     status='scratch',recl=lenbuf*4,
     2     access='direct')
C
C INITIALIZE VARIABLES
C
      junit=iunit
      nummax=0
      do 20 j=1,nbuf
        jboff(j)=-1
        jbnum(j)=j
        jbmod(j)=0
        do 10 i=1,lenbuf
          buffer(i,j)=0
   10   continue
   20 continue
C
C RETURN TO CALLING PROGRAM
C
  199 return
      end
C
C
C________________________________________________________________
C
C     SUBROUTINE  M E G A C L
C________________________________________________________________
C
C
      subroutine megacl
C
C VARIABLE DECLARATION
C
      parameter(lenbuf=450,nbuf=8)
      integer*4 junit,nummax,jboff(nbuf),jbnum(nbuf),jbmod(nbuf)
      INTEGER*4 buffer(lenbuf,nbuf)
      common /bufcom/ junit,nummax,jboff,jbnum,jbmod,buffer
C
C WRITE TO DISK ALL SPECIFIED BUFFERS
C
      do 10 i=1,nbuf
        if(jbmod(i).ne.0) call bfswap(junit,jboff(i)+1,0,nummax,
     1   buffer(1,jbnum(i)),lenbuf)
   10 continue
C
C CLOSE SCRATCH UNIT AND RETURN TO CALLING ROUTINE
C
      close(unit=junit)
  299 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE M E G A W R
C________________________________________________________________
C
C
      subroutine megawr(nbegin,length,array)
C
C VARIABLE DECLARATION
C
      INTEGER*4 array(1)
      parameter(lenbuf=450,nbuf=8)
      integer*4 junit,nummax,jboff(nbuf),jbnum(nbuf),jbmod(nbuf)
      INTEGER*4 buffer(lenbuf,nbuf)
      common /bufcom/ junit,nummax,jboff,jbnum,jbmod,buffer
C
C INITIALIZE OFFSET POINTERS
C
      inoff=nbegin-1
      iaoff=0
      iboff=(inoff)/lenbuf
      iwoff=inoff-iboff*lenbuf
C
C CHECK WHETHER FIRST BUFFER IS BEING REFERENCED
C
      if(iboff.eq.jboff(1)) then
        laspta=length-iaoff
        goto 301
      end if
C
C CHECK WHETHER 2ND - LAST BUFFERS ARE BEING REFERENCED
C
  303 laspta=length-iaoff
      do 10 i=2,nbuf
        if(iboff.eq.jboff(i)) goto 302
   10 continue
C
C OUTPUT LEAST ACTIVE BUFFER AND READ IN BUFFER BEING REFERENCED
C
      i=i-1
      if(jbmod(i).eq.0) then
        numout=0
      else
        numout=jboff(i)+1
        jbmod(i)=0
      end if
      if(iwoff.eq.0.and.laspta.ge.lenbuf) then
        numin=0
      else
        numin=iboff+1
      endif
      call bfswap(junit,numout,numin,nummax,buffer(1,jbnum(i)),
     1            lenbuf)
C
C ROTATE REFERENCED BUFFER TO FIRST POSITION
C
  302 ibnum=jbnum(i)
      ibmod=jbmod(i)
      do 20 j=i,2,-1
        j1=j-1
        jboff(j)=jboff(j1)
        jbnum(j)=jbnum(j1)
        jbmod(j)=jbmod(j1)
   20 continue
      jboff(1)=iboff
      jbnum(1)=ibnum
      jbmod(1)=ibmod
C
C BEGIN WRITING ARRAY INTO BUFFERS
C
  301 j=iwoff-iaoff
      jbmod(1)=-1
      if(iwoff+laspta.le.lenbuf) then
C
C REMAINDER OF ARRAY FITS ENTIRELY INTO PRESENT BUFFER
C
        do 30 i=iwoff+1,iwoff+laspta
          buffer(i,jbnum(1))=array(i-j)
   30   continue
      else
C
C REMAINDER OF ARRAY EXTENDS BEYOND PRESENT BUFFER
C
        do 40 i=iwoff+1,lenbuf
          buffer(i,jbnum(1))=array(i-j)
   40   continue
        iaoff=iaoff+lenbuf-iwoff
        iboff=iboff+1
        iwoff=0
        goto 303
      end if
C
C RETURN TO CALLING ROUTINE
C
  399 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE M E G A RD
C________________________________________________________________
C
C
      subroutine megard(nbegin,length,array)
C
C VARIABLE DECLARATION
C
      INTEGER*4 array(1)
      parameter(lenbuf=450,nbuf=8)
      integer*4 junit,nummax,jboff(nbuf),jbnum(nbuf),jbmod(nbuf)
      INTEGER*4 buffer(lenbuf,nbuf)
      common /bufcom/ junit,nummax,jboff,jbnum,jbmod,buffer
C
C INITIALIZE OFFSET POINTERS
C
      inoff=nbegin-1
      iaoff=0
      iboff=(inoff)/lenbuf
      iwoff=inoff-iboff*lenbuf
C
C CHECK WHETHER FIRST BUFFER IS BEING REFERENCED
C
      if(iboff.eq.jboff(1)) then
        laspta=length-iaoff
        goto 401
      end if
C
C CHECK WHETHER 2ND - LAST BUFFERS ARE BEING REFERENCED
C
  403 laspta=length-iaoff
      do 10 i=2,nbuf
        if(iboff.eq.jboff(i)) goto 402
   10 continue
C
C OUTPUT LEAST ACTIVE BUFFER AND READ IN BUFFER BEING REFERENCED
C
      i=i-1
      if(jbmod(i).eq.0) then
        numout=0
      else
        numout=jboff(i)+1
      end if
      numin=iboff+1
      call bfswap(junit,numout,numin,nummax,buffer(1,jbnum(i)),
     1            lenbuf)
C
C ROTATE REFERENCED BUFFER TO FIRST POSITION
C
  402 ibnum=jbnum(i)
      ibmod=jbmod(i)
      do 20 j=i,2,-1
        j1=j-1
        jboff(j)=jboff(j1)
        jbnum(j)=jbnum(j1)
        jbmod(j)=jbmod(j1)
   20 continue
      jboff(1)=iboff
      jbnum(1)=ibnum
      jbmod(1)=ibmod
C
C BEGIN WRITING BUFFERS INTO ARRAYS
C
  401 j=iwoff-iaoff
      if(iwoff+laspta.le.lenbuf) then
C
C REMAINDER OF ARRAY CONTAINED ENTIRELY BY PRESENT BUFFER
C
        do 30 i=iwoff+1,iwoff+laspta
          array(i-j)=buffer(i,jbnum(1))
   30   continue
      else
C
C REMAINDER OF ARRAY EXTENDS BEYOND PRESENT BUFFER
C
        do 40 i=iwoff+1,lenbuf
          array(i-j)=buffer(i,jbnum(1))
   40   continue
        iaoff=iaoff+lenbuf-iwoff
        iboff=iboff+1
        iwoff=0
        goto 403
      end if
C
C RETURN TO CALLING ROUTINE
C
  499 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  R M E G A W
C________________________________________________________________
C
C
      subroutine rmegaw(nbegin,length,array)
C
C VARIABLE DECLARATION
C
      REAL*4 array(1)
      parameter(lenbuf=450,nbuf=8)
      integer*4 junit,nummax,jboff(nbuf),jbnum(nbuf),jbmod(nbuf)
      INTEGER*4 buffer(lenbuf,nbuf)
      common /bufcom/ junit,nummax,jboff,jbnum,jbmod,buffer
      REAL*4 RUFFER(LENBUF,NBUF)
      EQUIVALENCE (BUFFER(1,1),RUFFER(1,1))
C
C INITIALIZE OFFSET POINTERS
C
      inoff=nbegin-1
      iaoff=0
      iboff=(inoff)/lenbuf
      iwoff=inoff-iboff*lenbuf
C
C CHECK WHETHER FIRST BUFFER IS BEING REFERENCED
C
      if(iboff.eq.jboff(1)) then
        laspta=length-iaoff
        goto 301
      end if
C
C CHECK WHETHER 2ND - LAST BUFFERS ARE BEING REFERENCED
C
  303 laspta=length-iaoff
      do 10 i=2,nbuf
        if(iboff.eq.jboff(i)) goto 302
   10 continue
C
C OUTPUT LEAST ACTIVE BUFFER AND READ IN BUFFER BEING REFERENCED
C
      i=i-1
      if(jbmod(i).eq.0) then
        numout=0
      else
        numout=jboff(i)+1
        jbmod(i)=0
      end if
      if(iwoff.eq.0.and.laspta.ge.lenbuf) then
        numin=0
      else
        numin=iboff+1
      endif
      call bfswap(junit,numout,numin,nummax,buffer(1,jbnum(i)),
     1            lenbuf)
C
C ROTATE REFERENCED BUFFER TO FIRST POSITION
C
  302 ibnum=jbnum(i)
      ibmod=jbmod(i)
      do 20 j=i,2,-1
        j1=j-1
        jboff(j)=jboff(j1)
        jbnum(j)=jbnum(j1)
        jbmod(j)=jbmod(j1)
   20 continue
      jboff(1)=iboff
      jbnum(1)=ibnum
      jbmod(1)=ibmod
C
C BEGIN WRITING ARRAY INTO BUFFERS
C
  301 j=iwoff-iaoff
      jbmod(1)=-1
      if(iwoff+laspta.le.lenbuf) then
C
C REMAINDER OF ARRAY FITS ENTIRELY INTO PRESENT BUFFER
C
        do 30 i=iwoff+1,iwoff+laspta
          Ruffer(i,jbnum(1))=array(i-j)
   30   continue
      else
C
C REMAINDER OF ARRAY EXTENDS BEYOND PRESENT BUFFER
C
        do 40 i=iwoff+1,lenbuf
          Ruffer(i,jbnum(1))=array(i-j)
   40   continue
        iaoff=iaoff+lenbuf-iwoff
        iboff=iboff+1
        iwoff=0
        goto 303
      end if
C
C RETURN TO CALLING ROUTINE
C
  399 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  R M E G A R
C________________________________________________________________
C
C
      subroutine rmegar(nbegin,length,array)
C
C VARIABLE DECLARATION
C
      REAL*4 array(1)
      parameter(lenbuf=450,nbuf=8)
      integer*4 junit,nummax,jboff(nbuf),jbnum(nbuf),jbmod(nbuf)
      INTEGER*4 buffer(lenbuf,nbuf)
      common /bufcom/ junit,nummax,jboff,jbnum,jbmod,buffer
      REAL*4 RUFFER(LENBUF,NBUF)
      EQUIVALENCE (BUFFER(1,1),RUFFER(1,1))
C
C INITIALIZE OFFSET POINTERS
C
      inoff=nbegin-1
      iaoff=0
      iboff=(inoff)/lenbuf
      iwoff=inoff-iboff*lenbuf
C
C CHECK WHETHER FIRST BUFFER IS BEING REFERENCED
C
      if(iboff.eq.jboff(1)) then
        laspta=length-iaoff
        goto 401
      end if
C
C CHECK WHETHER 2ND - LAST BUFFERS ARE BEING REFERENCED
C
  403 laspta=length-iaoff
      do 10 i=2,nbuf
        if(iboff.eq.jboff(i)) goto 402
   10 continue
C
C OUTPUT LEAST ACTIVE BUFFER AND READ IN BUFFER BEING REFERENCED
C
      i=i-1
      if(jbmod(i).eq.0) then
        numout=0
      else
        numout=jboff(i)+1
      end if
      numin=iboff+1
      call bfswap(junit,numout,numin,nummax,buffer(1,jbnum(i)),
     1            lenbuf)
C
C ROTATE REFERENCED BUFFER TO FIRST POSITION
C
  402 ibnum=jbnum(i)
      ibmod=jbmod(i)
      do 20 j=i,2,-1
        j1=j-1
        jboff(j)=jboff(j1)
        jbnum(j)=jbnum(j1)
        jbmod(j)=jbmod(j1)
   20 continue
      jboff(1)=iboff
      jbnum(1)=ibnum
      jbmod(1)=ibmod
C
C BEGIN WRITING BUFFERS INTO ARRAYS
C
  401 j=iwoff-iaoff
      if(iwoff+laspta.le.lenbuf) then
C
C REMAINDER OF ARRAY CONTAINED ENTIRELY BY PRESENT BUFFER
C
        do 30 i=iwoff+1,iwoff+laspta
          array(i-j)=Ruffer(i,jbnum(1))
   30   continue
      else
C
C REMAINDER OF ARRAY EXTENDS BEYOND PRESENT BUFFER
C
        do 40 i=iwoff+1,lenbuf
          array(i-j)=Ruffer(i,jbnum(1))
   40   continue
        iaoff=iaoff+lenbuf-iwoff
        iboff=iboff+1
        iwoff=0
        goto 403
      end if
C
C RETURN TO CALLING ROUTINE
C
  499 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  B F S W A P
C________________________________________________________________
C
C
      subroutine bfswap(iunit,numout,numin,nummax,buffer,lenbuf)
C
C VARIABLE DECLARATION
C
      INTEGER*4 buffer(lenbuf)
C
C WRITE OUT BUFFER AND INCREMENT MAX NUMBER OF RECORD COUNTER
C
      if(numout.gt.0) then
        write(iunit,rec=numout) buffer
        if(numout.gt.nummax) nummax=numout
      end if
C
C READ IN NEW BUFFER IF WITHIN RANGE OF RECORDS IN FILE
C
      if(numin.eq.0) goto 599
      if(numin.le.nummax) then
        read(iunit,rec=numin) buffer
C
C FILL BUFFER WITH ZEROS
C
      else
        do 10 i=1,lenbuf
          buffer(i)=0
   10   continue
      end if
C
C RETURN TO CALLING ROUTINE
C
  599 return
      end
C
C****************************************************************
