c
      subroutine vecterp(fy,f,ny,fylin,flin,nlin,spval,iprnt)
      real fy(ny),f(ny),fylin(nlin),flin(nlin)
c
      if (ny.lt.2) then
        write(6,"('>>> vecterp: ny must be >= 2: ny=',i3)") ny
        return
      endif
c
c Determine if fy array increases or decreases:
c
      inc = 1
      do k=2,ny
        if (fy(k).ne.spval.and.fy(k-1).ne.spval) then
          if (fy(k)-fy(k-1).lt.0.) inc = -1
          goto 100
        endif
      enddo
      if (iprnt.gt.0)
     +  write(6,"('>>> warning vecterp: could not determine ',
     +    'if fy is increasing or decreasing...')")
 100  continue
c
c Use only that part of fy .ne. spval.
c kbot = first non-spval in fy:
c
      kbot = 0
      do k=1,ny
        if (fy(k).ne.spval) then
          kbot = k
          goto 101
        endif 
      enddo
      if (iprnt.gt.0) then
        write(6,"('>>> warning vecterp: fy all spval?')")
        do k=1,nlin
          flin(k) = spval
        enddo
        return
      endif
 101  continue
c
c Assume that all spval's are either at beginning and/or at end
c (no isolated spval's in middle of array)
c nk = number of fy values to use:
c
      nk = 0
      do k=kbot,ny
        if (fy(k).ne.spval) nk = nk+1
      enddo
c
c     subroutine bracket(x,xx,nx,inc,n1,n2,ier)
c Bracket x in xx(nx), returning lower index in n1 and upper index in n2
c
      do k=1,nlin
        call bracket(fylin(k),fy(kbot),nk,inc,k0,k1,ier)
        if (ier.gt.0) then
          if (iprnt.gt.0) 
     +    write(6,"('>>> vecterp: err from bracket: k=',i2,' nlin=',i3,
     +      ' fylin(k)=',
     +      e12.4,' ny=',i2,' inc=',i2,' fy=',/(6e12.4))") 
     +      k,nlin,fylin(k),ny,inc,(fy(i),i=1,ny)
          flin(k) = spval
        else
          if (fy(k1)-fy(k0).gt.1.e-20.or.fy(k1)-fy(k0).lt.-1.e-20) then
            flin(k) = f(k0) + (f(k1)-f(k0))*(fylin(k)-fy(k0))/
     +        (fy(k1)-fy(k0))
          else
            flin(k) = spval
          endif
        endif
      enddo
      return
      end