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,igrid)
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*1 ctype
        common/captur/icapt,iwhite,ctype
        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(igrid) 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
        if(ctype.eq.'l') then
          icol = 1
        else
          icol = -1
        endif

        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,igrid,igrid,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,icol*1,sizel,0.,0.,0.)
        else if(wq.lt.fmax) then
        call fmod(wq,nchar,fmtc)
        write(chars,fmtc) wq
        if(nchar.gt.0) call vchar(x,y,lchars,nchar,icol*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,icol*1,sizel,0.,0.,0.)
        else if(wq.lt.fmax) then
        call fmod(wq,nchar,fmtc)
        write(chars,fmtc) wq
        if(nchar.gt.0) call vchar(x,y,lchars,nchar,icol*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,igrid,igrid,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  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  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 fmod(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      print *,'lenfmt,lendec,lencut',lenfmt,lendec,lencut
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
c      print *,'fcont2,nleft,lenfmt,nchar',fcont2,nleft,lenfmt,nchar
      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
