C________________________________________________________________
C
C     PROGRAM  PC  C O N T O U R
C________________________________________________________________
C
c
c  general contouring program
c
c  originally developed and coded by..
c    gerald ian evenden
c    u. s. geological survey
c    denver, colorado  80225
c
c  modified for the honeywell/multics computer
c   by r.h. godson - april,1977
c  station posting, hi-low symbols, lat-lon ticks, active
c  defaults
c   by mike webring, USGS
c  modified for the DEC/VAX by
c   Cindy Cooper, EDS and Ray Watts, USGS
c  converted to pc's with Microsort Fortran v4.01
c   by r.h. godson june,1988
c
c
c  grids of up to about 4500 columns can be contoured with nwork=16000.
c  in practice, 30 to 40 rows per tier is prefered so that 300
c  to 400 columns is a desirable limit.  the formula for rows
c  per tier (rpt) is ......
c     rpt=4*(nwork-ncval-ncol)/(5*ncol+4)
c  if larger grids are contoured, recompile with nwork and
c  dimension of work increased and appropriate [huge] inserts.
c
      common /titls/ title1,title2,title3,fmtx,fmty
      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/chksiz/conlim
      common /concom/ ncol,nrow,bmin,bmax,
     & grad,gflg,ijs,ije,ijsi,ijei,prime,ij4,cont,
     & fltmax,lmult(0:3),idashs,linet,sigma,nsig
      common /llp/ latm(3),latx(3),longm(3),longx(3),cm(3),
     &             baslat(3),
     & iproj,xxx(2),yyy(2),sizep,unit,ip,neat,tint,itpost,ibound,
     & sizet,phi1,phi2
      COMMON /LXCB/ LXFILE,LXFMT,LXPROJ,LXUNIT,LXTYPE
      common /vect/ ispost,lowhi,nid,ncharv,szpost,vmin,vmax,
     & lifmtv(14),lfmtv(4),ich(20),lchid(20),DGRAD,sizehi,szlab
      common /labcom/ lchars(6),lfmtc(4),nchar,size,ifill(3),
     1 delb
      common /contrc/ cmin,cmax,dcval,ncval,nsec,gradi,
     1                acdel(200),acgrad(200),jcdash(200),
     2                ncval2,acsize(200),acmin(200),acmax(200)
      common /contur/ acval(200)
      COMMON /HACOM/ HACHLN,HACHSP,HACHGP,HACHMX,HACHLM,HACHVB
      common/captur/icapt
c
      dimension work(16000),idim(2),cr(2),iwork(16000)
c      character pgm*8
      character*56 title1,title2,title3,LXFILE
      character*16 fmtc,fmtx,fmty,LXFMT
      CHARACTER*16 FMTC2
c      integer asciino(6)
      real latm,latx,longm,longx
      logical prime
      equivalence (lfmtc,fmtc),(work(1),iwork(1))
      data nwork/16000/,idim(1)/0/
      data iprint,igrid,ista/6,10,11/
c
c     write introduction
c
      call header
c
c
c  get control paramters
      call name(*350)
c
c  convert map scale to data units/inch
      if(mscale.gt.0) then
        if(unit.gt.2.5 .or. unit.lt.-.5)
     1   stop ' invalid ''UNIT'' parm'
        unitc=1.0
        if(unit.gt.0.5) unitc=39.370079
        if(unit.gt.1.5) unitc=39370.079
        xscale=float(mscale)/unitc
        yscale=xscale
      endif
c
c  given limits instead of a grid - post stations
      if(abs(xxx(1))+abs(xxx(2)).eq.0. .and.
     1 abs(yyy(1))+abs(yyy(2)).eq.0.) go to 31
        xx(1)=xxx(1)
        xx(2)=xxx(2)
        yy(1)=yyy(1)
        yy(2)=yyy(2)
      call setnic(ier)
      call scale(xxx,yyy,xp1,yp1,3,ier)
      if(ispost.ne.0) call vector(xxx,yyy,szlab)
      idim(2)=2
      call endpt(idim)
      close(ista)
      stop
c
c  default contour interval or labels
   31 if (ncval.eq.0 .and. dcval.eq.0.0) go to 35
      if(dcval.gt.0. .and. ncval.gt.0) ncval=0
      if(dcval.lt.0. .or. nchar.lt.0 .or. fmtc(1:6).eq.'      ')
     1   then
        cr(1)=cmin
        cr(2)=cmax
        if(ncval.gt.0) then
          cr(1)=acval(1)
          cr(2)=acval(ncval)
        else
          read(igrid) title1,ip1,ip2,nc,nr,nz,xo,dx,yo,dy
          na=nc*nz
          call gmax(igrid,fltmax,na,work,cr(1),cr(2))
        endif
        if(dcval.lt.0.) dcval=0.
        if(cr(1).eq.cr(2)) dcval=1.
        call setax(cr,dcval,20,NCHAR2,FMTC2)
      endif
C
C ADDITION TO PROGRAM MADE 17FEB87 TO ALLOW USER SELECTION OF
C CONTOUR LABEL FORMAT WHEN ACVALS ARE BEING USED
C
      IF(NCHAR.LT.0 .OR. FMTC(1:6).EQ.'      ') THEN
        NCHAR=NCHAR2
        FMTC=FMTC2
      END IF
c
   35 if(ncval.gt.0) then
        do 40 i=1,ncval
   40   work(i)=acval(i)
      endif
      nwkres=nwork
      if (nsec.lt.0) nsec=0
c      if (gradi.lt.0.) gradi=0.
      if (nchar.lt.0) nchar=0
      if (size.lt.0.) size=0.
      if (nsig.ne.0) nsig=4
      sigma=.001*sigma
c
c  contour values...
   50 if(dcval.gt.0.0 .or. ncval.gt.0 .or. ispost.ne.0) go to 52
      write(iprint,51)
   51 format(' %neither contours or station posting specified')
      stop
   52 if (ncval.gt.0) go to 80
c
c  incremental contours mode
      if (dcval.ge.0.0) go to 70
      write(iprint,60)
   60 format(' %dcval less than zero')
      stop
   70 ixad=1
      go to 120
c
c  specified contours mode.
   80 if (ncval.eq.1) go to 110
c
c  check ascendency.
      do 100 i=2,ncval
      if (work(i).gt.work(i-1)) go to 100
      IF(acgrad(1).NE.0) THEN
        IF(acgrad(I).NE.acgrad(I-1)) GOTO 100
      END IF
      write(iprint,90)
90    format(' %non-ascending specified contours')
      stop
  100 continue
  110 ixad=ncval+1
c
c  open grid file for preliminary check.
  120 nwkres=nwkres-ncval
      call openck(igrid,xo,yo,delx,dely,
     1 work(ixad),nwkres,iquad)
      if (ncol.eq.0) go to 380
      xxx(1)=xx(1)
      xxx(2)=xx(2)
      yyy(1)=yy(1)
      yyy(2)=yy(2)
c
c  slice up core
      if(iquad.eq.3) go to 150
      nrowt=4*(nwkres-ncol)/(5*ncol+4)
      if(nrowt.ge.2) go to 140
c
      write(iprint,130)
  130 format(' %insufficient memory')
      go to 340
c
  140 idad=ixad+ncol
      iyad=nwork-nrowt+1
      ifad=idad+nrowt*ncol
      go to 160
c
  150 nrowt=4*nwkres/(13*ncol)
      if(nrowt.lt.2) go to 120
      k=nrowt*ncol
      idad=ixad+k
      iyad=nwork-k+1
      ifad=idad+k
c
c  ok so far, scale and annotate
  160 call setnic(i)
      if(i.ne.0) go to 340
      if(iquad.eq.3) go to 200
      k=ixad+ncol-1
      if(delx.ne.0) go to 180
      do 170 i=ixad,k
  170 work(i)=work(i)*xxscal
      go to 200
  180 xo=xo*xxscal
      delx=delx*xxscal
      do 190 i=ixad,k
      work(i)=xo
  190 xo=xo+delx
c
  200 npass=(nrow+nrowt-3)/(nrowt-1)
      jres=1+nrow-nrowt-(npass-2)*(nrowt-1)
c
c   branch around contour trace
c
      IF(HACHLM.LT.0.) HACHLM=HACHSP*7.
      CALL hacher(npass,NROWT,YO*YYSCAL,DELY*YYSCAL,
     1 hachln,hachsp,HACHGP,HACHMX,HACHLM,HACHVB,CONLIM,delb)
c      if(dcval.eq.0. .and. ncval.eq.0) go to 205
      if(ncval.eq.1.and.acval(1).ge.1.e37) go to 205
      call conttr(work(ixad),work(iyad),work(idad),
     1 iwork(ifad),work,nrowt,ncol,xxscal,yyscal,gflg,
     1 yo*yyscal,dely*yyscal,iquad,npass,jres,igrid,
     1 ixad,iyad)
      CALL gather
c
c  rescale in grid data units
c
  205 call scale(xxx,yyy,xp1,yp1,3,ier)
      if(ispost.eq.0) go to 210
c
c  post stations
      call vector(xxx,yyy,szlab)
      close(ista)
210   if(lowhi.eq.0) go to 211
c
c  plot high and low symbols
      if(dcval.eq.0. .and. ncval.eq.0) go to 211
      IF(LOWHI.GE.0) THEN
        NOHI=0
      ELSE
        NOHI=1
        LOWHI=-1*LOWHI
      END IF
      nc2=ncol+2*lowhi
      iw=2*lowhi+1
      if(nc2*(iw+4)+iw .gt. nwork) go to 211
      rewind igrid
      call hilow(work,work(nc2+1),work(2*nc2+1),work(3*nc2+1),
     & work(4*nc2+1),work(nc2*(iw+4)+1),nc2,iw,sizehi,DGRAD,NOHI)
211   continue
c
c  close plot
      if(icapt.eq.0) then
      idim(2)=2
      call endpt(idim)
      else
c      call setdef
      call setmod(2)
      print *,'files contour.cor and contour.val created'
      endif
  340 close(igrid)
      go to 390
c
c  command file end trap
  350 write(iprint,360)
  360 format(' %odd eof on command segment')
  380 continue
  390 stop
      end
C
C________________________________________________________________
C
C     SUBROUTINE  N A M E
C________________________________________________________________
C
      subroutine name(*)
c
c     namelist input and common block initialization
c
      common /titls/ title1,title2,title3,fmtx,fmty
      common /setxy/ xx(2),yy(2),xp1(4),yp1(4),iplotr,sizel,
     1 ncharx,nchary,sizex,sizey,adelx,adely,pllx,plly,
     & lintx,linty,xscale,yscale,xxscal,yyscal,mscale,delx,dely
      common/chksiz/conlim
      common /concom/ ncol,nrow,bmin,bmax,
     1 grad,gflg,ijs,ije,ijsi,ijei,prime,ij4,cont,
     1 fltmax,lmult(0:3),idashs,linet,sigma,nsig
      common /llp/ latm(3),latx(3),longm(3),longx(3),cm(3),
     1             baslat(3),
     1 iproj,xxx(2),yyy(2),sizep,unit,ip,neat,tint,itpost,ibound,
     1 sizet,phi1,phi2
      COMMON /LXCB/ LXFILE,LXFMT,LXPROJ,LXUNIT,LXTYPE
      common/txcb/txfile,iverb
      common /vect/ ispost,lowhi,nid,ncharv,szpost,vmin,vmax,
     & lifmtv(14),lfmtv(4),ich(20),lchid(20),DGRAD,sizehi,szlab
      common /labcom/ lchars(6),fmtc,nchar,size,idum(3),delb
      common /contrc/ cmin,cmax,dcval,ncval,nsec,gradi,
     1                acdel(200),acgrad(200),jcdash(200),
     2                ncval2,acsize(200),acmin(200),acmax(200)
      common /contur/ acval(200)
      COMMON /HACOM/ HACHLN,HACHSP,HACHGP,HACHMX,HACHLM,HACHVB
      common /pltcom/ npens
      common/files/ifile,ifile2,acfile,mxhach
      common/captur/icapt
c
      character*56 title1,title2,title3,cfile,ifile,ifile2,
     1             iblank,ifmtv
      character*16 fmtc,fmtx,fmty,jblank,fmtv
      CHARACTER*56 LXFILE,ACFILE,TXFILE
      CHARACTER*16 LXFMT
      character*4 chid(20)
      character*8 p(88)
      double precision p1,p2
      real latm,latx,longm,longx
      logical prime
c      ,ONGRAD,ONDASH
      REAL*4 ACFARR(7,200),JCFARR(7,200)
      INTEGER*4 LNKACF(210),IKEY(2)
      character*20 fields(40)
      EQUIVALENCE (ACFARR,JCFARR)
      equivalence (lchid,chid),(lifmtv,ifmtv),(lfmtv,fmtv)
      data (p(i),i=1,41)/'ifile= ','ifile2= ','acfile= ',
     & 'lxfile= ','title2= ','title3= ','fmtc= ','fmtx= ','fmty= ',
     & 'fmtv= ','ifmtv= ','lxfmt= ','baslat= ','cm= ','latm= ',
     & 'latx= ','longm= ','longx= ','xxx= ','yyy= ','adelx= ',
     & 'adely= ','cmax= ','cmin= ','conlim= ','dcval= ','delb= ',
     & 'dgrad= ','gradi= ','hachgp= ','hachlm= ','hachln= ',
     & 'hachsp= ','ibound= ','idashs= ','iplotr= ','iproj= ',
     & 'ispost= ','itpost= ','lintx= ','linty= '/
c
      data (p(i),i=42,88)/'lowhi= ','lxproj= ','lxunit= ',
     & 'lxtype= ','mscale= ','mxhach= ','nchar= ','ncharv= ',
     & 'ncharx= ','nchary= ','ncval= ','neat= ','nid= ','nsec= ',
     & 'nsig= ','phi1= ','phi2= ','pllx= ','plly= ','sigma= ',
     & 'size= ','sizel= ','sizep= ','sizex= ','sizey= ',
     & 'szpost= ','tint= ','unit= ','vmax= ','vmin= ','xscale= ',
     & 'yscale= ','chid= ','ich= ','acdel= ','acgrad= ',
     & 'acval= ','jcdash= ','icapt= ','hachvb= ','sizehi= ',
     & 'szlab= ','acsize= ','acmin= ','acmax= ','txfile= ',
     & 'sizet= '/
c
      data itty,iprint,icmd,igrid,ista/5,6,9,10,11/
      data iblank/' '/,jblank/' '/
      data dval/1.e37/,ddval/1.7e38/
c
      lmult(0)=256*256*256
      lmult(1)=256*256
      lmult(2)=256
      lmult(3)=1
      ibound=999
      conlim=0.
      lowhi=0
      sizehi=.08
      DGRAD=0.
      ispost=0
      xxx(1)=0.
      xxx(2)=0.
      yyy(1)=0.
      yyy(2)=0.
      do 9 i=1,20
      ich(i)=4
    9 chid(i)=' '
      nid=0
      szlab=.08
      szpost=.1
      vmin=0.
      vmax=1.
      ncharv=0
      fmtv='(f7.2)'
      ifmtv=' '
      cmin=0.
      cmax=0.
      dcval=-1.
      nchar=-1
      ncval=0
      nsec=5
      sizel=0.
      gradi=30.
      acdel(1)=0.
      acgrad(1)=0.
      gflg=1.0e37
      fltmax=1.0e37
      iplotr=9
      idashs=0
      JCDASH(1)=999
      acsize(1)=-1.
      acmin(1)=ddval
      acmax(1)=-ddval
      ACFILE=IBLANK
      linet=0
      ncharx=0
      nchary=0
      lintx=5
      linty=5
      size=.1
      sizex=.1
      sizey=.1
      adelx=0.
      adely=0.
      pllx=0.
      plly=0.
      xscale=0.
      yscale=0.
      mscale=0
      xp1(2)=0.
      yp1(2)=0.
      sigma=7.3
      nsig=0
      delb=7.0
      sizep=.1
      sizet=.1
      phi1=33.
      phi2=45.
      iproj=999
      unit=2
      neat=0
      tint=0.
      itpost=2
      cm(1)=999.
      cm(2)=0.
      cm(3)=0.
      baslat(1)=999.
      baslat(2)=0.
      baslat(3)=0.
      do 10 i=1,3
      latm(i)=0.
      latx(i)=0.
      longm(i)=0.
   10 longx(i)=0.
      title2=iblank
      title3=iblank
      fmtc=jblank
      fmtx=jblank
      fmty=jblank
      npens=6
      iverb=0
      ifile=iblank
      ifile2=iblank
      HACHSP=.15
      HACHLN=.05
      HACHGP=0.
      MXHACH=0
      HACHLM=-1.
      HACHVB=0.
      TXFILE=IBLANK
      LXFILE=IBLANK
      LXFMT=jblank
      LXPROJ=999
      LXUNIT=2
      LXTYPE=0
      ICAPT=0
   19 write(iprint,20)
   20 format(' enter command filename :'$)
      read(itty,30,err=21) cfile
      go to 22
   21 print *,'file not found, try again'
      go to 19
   30 format(a56)
   22 if(cfile.eq.iblank) go to 40
      open(unit=icmd,file=cfile,form='formatted',
     1 status='old')
c     read(icmd,parms,end=350)
      call namel(icmd)
      close (icmd)
      HACHMX=FLOAT(MXHACH)
      if(iverb.eq.1) then
      write(*,1000) p(1),ifile,p(2),ifile2,p(3),acfile,p(4),lxfile,
     & p(87),txfile,
     & p(5),title2,p(6),title3,p(7),fmtc,p(8),fmtx,p(9),fmty,
     & p(10),fmtv,p(11),ifmtv,p(12),lxfmt
 1000 format(1x,a8,a20,5x,a8,a20,/,3(1x,a8,a15,2x),/,
     & 2(1x,a8,a56,/),
     & 1x,3(a8,a16,2x),/,1x,3(a8,a16,2x))
      write(*,1100) p(13),baslat,p(14),cm,p(15),latm,p(16),latx,
     & p(17),longm,p(18),longx,p(19),xxx,p(20),yyy
 1100 format(3(2(1x,a8,f5.0,1x,f3.0,1x,f3.0),/),
     & 2(1x,a8,2f10.2,2x))
      write(*,1200) p(21),adelx,p(22),adely,p(23),cmax,p(24),cmin,
     & p(25),conlim,p(26),dcval,p(27),delb,p(28),dgrad,p(29),gradi,
     & p(30),hachgp,p(31),hachlm,p(32),hachln,p(33),hachsp,
     & p(81),hachvb
 1200 format(4(1x,a8,f9.2))
      write(*,1400) p(34),ibound,p(80),icapt,p(35),idashs,p(36),iplotr,
     & p(37),iproj,p(38),ispost,p(39),itpost,p(40),lintx,
     & p(41),linty,p(42),lowhi,p(43),lxproj,p(44),lxunit,
     & p(45),lxtype,p(46),mscale,p(47),mxhach,p(48),nchar,
     & p(49),ncharv,p(50),ncharx,p(51),nchary,p(52),ncval,
     & p(53),neat,p(54),nid,p(55),nsec,p(56),nsig
 1400 format(6(1x,a8,i3,1x))
      write(*,1500) p(57),phi1,p(58),phi2,p(59),pllx,p(60),plly,
     & p(61),sigma,p(62),size,p(82),sizehi,p(63),sizel,p(64),sizep,
     & p(88),sizet,
     & p(65),sizex,p(66),sizey,p(83),szlab,p(67),szpost,p(68),tint,
     & p(69),unit,p(70),vmax,p(71),vmin,p(72),xscale,p(73),yscale
 1500 format(4(1x,a8,f10.2))
      pause 'Hit Return to Continue ;Otherwise Hit Cntrl C'
      if(chid(1).ne.' ') then
      write(*,1600) p(74),chid,p(75),ich
 1600 format(1x,a8,/2(10(2x,a4)),/,1x,a8,/2(10(2x,i4)),/)
      pause 'Hit Return to Continue ;Otherwise Hit Cntrl C'
      endif
      if(acdel(1).ne.0.) then
      write(*,1700) p(76),(acdel(i),i=1,ncval)
      pause 'Hit Return to Continue ;Otherwise Hit Cntrl C'
      endif
      if(acgrad(1).ne.0.) then
      write(*,1700) p(77),(acgrad(i),i=1,ncval)
      pause 'Hit Return to Continue ;Otherwise Hit Cntrl C'
      endif
      if(acval(1).ne.0.) then
      write(*,1700) p(78),(acval(i),i=1,ncval)
      pause 'Hit Return to Continue ;Otherwise Hit Cntrl C'
      endif
      if(acsize(1).ne.-1.) then
      write(*,1700) p(84),(acsize(i),i=1,ncval)
      pause 'Hit Return to Continue ;Otherwise Hit Cntrl C'
      endif
      if(acmin(1).ne.ddval) then
      write(*,1700) p(85),(acmin(i),i=1,ncval)
      pause 'Hit Return to Continue ;Otherwise Hit Cntrl C'
      endif
      if(acmax(1).ne.-ddval) then
      write(*,1700) p(86),(acmax(i),i=1,ncval)
      pause 'Hit Return to Continue ;Otherwise Hit Cntrl C'
      endif
      if(jcdash(1).ne.999) then
      write(*,1800) p(79),(jcdash(i),i=1,ncval)
      pause 'Hit Return to Continue ;Otherwise Hit Cntrl C'
      endif
      endif
 1700 format(1x,a8,/,7(f10.2))
 1800 format(1x,a8,/,15(i5))
C
C IF ACFILE IS SPECIFIED THEN THIS READS IT FREE FORMAT
C    FIELD 1 - ACVAL
C    FIELD 2 - THE CONTOUR INTERVAL APPLYING TO THE ACVAL
C              (THE MAXIMUM NUMBER OF CONTOUR LINES PER INCH
C               IS CALCULATED BASED ON THIS CONTOUR INTERVAL)
C    FIELD 3 - THE MAXIMUM GRADIENT FOR THE ACVAL IN UNITS OF
C              CONTOUR LINES PER INCH (CONTOUR INTERVAL FOR
C              CALCULATING THIS IS THE INTERVAL IN FIELD 2)
C    FIELD 4 - THE TYPE OF LINE TO BE DRAWN FOR THIS ACVAL
C    FIELD 5 - THE SIZE OF THE LABELS
C    FIELD 6 - THE MINIMUM CONTOUR LEVEL WHEN ACDEL IS NEGATIVE.
C    FIELD 7 - THE MAXIMUM CONOTUR LEVEL WHEN ACDEL IS NEGATIVE.
C
      if(acfile.eq.iblank) goto 35
      open(unit=icmd,file=acfile,status='old',form='formatted')
      open(24,status='scratch')
      ncval=0
   31 ncval=ncval+1
   36 call freerd(icmd,nfield,fields,0,*32,*32)
      if(nfield.lt.1) goto 36
      do 37 i=1,nfield
        k=ichar(fields(i)(1:1))
        if(k.lt.43.or.k.gt.57.or.k.eq.47) goto 38
   37 continue
   38 nfield=i-1
      if(nfield.lt.1) goto 36
      if(nfield.lt.7) then
        do 39 i=nfield+1,7
          fields(i)=iblank
   39   continue
      endif
      goto (71,72,73,74,75,76,77),nfield
   77 read(24,*) acval(ncval),acdel(ncval),acgrad(ncval),
     & jcdash(ncval),acsize(ncval),acmin(ncval),acmax(ncval)
      go to 711
   76 read(24,*) acval(ncval),acdel(ncval),acgrad(ncval),
     & jcdash(ncval),acsize(ncval),acmin(ncval)
      go to 711
   75 read(24,*) acval(ncval),acdel(ncval),acgrad(ncval),
     & jcdash(ncval),acsize(ncval)
      go to 711
   74 read(24,*) acval(ncval),acdel(ncval),acgrad(ncval),
     & jcdash(ncval)
      go to 711
   73 read(24,*) acval(ncval),acdel(ncval),acgrad(ncval)
      go to 711
   72 read(24,*) acval(ncval),acdel(ncval)
      go to 711
   71 read(24,*) acval(ncval)
  711 if(ncval.eq.1) then
        goto (51,52,53,54,55,56,57),nfield+1
   51   goto 58
   52   acdel(ncval)=1.e+7
   53   acgrad(ncval)=abs(gradi)
   54   jcdash(ncval)=0.
   55   acsize(ncval)=size
   56   acmin(ncval)=0.
   57   acmax(ncval)=0.
   58   goto 31
      endif
      goto (61,62,63,64,65,66,67),nfield+1
   61 goto 68
   62 acdel(ncval)=1.e+7
   63 acgrad(ncval)=abs(gradi)
   64 jcdash(ncval)=0.
   65 acsize(ncval)=size
   66 acmin(ncval)=0.
   67 acmax(ncval)=0.
   68 goto 31
C
C      IF(ACFILE.EQ.IBLANK) GOTO 35
C      OPEN(UNIT=ICMD,FILE=ACFILE,FORM='FORMATTED',STATUS='OLD',
C     1     READONLY,CARRIAGECONTROL='LIST')
C      NCVAL=0
C   31 NCVAL=NCVAL+1
C      READ(ICMD,*,END=32) ACVAL(NCVAL),ACDEL(NCVAL),
C     1 ACGRAD(NCVAL),JCDASH(NCVAL)
C      GOTO 31
C
   32 NCVAL=NCVAL-1
      CLOSE(UNIT=ICMD)
      close(24)
C
C SORT THE ACVAL, ACDEL, ACGRAD, JCDASH ,ACSIZE, ACMIN,
C AND ACMAX ARRAYS
C
   35 if(ncval.gt.0) then
        do 41 j=1,ncval
          if(acdel(j).lt.0) acval(j)=acdel(j)
   41   continue
        if(acdel(1).eq.0.) then
          do 42 j=1,ncval
            acdel(j)=1.e+7
   42   continue
        end if
        if(acgrad(1).eq.0.) then
          do 43 j=1,ncval
            acgrad(j)=abs(gradi)
   43   continue
        end if
        do 44 j=1,ncval
          acgrad(j)=acdel(j)*abs(acgrad(j))
   44   continue
        if(jcdash(1).eq.999) then
          do 45 j=1,ncval
            jcdash(j)=iabs(idashs)
   45   continue
        end if
        if(acsize(1).lt.0.) then
          do 46 j=1,ncval
            acsize(j)=size
   46   continue
        end if
        if(acmin(1).gt.dval) then
          do 47 j=1,ncval
            acmin(j)=0.
   47   continue
        end if
        if(acmax(1).lt.-dval) then
          do 48 j=1,ncval
            acmax(j)=0.
   48   continue
        end if
        nkey=2
        ikey(1)=3
        ikey(2)=1
        do 49 j=1,ncval
          acfarr(1,j)=acval(j)
          acfarr(2,j)=acdel(j)
          acfarr(3,j)=acgrad(j)
          jcfarr(4,j)=iabs(jcdash(j))
          acfarr(5,j)=acsize(j)
          acfarr(6,j)=acmin(j)
          acfarr(7,j)=acmax(j)
   49   continue
        call sortil(acfarr,7,200,ncval,ikey,nkey,lnkacf,210,
     1   ier2)
        ll=lnkacf(ncval+1)
        do 50 j=1,ncval
          acval(j)=acfarr(1,ll)
          acdel(j)=acfarr(2,ll)
          acgrad(j)=acfarr(3,ll)
          jcdash(j)=jcfarr(4,ll)
          acsize(j)=acfarr(5,ll)
          acmin(j)=acfarr(6,ll)
          acmax(j)=acfarr(7,ll)
          ll=lnkacf(ll)
   50   continue
        do 81 j=1,ncval
          if(acdel(j).ge.0) goto 33
          acfarr(1,j)=acval(j)
          acfarr(2,j)=acdel(j)
          acfarr(3,j)=acgrad(j)
          acfarr(4,j)=jcdash(j)
          acfarr(5,j)=acsize(j)
          acfarr(6,j)=acmin(j)
          acfarr(7,j)=acmax(j)
   81   continue
   33   ncval2=j-1
        if(ncval2.le.0) goto 34
        nkey=1
        ikey(1)=2
        call sortil(acfarr,7,200,ncval2,ikey,nkey,lnkacf,210,
     1   ier2)
        ll=lnkacf(ncval2+1)
        do 82 j=1,ncval2
          acval(j)=acfarr(1,ll)
          acdel(j)=abs(acfarr(2,ll))
          acgrad(j)=abs(acfarr(3,ll))
          jcdash(j)=jcfarr(4,ll)
          acsize(j)=acfarr(5,ll)
          acmin(j)=acfarr(6,ll)
          acmax(j)=acfarr(7,ll)
          ll=lnkacf(ll)
   82   continue
   34   continue
      end if
C
   40 ip=iplotr
      p1=phi1*1.745329252d-2
      p2=phi2*1.745329252d-2
      if(iproj.eq.4) call setlam(p1,p2)
      nogrid=0
      if(abs(xxx(1))+abs(xxx(2)).ne.0. .and.
     1 abs(yyy(1))+abs(yyy(2)).ne.0.) nogrid=1
c
      if(nogrid.eq.0 .and. ifile.eq.iblank) then
        print 60
   60   format(' enter grid filename :'$)
        read(itty,69) ifile
   69   format(a56)
      endif
      if(nogrid.eq.0) open(unit=igrid,file=ifile,
     1                status='old',form='unformatted')
c
      if(ispost.eq.0) return
      if(ifile2.eq.iblank) then
        print 70
   70   format(' enter station filename :'$)
        read(itty,69) ifile2
      endif
      if(ifmtv.eq.iblank) then
        print 80
   80   format(' enter format, car ret if binary :'$)
        read(itty,69) ifmtv
      endif
      if(ifmtv.eq.iblank) then
        open(unit=ista,file=ifile2,form='unformatted',
     1  status='old')
      else
        open(unit=ista,file=ifile2,form='formatted',
     1  status='old')
      endif
      return
  350 close(icmd)
      return 1
      end
C
C________________________________________________________________
C
C     SUBROUTINE  O P E N C K
C________________________________________________________________
C
      subroutine openck(igrid,xo,yo,delx,dely,
     & xdata,nwork,iquad)
c
c  initialize grid input file.
c
      common /titls/ title1,title2,title3,fmtx,fmty
      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,delxx,delyy
      common /concom/ ncol,nrow,bmin,bmax,
     & grad,gflg,ijs,ije,ijsi,ijei,prime,ij4,cont,
     & fltmax,lmult(0:3),idashs,linet,sigma,nsig
      character*56 title1,title2,title3
      character*16 fmtx,fmty
      dimension xdata(1),ipgm(2)
      data iprint/6/
c
      read(igrid,end=90) title1,ipgm,ncol,nrow,iquad,xo,delx,yo,
     1                   dely
      delxx=delx
      delyy=dely
      if (ncol.lt.3.or.nrow.lt.3) go to 110
      if(iquad.eq.3.and.delx.eq.0..and.dely.eq.0.) go to 80
      if (ncol.gt.nwork) go to 150
      if (iquad.ne.1) go to 130
      if (delx.ne.0) go to 30
      read(igrid,end=90) (xdata(i),i=1,ncol)
      xx(1)=xdata(1)
      xx(2)=xdata(ncol)
      go to 40
   30 xx(1)=xo
      xx(2)=xo+(ncol-1)*delx
   40 if(dely.ne.0) go to 50
      call yread(igrid,yy,nrow)
      if(nrow.lt.3) go to 110
      go to 60
   50 yy(1)=yo
      yy(2)=yo+(nrow-1)*dely
   60 if(delx.ne.0..and.dely.ne.0.) return
   70 rewind igrid
      read(igrid)
      if(delx.eq.0.) read(igrid)
      return
   80 i=ncol*3
      if(i.gt.nwork) go to 150
      read(igrid)
      call xyread(igrid,xx,yy,xdata,i)
      if(nrow.lt.3) go to 110
      go to 70
   90 write(iprint,100)
  100 format(' %end of file while processing header of input',
     1       ' file')
      go to 170
  110 write(iprint,120)
  120 format(' %no. rows or columns less than 3')
      go to 170
  130 write(iprint,140)
  140 format(' %nz on header of standard file is greater than 1')
      go to 170
  150 write(iprint,160)
  160 format(' %core exceeded for x values')
c
  170 close (igrid)
      ncol=0
      return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  Y R E A D
C________________________________________________________________
C
      subroutine yread(igrid,yr,nrow)
c
c     scans specified row grid
c
      dimension yr(2)
      nrow=0
      read(igrid,end=20) yr(1)
   10 nrow=nrow+1
      read(igrid,end=20) a
      yr(2)=a
      go to 10
   20 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  X Y R E A D
C________________________________________________________________
C
      subroutine xyread(igrid,xr,yr,work,ncol3)
c
c     scans quadrilateral grid for x-y range
c
      dimension xr(2),yr(2),work(ncol3)
      common /concom/ ncol,nrow,bmin,bmax,
     & grad,gflg,ijs,ije,ijsi,ijei,prime,ij4,cont,
     & fltmax,lmult(0:3),idashs,linet,sigma,nsig
c
c     initialize
c
      nrow=0
      read(igrid,end=70) dummy,work
      xmin=fltmax
      xmax=-fltmax
      ymin=fltmax
      ymax=-fltmax
c
c scan loop
c
   10 nrow=nrow+1
      if(work(i+2).gt.gflg) go to 55
      do 50 i=1,ncol3,3
      if(work(i).le.xmax) go to 20
      xmax=work(i)
      go to 30
   20 if(work(i).ge.xmin) go to 30
      xmin=work(i)
   30 if(work(i+1).ge.ymin) go to 40
      ymin=work(i+1)
      go to 50
   40 if(work(i+1).le.ymax) go to 50
      ymax=work(i+1)
   50 continue
   55 continue
c
      read(igrid,end=60) dummy,work
      go to 10
c
   60 xr(1)=xmin
      xr(2)=xmax
      yr(1)=ymin
      yr(2)=ymax
c
   70 return
      end
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)
c
c     basic input and tier control
c
      dimension x(1),y(1),z(1),c(1)
      integer f(1)
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)
      CALL grdegg
  180 if(last.ne.0) return
  190 continue
      return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  C O N T R
C________________________________________________________________
C
      subroutine contr(grid,flags,work,acval,nrowf)
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
        PRIME=.TRUE.
        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)
  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)
          END IF
    6 continue
        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.
      do 40 i=ii,ncval
      cont=acval(i)
      if (cont.gt.bmax) go to 140
      IF(JCDASH(1).NE.999) LINET=IABS(JCDASH(I))
      call setlab
   40 call scan(grid,flags,work)
      go to 140
c
c  execution for delta contour levels.
   50 if (cmin.eq.0..and.cmax.eq.0.) go to 60
      icont=amax1(bmin,cmin)/dcval
      amx=amin1(bmax,cmax)
      go to 70
   60 icont=bmin/dcval
      amx=bmax
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 linet=iabs(mod(idashs,8))
      go to 135
   90 linet=0
      go to 135
  100 if (idashs) 110,130,120
  110 linet=iabs(idashs/8)*8
      go to 130
  120 linet=idashs
  130 call setlab
C
  135 call scan(grid,flags,work)
      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)
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)
      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)
      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)
      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 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.This was not a
c     serious problem and it was a random thing that
c     was dependent on the size of the work array.
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)
      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)
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)
      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
c        call bsplin(x,y,npts,xx,yy,nnt,step,sigma)
        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
c          call bsplin(x,y,npts,xx,yy,nnt,step,sigma)
          call bsplin(x,y,npts,xx,yy,nnt,ninc,sigma)
          call feed(xx,yy,nnt,iccc,-999,100+kk)
c          ninc=nint(1./step)
          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 = 2 +(linet-mod(linet,100))
                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
      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________________________________________________________________
C
C     SUBROUTINE  S E T A X
C________________________________________________________________
C
      subroutine setax(x,dx,maxint,nch,fmt)
c  adjust interval and labeling format
c  input x-min max range, optional dx-spacing and
c  maxint-intervals
c  returns dx, nch-number of significant figures, fmt-labeling
c  format
      dimension x(2)
      character fmt*(*)
      ixpn(r)=int(alog10(abs(r))+100.)-100
      fmt='(1pe13.6)'
      nch=13
      if(dx.ne.0.) go to 5
      if(maxint.le.0) maxint=20
      t=(x(2)-x(1))/float(maxint)
      if(t.lt.1.e-20) return
      p10=sign(1.0,t)*10.**ixpn(t)
      t1=t/p10
      if(t1.le.1.0) dx=p10
      if(t1.gt.1.0) dx=2.*p10
      if(t1.gt.2.0) dx=5.*p10
      if(t1.gt.5.0) dx=10.*p10
    5 idecm=0
      n10=ixpn(dx)
      if(n10.lt.0) idecm=iabs(n10)
      n10=ixpn( amax1(abs(x(1)),abs(x(2))) )
      iw=4+idecm
      if(n10.ge.0) iw=3+n10+idecm
      if(iw.gt.9) go to 11
      write(fmt,10) iw,idecm
   10 format('(f',i1,'.',i1,')' )
      nch=iw
      if(idecm.eq.0) nch=nch-1
   11 return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  G M A X
C________________________________________________________________
C
      subroutine gmax(igrid,fltmax,na,z,zmin,zmax)
      dimension ihdr(16),z(na)
      zmin=fltmax
      zmax=-fltmax
      rewind(igrid)
      read(igrid) ihdr,nc,nr,nz,xo,dx,yo,dy
      if(nz.gt.1) go to 30
      if(dx.eq.0.) read(igrid)
      do 20 j=1,nr
      read(igrid,end=30) y,z
      do 10 i=1,nc
      if(z(i).gt.fltmax) go to 10
      if(z(i).gt.zmax) zmax=z(i)
      if(z(i).lt.zmin) zmin=z(i)
   10 continue
   20 continue
   30 rewind(igrid)
      return
      end
C
C________________________________________________________________
C
C     SUBROUTINE  S E T N I C
C________________________________________________________________
C
      subroutine setnic(ier)
c
c  general setup of contour niceties.
c
      common /titls/ title1,title2,title3,fmtx,fmty
      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 /llp/ latm(3),latx(3),longm(3),longx(3),cm(3),
     1             baslat(3),
     & iproj,xxx(2),yyy(2),sizep,unit,ip,neat,tint,itpost,ibound,
     & sizet,phi1,phi2
      COMMON /LXCB/ LXFILE,LXFMT,LXPROJ,LXUNIT,LXTYPE
      common/txcb/txfile,iverb
      common/captur/icapt
c
      dimension idim(1),it1(14),it2(14),it3(14),ifmtx(4),ifmty(4)
      character*56 title1,title2,title3,blankt,LXFILE,txfile
      character*16 fmtx,fmty,blankf,LXFMT
      real latm,latx,longm,longx
      equivalence (it1,title1),(it2,title2),(it3,title3)
      equivalence (ifmtx,fmtx),(ifmty,fmty)
      data iprint/6/,idim(1)/0/,blankf/' '/,blankt/' '/
c
      tspace=0.
      if(title1.ne.blankt) tspace=tspace+1.
      if(title2.ne.blankt) tspace=tspace+1.
      if(title3.ne.blankt) tspace=tspace+1.
      if(adelx.eq.0.0 .and. ncharx.eq.0) go to 1
      if(fmtx.ne.blankf .and. ncharx.gt.0) go to 1
      call setax(xxx,adelx,20,ncharx,fmtx)
    1 if(adely.eq.0.0 .and. nchary.eq.0) go to 2
      if(fmty.ne.blankf .and. nchary.gt.0) go to 2
      call setax(yyy,adely,20,nchary,fmty)
    2 if (sizex.lt.0.) sizex=.08
      if (sizey.lt.0.) sizey=.08
      if (sizel.lt.0.) sizel=.08
      if (sizep.lt.0.) sizep=.08
      if (sizet.lt.0.) sizet=.08
      if (ncharx.lt.0) ncharx=0
      if (nchary.lt.0) nchary=0
      if (lintx.lt.1) lintx=5
      if (linty.lt.1) linty=5
      sizelb=1.5*sizel
      if(pllx.le.0. .and. (iproj.ne.0 .or. iproj.ne.999))
     1  pllx=12.*sizep
      if(plly.le.0. .and. (iproj.ne.0 .or. iproj.ne.999))
     1 plly=5.*sizep+sizelb*(tspace+2.)
      if(pllx.eq.0. .and. plly.eq.0.) sizep=0.
c
c  margin requirement setup.
c
   20 if (pllx.le.0.) go to 30
      xp1(3)=pllx
      go to 40
   30 xp1(3)=(ncharx+1)*sizex
   40 if (plly.le.0.) go to 50
      yp1(3)=plly
      go to 60
   50 yp1(3)=(nchary+1)*sizey+(tspace+1.)*sizelb
c
   60 dx=xx(2)-xx(1)
      dy=yy(2)-yy(1)
      call pltset(iplotr,xp1(4),yp1(4),idim)
c     if(iplotr.eq.5) write(iprint,65) xp1(4),yp1(4)
c  65 format(/,' plotter size(inches): x=',f7.2,'  y=',f7.2)
c
c  check scaling
c
      if (xscale.le.0.) go to 70
      if (yscale.le.0.) yscale=xscale
      go to 80
   70 if (yscale.le.0.) go to 90
      xscale=yscale
   80 xxscal=sign(1./xscale,dx)
      yyscal=sign(1./yscale,dy)
c
c  fixed scaling.
c
      xp1(1)=abs(dx*xxscal)
      yp1(1)=abs(dy*yyscal)
      xp1(4)=xp1(1)+xp1(3)+sizex+13.*sizep
      yp1(4)=yp1(1)+yp1(3)+sizey+3.*sizep
      go to 120
c
c  relative scaling, eg. xxscal=yyscal=0.
c
   90 xp1(4)=amin1(xp1(4),10.)
      yp1(4)=amin1(yp1(4),8.)
      xp1(1)=xp1(4)-xp1(3)-sizex-6.*sizep
      yp1(1)=yp1(4)-yp1(3)-sizey-3.*sizep
      if (xp1(1).gt.0..and.yp1(1).gt.0.) go to 110
      write(iprint,100)
  100 format(' %margin requires all plot area')
      ier=2
      go to 140
  110 xxscal=xp1(1)/dx
      yyscal=yp1(1)/dy
      xy=amin1(abs(xxscal),abs(yyscal))
      xxscal=sign(xy,xxscal)
      yyscal=sign(xxscal,yyscal)
      xp1(1)=abs(dx*xxscal)
      yp1(1)=abs(dy*yyscal)
c
c  initial labelling scale call
c
  120 call scale(xx,yy,xp1,yp1,4,ier)
      if (ier.ne.0) go to 140
c
c     don't plot if contour capturing is requested
c
      if(icapt.eq.0) then
c     if(iplotr.ne.1 .and. iplotr.ne.4) write(iprint,125) xp1(4),
c    1 yp1(4)
  125 format(/,' plot size (inches): x=',f7.2,'  y=',f7.2)
c
c  plot titles
c
      if (sizel.le.0.) go to 130
      xcent=0.5*xp1(1)+xp1(3)
      yt=sizelb
      if(title3.eq.blankt) go to 127
      n=leftj(title3)
      xt=abs(xcent-float(n/2)*.75*sizel)
      call vchar(xt,yt,it3,56,3,.75*sizel,0.,0.,0.)
      yt=sizelb+sizelb
127   if(title2.eq.blankt) go to 128
      n=leftj(title2)
      xt=abs(xcent-float(n/2)*.75*sizel)
      call vchar(xt,yt,it2,56,3,.75*sizel,0.,0.,0.)
      yt=yt+sizelb
128   if(title1.eq.blankt) go to 130
      n=leftj(title1)
      xt=abs(xcent-float(n/2)*sizel)
      call vchar(xt,yt,it1,56,3,sizel,0.,0.,0.)
c
c  plot axis
c
  130 if(adelx.eq.0.) go to 131
      call xaxis(xx,yy,xp1,adelx,lintx,sizex,ifmtx,ncharx)
  131 if(adely.eq.0.) go to 132
      call yaxis(yy,xx,yp1,adely,linty,sizey,ifmty,nchary)
  132 if(neat.eq.0) call neatl
c
c  plot lat-long
      if(iproj.gt.0 .and. iproj.lt.999) call llpost
c
c  subroutine state requires the world data bank 2
  135 if(ibound .ne. 999) call state
      IF(LXFILE.NE.BLANKT)
     1 CALL LXPLOT(LXFILE,LXFMT,LXPROJ,LXUNIT,LXTYPE)
      if(txfile.ne.blankt) call txplot
      else
      xxscal=1.
      yyscal=1.
      endif
c
c  rescale plotter for grid.
c
      xx(1)=xx(1)*xxscal
      xx(2)=xx(2)*xxscal
      yy(1)=yy(1)*yyscal
      yy(2)=yy(2)*yyscal
      call scale(xx,yy,xp1,yp1,3,ier)
  140 return
      end
C
C
C________________________________________________________________
C
C     SUBROUTINE  H I L O W
C________________________________________________________________
C
        subroutine hilow(jx,zmax,jm,zmin,w,jref,nc2,iw,sizel,
     1                   DGRAD,NOHI)
c
c  routine to plot h or l symbols at local maxima/minima on
c  contours.
c  'local' is defined by a search radius in grid units,
c  parameter 'iw'
c  is the window size ... 2*radius+1.
c  coded by Mike Webring 1/80.
c
        common/labcom/lchars(6),fmtc,nchar,dum(5)
        real w(nc2,iw),jx(nc2),jm(nc2),zmax(nc2),zmin(nc2),
     1            jref(iw)
        character id*56,p*8,fmtc*16,chars*24
        equivalence (lchars,chars)
c  dv is the maximum floating point number on the VAX
        data low/76/,ihigh/72/,dv/1.e38/,fmax/1.e37/
        read(10) id,p,nc,nr,nz,xo,dx,yo,dy
        if(nc.lt.iw .or. nz.gt.1 .or. dx.eq.0. .or.
     &  dy.eq.0.) return
c
        iht=0
        ilt=0
        ihw=iw/2
        ihw1=ihw+1
C
C RE-EVALUATE DGRAD SO THAT THE VALUE PASSED TO THIS SUBROUTINE
C IS THE DIFFERENCE IN Z BETWEEN THE CENTER OF THE WINDOW AND
C THE EDGE OF THE WINDOW.
C
      DGRAD=2.0*DGRAD/(FLOAT(IHW)*(DX+DY)/2.0)**2
C
        jwt=0
        midl=ihw1
        jwl=iw+1
        lmtx=nc2-ihw+1
        do 4 j=1,iw
        do 1 i=1,nc2
1       w(i,j)=dv
4       jref(j)=j
        ir=0
        do 2 j=ihw1,iw
        call rowio(nc,w(ihw1,j),-1,10,10,ie)
2       ir=ir+1
        do 3 i=1,nc2
        zmax(i)=-1.e20
        zmin(i)=1.e20
        jx(i)=0
3       jm(i)=0
        iptr=1
        inew=1
        istop=0
C        IY=0
        y=yo
c
100     do 10 i=ihw1,nc2-ihw
        if(jx(i).gt.jwt) go to 12
        zmax(i)=-1.e20
        do 11 j=1,iw
        w2=w(i,j)
        if(w2.lt.zmax(i)) go to 11
        if(w2.gt.fmax) go to 11
        zmax(i)=w2
        jx(i)=jref(j)
11      continue
        go to 10
12      w2=w(i,inew)
        if(w2.lt.zmax(i)) go to 10
        if(w2.gt.fmax) go to 10
        zmax(i)=w2
        jx(i)=jref(inew)
10      continue
c
        do 15 i=ihw1,nc2-ihw
        if(jm(i).gt.jwt) go to 17
        zmin(i)=1.e20
        do 16 j=1,iw
        w2=w(i,j)
        if(w2.gt.zmin(i)) go to 16
        zmin(i)=w2
        jm(i)=jref(j)
16      continue
17      w2=w(i,inew)
        if(w2.gt.zmin(i)) go to 15
        zmin(i)=w2
        jm(i)=jref(inew)
15      continue
c
        m=ihw1
101     iwl=m+ihw
        iwt=m-ihw
        if(jx(m).ne.midl) go to 18
        do 19 i=iwt,iwl
        if(zmax(i).gt.zmax(m)) go to 22
19      continue
        x=xo+float(m-ihw1)*dx
          CALL LHDGRD(W,NC2,IW,M,IHW1,DX,DY,DGRADN,DGRADX,
     &    dirlh,iptr)
        IF(DGRAD.GT.0.) THEN
          IF(ABS(DGRADX).LT.DGRAD.OR.DGRADX.GE.FMAX) GOTO 32
        END IF
        wq=-dv
        do 31 jq=1,iw
        if(w(m,jq).gt.wq.and.w(m,jq).lt.fmax) wq=w(m,jq)
   31   continue
        if(nohi.eq.0) then
        call vchar(x,y,ihigh,1,1,sizel,0.,0.,0.)
        else if(wq.lt.fmax) then
        call fmtmod(wq,nchar,fmtc)
        write(chars,fmtc) wq
        if(nchar.gt.0) call vchar(x,y,lchars,nchar,2,sizel,
     &  dirlh,-float(nchar/2)*sizel,0.)
        endif
32      m=m+ihw1
        go to 20
18      if(jm(m).ne.midl) go to 22
        do 21 i=iwt,iwl
        if(zmin(i).lt.zmin(m)) go to 22
21      continue
        x=xo+float(m-ihw1)*dx
          CALL LHDGRD(W,NC2,IW,M,IHW1,DX,DY,DGRADN,DGRADX,
     &    dirlh,iptr)
        IF(DGRAD.GT.0.) THEN
          IF(ABS(DGRADN).LT.DGRAD.OR.DGRADN.GE.FMAX) GOTO 34
        END IF
        wq=dv
        do 33 jq=1,iw
        if(w(m,jq).lt.wq.and.w(m,jq).lt.fmax) wq=w(m,jq)
   33   continue
        if(nohi.eq.0) then
        call vchar(x,y,low,1,1,sizel,0.,0.,0.)
        else if(wq.lt.fmax) then
        call fmtmod(wq,nchar,fmtc)
        write(chars,fmtc) wq
        if(nchar.gt.0) call vchar(x,y,lchars,nchar,2,sizel,
     &  dirlh,-float(nchar/2)*sizel,0.)
        endif
34      m=m+ihw1
        go to 20
c
22      m=m+1
20      if(m.lt.lmtx) go to 101
c
        if(ir.eq.nr) go to 98
        call rowio(nc,w(ihw1,iptr),-1,10,10,ie)
        if(ie.ne.0) go to 99
        ir=ir+1
        go to 96
98      if(istop.eq.ihw) go to 99
        do 97 i=1,nc2
97      w(i,iptr)=dv
        istop=istop+1
c
96      inew=iptr
        jref(iptr)=jwl
        iptr=iptr+1
        if(iptr.gt.iw) iptr=1
        jwt=jwt+1
        midl=midl+1
        jwl=jwl+1
C        IY=IY+1
C        Y=DY*IY+Y0
        y=y+dy
        go to 100
99      return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  R O W I O
C________________________________________________________________
C
        subroutine rowio(n,z,iop,idev,jdev,iend)
        dimension z(n)
        iend=0
        if(iop)1,2,1
1       read(idev,end=10) y,z
        if(iop)9,9,2
2       write(jdev) y,z
9       return
10      iend=1
        return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  V E C T O R
C________________________________________________________________
C
        subroutine vector(xx,yy,sizel)
c
c       plots random data values and positions
c       in either vector or scalar form
c
        common /vect/ ispost,lowhi,nid,ncharv,szpost,vmin,vmax,
     1    lifmtv(14),lfmtv(4),ich(20),lchid(20),DGRAD,vdum1(2)
        dimension xx(2),yy(2),bz(6)
        character*4 chid(20)
        character*56 ifmtv,blank
        character ivl*3,sta*24,fmtv*16,bid*8
        logical free
        integer ista(6)
        equivalence(ivl,jvl),(sta,ista(1))
        equivalence (lifmtv,ifmtv),(lfmtv,fmtv),(lchid,chid)
        data con/.017453292/,pi/3.14159265/,a/.9423/,
     1    blank/' '/
        free=.false.
        if(ifmtv(1:1).eq.'*') free=.true.
c
        if(nid.gt.19) nid=19
        if(ncharv.gt.24) ncharv=24
        jch=ich(1)
        if(ispost.gt.1 .or. ispost.lt.0) go to 100
        dmax=1.e-10
1       if(ifmtv.ne.blank)
     1    read(11,ifmtv,end=3)x,y,vm,vi,vd
        if(ifmtv.eq.blank)
     1    read(11,end=3) x,y,vm,vi,vd
        if(vm.gt.dmax) dmax=vm
        go to 1
3       continue
        sca=(vmax-vmin)/alog(dmax)
        pima=pi-a
        pia=pi+a
        pima2=2.*pi-a
        rewind(11)
10       if(ifmtv.ne.blank)
     1    read(11,ifmtv,end=99)x,y,vm,vi,vd
        if(ifmtv.eq.blank)
     1    read(11,end=99) x,y,vm,vi,vd
        incl=int(vi)
        v=alog(vm)
        if(v.lt.0.) v=0.
        vlen=sca*v+vmin
        th=-vd*con
        if((th.gt.a .and. th.lt.pima)  .OR.
     1    (th.gt.pia .and. th.lt.pima2)) THEN
         XL = SIZEL
         YL = 0.0
        ELSE IF((TH.le.a)  .OR.
     1    (th.ge.pi .and. th.le.PIA)) THEN
         XL = 0.0
         YL = -SIZEL
         ELSE
         XL = 0.0
         YL = SIZEL
        ENDIF
c        WRITE (IVL,FMTV) INCL
        WRITE (IVL,'(i3)') INCL
c        CALL VCHAR(x,y,3,1,1,.05,0.,0.,0.)
        if(szpost.gt.0.)
     1      CALL VCHAR(x,y,3,1,1,szpost,0.,0.,0.)
c        CALL VCHAR(x(i),y(i),arrow,1,1,vlen,th,0.,0.)
        CALL VCHAR(x,y,94,1,1,vlen,th,0.,0.)
        if(ncharv.ne.0)
     1      CALL VCHAR(x,y,jvl,3,2,sizel,0.,xl,yl)
c10      continue
        go to 10
c        go to 99
c
100     m1=1
        m2=2
        if(xx(1).lt.xx(2)) go to 140
        m1=2
        m2=1
140     m3=1
        m4=2
        if(yy(1).lt.yy(2)) go to 141
        m3=2
        m4=1
141     yoff=-sizel*.5
        xoff=2.*sizel
c
c     put in by r.godson to plot different intervals
c     of number of input points
c
      ic=0
      if(ispost.gt.-31.and.ispost.lt.3) then
      icnt=1
      else if(ispost.gt.2) then
      icnt=ispost/2
      else
      icnt=iabs(ispost)/31
      ispost=ispost-(-31*icnt)
      endif
        iz=mod(iabs(ispost),10)
        idbr=iabs(ispost)/10 - 1
         ibr=-1
         if(ifmtv.eq.blank) ibr=0
         if(ispost.lt.0) ibr=1
        if(ibr) 101,102,103
101     if(free) then
        read(11,*,end=99) sx,sy,sz
        else
        read(11,ifmtv,end=99) sx,sy,sz
        endif
        go to 104
102     read(11,end=99) sx,sy,sz
        go to 104
103     read(11,end=99) bid,sx,sy,bz
        sz=bz(iz)
104     if(sx.lt.xx(m1) .or. sx.gt.xx(m2)) go to 122
        if(sy.lt.yy(m3) .or. sy.gt.yy(m4)) go to 122
        if(nid.ne.0) then
         call match(nid,chid,bid,isym)
         jch=ich(isym)
        endif
      if(mod(ic,icnt).eq.0) then
      call vchar(sx,sy,jch,1,1,szpost,0.,0.,0.)
      endif
      ic=ic+1
        if(idbr.ne.0 .and. ncharv.gt.0) then
         write(sta,fmtv)sz
         n=leftj(sta)
         call vchar(sx,sy,ista,n,2,sizel,0.,xoff,yoff)
        endif
        if(idbr.ge.0) then
      if(bid.ne.'       ') then
         sta=bid
         n=leftj(sta)
         call vchar(sx,sy,ista,n,2,sizel,.785,xoff,0.)
      endif
        endif
122     if(ibr) 101,102,103
99      return
        end
C
C________________________________________________________________
C
C     FUNCTION  L E F T J
C________________________________________________________________
C
        function leftj(a)
c  left justifies a string and returns the position
c  of the last nonblank character
        character a*(*)
        n=len(a)
        if(a(1:1).ne.' ') go to 15
        do 1 m=2,n
1       if(a(m:m).ne.' ') go to 5
        leftj=0
        return
5       i2=1
        do 10 i=m,n
        a(i2:i2)=a(i:i)
        a(i:i)=' '
10      i2=i2+1
        n=n-m+1
15      do 20 leftj=n,1,-1
20      if(a(leftj:leftj).ne.' ') go to 25
25      return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  M A T C H
C________________________________________________________________
C
        subroutine match(n,chk,test,isym)
c       compares two character strings and returns
c       a symbol index number
        character*8 test
        character*4 chk(20)
        do 10 isym = 1,n
        if(chk(isym).eq.' ') go to 10
        length = index(chk(isym)//' ', ' ')-1
        if (index(test,chk(isym)(1:length)) .gt. 0) return
10      continue
        isym=n+1
        return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  L L P O S T
C________________________________________________________________
C
        subroutine llpost
        common /llp/ latm(3),latx(3),longm(3),longx(3),cm(3),
     1               baslat(3),
     1    iproj,xxx(2),yyy(2),sizep,unit,iplotr,neat,tint,itpost,
     1    ib,sizet,phi1,phi2
        dimension scaf(3),range(2),xx(2),yy(2)
        real latm,latx,longm,longx
        character fmt*16
        logical nos,nom
        data d2r/1.7453292e-2/, conm/1.666667e-2/,
     1       cons/2.777778e-4/,
     1    scaf/39.370079,1.,.001/, icross/6/,
     1    nmax/10/
        dms2d(a,b,c)=a+sign(b,a)*conm+sign(c,a)*cons
c
        if(iproj.lt.1 .or.iproj.gt.10) return
        sca=scaf(int(unit)+1)
        do 1 i=1,2
        xx(i)=xxx(i)/(1000.*sca)
1       yy(i)=yyy(i)/(1000.*sca)
c       changes made by r.godson 3/90 to substitute xxx & yyy
c       for xx & yy in several places below to allow proper
c       labelling when units are not in kilometers.
        hfx=.5*(xxx(1)+xxx(2))
        hfy=.5*(yyy(1)+yyy(2))
        if(iproj.gt.4) call setalb(iproj)
        bld=dms2d(baslat(1),baslat(2),baslat(3))
        cmd=dms2d(cm(1),cm(2),cm(3))
        if(abs(bld).gt.90. .or. abs(cmd).gt.360.) return
        call prjctl(bld,cmd,xb,yb,cmd,sca,iproj)
        ydm=dms2d(latm(1),latm(2),latm(3))
        ydx=dms2d(latx(1),latx(2),latx(3))
        xdx=dms2d(longm(1),longm(2),longm(3))
        xdm=dms2d(longx(1),longx(2),longx(3))
c  adjust cem for conic projection distance to equator
        cem=4.0e7
        if(iproj.eq.4 .and. abs(bld).lt.15.) cem=4.21e7
        if(iproj.gt.4 .and. abs(bld).lt.15.) cem=3.84e7
        ce=cem*sca
        dt=(yy(2)-yy(1))*360./(ce*5.)
c  min-max latitude of area
        if(ydm.ge.ydx) go to 9
        if(ydm.ne.0.0 .and. ydx.ne.0.0) go to 10
9       ydm=bld+(yy(1)*360.)/ce-dt
        ydx=bld+(yy(2)*360.)/ce+dt
        if(ydm.lt.-90.0) ydm=-90.0
        if(ydx.gt.90.0)  ydx=90.0
        latx(1)=ydx
        latm(1)=ydm
c  min-max longitude of area
10      if(xdm.ne.0.0 .and. xdx.ne.0.0) go to 20
        cosp=cos(d2r*amax1(abs(ydx),abs(ydm)))
        if(cosp.lt..17364) cosp=.17364
        xdm=cmd+(xx(1)*360.)/(cosp*ce)-dt
        xdx=cmd+(xx(2)*360.)/(cosp*ce)+dt
        longm(1)=xdx
        longx(1)=xdm
c  default tick interval for less than nmax ticks per side
c  result is multiple of 1, 2 or 5 deg, min, or, sec
20      if(tint.gt.0.0) go to 30
        range(1)=0.
        range(2)=amax1((xdx-xdm),(ydx-ydm))
        tint=0.
        dr=range(2)*60./float(nmax)
        if(dr.ge.60.) then
         call setax(range,tint,nmax,ndum,fmt)
         tint=tint*60.
        else if(dr.ge.1. .and. dr.lt.60.) then
         range(2)=range(2)*60.
         call setax(range,tint,nmax,ndum,fmt)
        else
         range(2)=range(2)*3600.
         call setax(range,tint,nmax,ndum,fmt)
         tint=tint/60.
        endif
c  final initialization
30      if(itpost.le.0) itpost=2
        dt=(aint(tint*60.+.0001))/3600.
        dt2=dt*float(itpost)
        dt3=dt2*60.0
        nos=.false.
        if(dt3-aint(dt3).eq.0.0) nos=.true.
        nom=.false.
        if(amod(dt3,60.0).eq.0.0) nom=.true.
        if(yy(2)-yy(1).lt.20.)  ydm=aint(60.*ydm)/60.
        if(yy(2)-yy(1).ge.20.)  ydm=aint(ydm+90.)-90.
        if(yy(2)-yy(1).gt.500.) ydm=5.*aint((ydm+90.)/5.) - 90.
        if(xx(2)-xx(1).lt.20.)  xdm=aint(60.*xdm)/60.
        if(xx(2)-xx(1).ge.20.)  xdm=aint(xdm+360.)-360.
c     utm has a limit of 3 degrees
        if(xx(2)-xx(1).gt.500. .and. iproj.ne.2)
     1    xdm=5.*aint((xdm+360.)/5.) - 360.
c  plot longitude labels
        if(sizep.le.0.0) return
        xd=xdm
100     if(xd.gt.xdx) go to 190
        call dms(xd,id,im,is)
        id=id
        yd=ydm
105     call prjctl(yd,xd,x,y,cmd,sca,iproj)
        y=y-yb
        call prjctl(yd-dt,xd,xm,ym,cmd,sca,iproj)
        ym=ym-yb
        xm=x-(x-xm)/(y-ym)*(y-yyy(1))
        if(x.ge.xxx(1).and.x.le.xxx(2) .and.
     1    y.ge.yyy(1).and.y.lt.hfy    .and.
     2    abs(y-yyy(1)).le.abs(y-ym)) go to 109
        yd=yd+dt
        if(yd-ydx) 105,110,110
109     call latlab(-1.,0,sizep,xm,yyy(1),id,im,is,nos,nom)
110     if(iplotr.eq.1 .or. iplotr.eq.4) go to 120
        yd=ydx
115     call prjctl(yd,xd,x,y,cmd,sca,iproj)
        y=y-yb
        call prjctl(yd+dt,xd,xp,yp,cmd,sca,iproj)
        yp=yp-yb
        xp=x+(xp-x)/(yp-y)*(yyy(2)-y)
        if(x.ge.xxx(1).and.x.le.xxx(2) .and.
     1    y.gt.hfy.and.y.le.yyy(2)    .and.
     2    abs(y-yyy(2)).le.abs(y-yp)) go to 119
        yd=yd-dt
        if(yd-ydm) 120,115,115
119     call latlab(1.,0,sizep,xp,yyy(2),id,im,is,nos,nom)
120     xd=xd+dt2
        go to 100
c  plot latitude labels
190     yd=ydm
200     if(yd.gt.ydx) go to 300
        call dms(yd,id,im,is)
        xd=xdm
205     call prjctl(yd,xd,x,y,cmd,sca,iproj)
        call prjctl(yd,xd-dt,xm,ym,cmd,sca,iproj)
        y=y+(ym-y)/(x-xm)*(x-xxx(1))
        y=y-yb
        if(y.ge.yyy(1).and.y.le.yyy(2) .and.
     1    x.ge.xxx(1).and.x.lt.hfx     .and.
     2    abs(x-xxx(1)).le.abs(x-xm)) go to 209
        xd=xd+dt
        if(xd-xdx) 205,210,210
209     call latlab(-1.,1,sizep,xxx(1),y,id,im,is,nos,nom)
210     if(iplotr.eq.1 .or. iplotr.eq.4) go to 220
        xd=xdx
215     call prjctl(yd,xd,x,y,cmd,sca,iproj)
        call prjctl(yd,xd+dt,xp,yp,cmd,sca,iproj)
        y=y+(yp-y)/(xp-x)*(xxx(2)-x)
        y=y-yb
        if(y.ge.yyy(1).and.y.le.yyy(2) .and.
     1    x.gt.hfx.and.x.le.xxx(2)     .and.
     2    abs(x-xxx(2)).le.abs(x-xp)) go to 219
        xd=xd-dt
        if(xd-xdm) 220,215,215
219     call latlab(1.,1,sizep,xxx(2),y,id,im,is,nos,nom)
220     yd=yd+dt2
        go to 200
c  plot interior lat-lon ticks
300     yd=ydm
        if(sizet.le.0.0) return
        do 360 j=1,int((ydx-yd)/dt)+1
        xd=xdm
        do 350 i=1,int((xdx-xd)/dt)+1
        call prjctl(yd,xd,x,y,cmd,sca,iproj)
        y=y-yb
        if(x.lt.xxx(1) .or. x.gt.xxx(2)) go to 350
        if(y.lt.yyy(1) .or. y.gt.yyy(2)) go to 350
        call prjctl(yd+dt,xd,xp,yp,cmd,sca,iproj)
        call prjctl(yd-dt,xd,xm,ym,cmd,sca,iproj)
        xp=xm-xp
        yp=yp-ym
        ang=atan2(xp,yp)
c        if(sizep.lt..00025) then
c          call vchar(x,y,icross,1,2,1000.0*sizep,ang,0.,0.)
c        else
          call vchar(x,y,icross,1,2,sizet,ang,0.,0.)
c        end if
350     xd=xd+dt
360     yd=yd+dt
        return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  D M S
C________________________________________________________________
C
        subroutine dms(d,id,im,is)
        is=int(abs(d)*3600.+.5)
        id=is/3600
        im=(is-id*3600)/60
        is=is-(id*3600+im*60)
        id=isign(id,int(d))
        return
        end
C
C________________________________________________________________
C
C     SUBROUTINE  L A T L A B
C________________________________________________________________
C
        subroutine latlab(side,lat,sz,x,y,jd,jm,js,nosec,nomin)
        dimension labint(3)
        character label*12,jd2*1,jm2*1,js2*1
        equivalence (labint,label)
        logical nosec,nomin
        jd2=char(24)
        jm2=char(39)
        js2=char(34)
        write(label,50) iabs(jd),jd2,jm,jm2,js,js2
50      format(i4,a1,2(i2,a1))
        if(jm.ne.0 .or. js.ne.0) nomin=.false.
        if(js.ne.0) nosec=.false.
        nch=11
        if(nosec) nch=8
        if(nomin) nch=5
        xl=nch*sz
        yof=0.
        if(lat.eq.1) then
        if(side.lt.0.) xof=-(2.*sz+xl)
        if(side.gt.0.) xof=2.*sz
        else
        xof=-xl/2.
        yof=2.*sz*side
        endif
        call vchar(x,y,labint,nch,2,sz,0.,xof,yof)
        return
        end
C
C________________________________________________________________
C
C   SUBROUTINE  L H D G R D
C________________________________________________________________
C
C SUBROUTINE LOWHI DGRADIENT WHEN CALLED BY CONTOUR RETURNS THE
C MINIMUM AND MAXIMUM CURVATURE (DELTA GRADIENT) FOUND IN THE
C WINDOW DETERMINED BY SUBROUTINE HILOW.  THE PROCEDURE IS TO
C FIT A LEAST SQUARES PARABOLA TO EACH OF 4 SLICES TAKEN THROUGH
C THE CENTER OF THE WINDOW AND THEN FIND THE MINIMUM AND MAXIMUM
C OF TWICE THE SECOND ORDER COEFFICIENT.
C
C SUBROUTINE LHDGRD WRITTEN BY ROB BRACKEN, USGS, 6MAR87
C
C
      subroutine lhdgrd(w,nc2,iw,ic,jc,dx,dy,dgradn,dgradx,
     & dirlh,iptr)
C
C VARIABLE DECLARATION
C
      real*4 w(nc2,iw),dgrad(4),dgdir(4)
      real*8 x(1000),y(1000)
      real*8 a,b,c
      data ddval/1.e38/,dval/1.e+37/
      data dgdir/0.,-1.570796,.7853982,-.7853982/
C
C SET UP HALF WINDOW WIDTHS
C
      ihw=iw/2
      ihw1=ihw+1
C
C SET UP LEFT, BOTTOM, AND TOP COORDINATES
C
      il=ic-ihw1
      jb=jc-ihw1
      jt=jc+ihw1
C
C FIND CURVATURE OF E-W SLICE
C
       k=iptr+ihw
       if(k.gt.iw) k=k-iw
       do 10 i=1,iw
        x(i)=dble(dx*i)
        y(i)=dble(w(il+i,k))
   10 continue 
      call lsqpbl(x,y,1,iw,a,b,c,ei)
      dgrad(1)=2.0*a
      if(ei.ne.0) dgrad(1)=ddval
C
C FIND CURVATURE OF S-N SLICE
C
       do 20 j=1,iw
        x(j)=dble(dy*j)
        k=iptr-1+j
        if(k.gt.iw) k=k-iw
        y(j)=dble(w(ic,k))
   20 continue
      call lsqpbl(x,y,1,iw,a,b,c,ei)
      dgrad(2)=2.0*a
      if(ei.ne.0) dgrad(2)=ddval
C
C FIND CURVATURE OF SW-NE SLICE
C
      diag=sqrt(dx*dx+dy*dy)
       do 30 ij=1,iw
        x(ij)=dble(diag*ij)
        k=iptr-1+ij
        if(k.gt.iw) k=k-iw
        y(ij)=dble(w(il+ij,k))
   30 continue
      call lsqpbl(x,y,1,iw,a,b,c,ei)
      dgrad(3)=2.0*a
      if(ei.ne.0) dgrad(3)=ddval
C
C FIND CURVATURE OF NW-SE SLICE
C
       do 40 ji=1,iw
        k=iptr-ji
        if(k.lt.1) k=k+iw
        y(ji)=dble(w(il+ji,k))
   40 continue
      call lsqpbl(x,y,1,iw,a,b,c,ei)
      dgrad(4)=2.0*a
      if(ei.ne.0) dgrad(4)=ddval
C
C FIND MINIMUM AND MAXIMUM CURVATURE
C
      dgradn=ddval
      dgradx=-1*ddval
      imx=1
      imn=1
       do 50 i=1,4
         if(dgrad(i).lt.dval) then
          if(dgrad(i).gt.dgradx) then
            dgradx=dgrad(i)
            imx=i
          endif
          if(dgrad(i).lt.dgradn) then
            dgradn=dgrad(i)
            imn=i
          endif
         end if
   50 continue
      if(dgradx.lt.-1*dval) dgradx=ddval
      if(dgradx.lt.0.) dirlh=dgdir(imx)
      if(dgradn.gt.0.) dirlh=dgdir(imn)
      return
      end
C
C________________________________________________________________
C
C   SUBROUTINE  L S Q P B L
C________________________________________________________________
C
C LEAST SQUARE PARABOLA  --  FINDS THE PARABOLA WHICH
C IS THE BEST FIT TO THE LOCUS OF POINTS PASSED TO
C THE SUBROUTINE.
C
C FOLLOWING IS A LIST OF THE PASSED VARIABLES:
C
C    X     - DOUBLE PRECISION ARRAY OF INDEPENTANT VARIABLES
C    Y     - DOUBLE PRECISION ARRAY OF DEPENDANT VARIABLES
C    IB    - INTEGER BEGINNING POSITION IN ARRAYS X AND Y FOR FIT
C    IE    - INTEGER ENDING POSITION IN ARRAYS X AND Y FOR FIT
C    A,B,C - DOUBLE PRECISION COEFFICIENTS FOR THE BEST FIT
c            PARABOLA ACCORDING TO THE EQUATION  y = Axx + Bx + C
C    EI    - REAL ERROR INDICATOR
C                0.0 - SUCCESSFUL CALCULATION
C                1.0 - DIVISION BY ZERO ERROR
C
C SUBROUTINE LSQPBL WRITTEN BY  ROB BRACKEN, USGS, 2MAY86
C
      subroutine lsqpbl( x,y, ib,ie,   a,b,c, ei )
C
C VARIABLE DECLARATION
C
      double precision x(1),y(1)
      double precision a,b,c
      double precision x1,x2,x3,x4,yx1,yx2,y1
      double precision dd,d10,d20,d11,d21,d12,d22,xs
      data dval/1.e+37/
C
C VARIABLE INITIALIZATION
C
      ei=1.0
      x1=0.d0
      x2=0.d0
      x3=0.d0
      x4=0.d0
      yx1=0.d0
      yx2=0.d0
      y1=0.d0
C
C FIND SUMS
C
      ndval=0
      do 10 i=ib,ie
        if(x(i).lt.dval.and.y(i).lt.dval) then
          x1=x1+x(i)
          xs=x(i)*x(i)
          x2=x2+xs
          x3=x3+xs*x(i)
          x4=x4+xs*xs
          yx1=yx1+y(i)*x(i)
          yx2=yx2+y(i)*xs
          y1=y1+y(i)
        else
          ndval=ndval+1
        end if
   10 continue
      n=ie-ib+1-ndval
C
C CHECK FOR DIVISION BY ZERO
C
      if(dabs(x2).lt.1.d-30) return
      if(dabs(x3).lt.1.d-30) return
      if(dabs(x4).lt.1.d-30) return
      if(n.lt.3) return
C
C FIND ELEMENTS OF 2 X 2 DETERMINANT
C
      d10=yx2/x4-yx1/x3
      d20=yx1/x3-y1/x2
      d11=x3/x4-x2/x3
      d21=x2/x3-x1/x2
      d12=x2/x4-x1/x3
      d22=x1/x3-dble(float(n))/x2
      dd=d11*d22-d12*d21
C
C FIND COEFFICIENTS
C
      if(dabs(dd).lt.1.d-30) return
      c=(d11*d20-d10*d21)/dd
      b=(d10*d22-d12*d20)/dd
      a=(y1-b*x1-c*n)/x2
C
C CALCULATION COMPLETED
C
      ei=0.0
      return
      end
C
C_______________________________________________________________
C
C     SUBROUTINE  B S P L I N E
C_______________________________________________________________
C
c      subroutine bsplin(x,y,n,xx,yy,nn,step,sigmap)
      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
c      s1=1.-step
      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
