subroutine lagrangeco(val,vec,veclen,first,pdeg,coef,inrange)
c
c purpose
c   this subroutine determines which points of the data vector vec
c   to use for a lagrange interpolation of degree le 3, and calculates
c   the coefficients for the terms that go into the polynomial sum.
c   variable inrange is used to indicate if the value lies outside the
c   range of values found in vec (see below).
c
c calling parameters
c	val	r8	i	value of the variable
c	vec	r8	i	table of variables where the fn. value is known
c	veclen	i4	i	length of vec
c	first	i4	o	index of first element of vec to use for the
c				interpolation (starts at 0)
c	pdeg	i4	o	degree of polynomial to be used
c	coef 	r8	o	table of coefficients to use for terms in
c				the interpolation
c	inrange	i4	o	/ -1 if val < vec(i) for all i
c				|  0 if val is amongst the vec(i)'s
c				\ +1 if val > vec(i) for all i
c
c revision history
c   24.iv.97 adapted from int_rang.f and lagrangeco.pro
c   6.ii.97 added feature to allow vec elements to be either monotone
c           increasing or monotone decreasing.  (eac)
c   vii.96  original written & tested (eac)
c
c---assumes: vec is ordered:  vec(i) .lt. vec(i+1) for all i =or=
c                             vec(i) .gt. vec(i+1) for all i
c
      integer*4  veclen, locinvec, first, pdeg, inrange, i, k
      integer*4  lastel
      real*8 s
      real*8     val, vec(0:veclen-1), coef(0:3)
      real*8 dsgn
      lastel= veclen-1
c     -- s indicates direction of independent variable in vec.
      s = dsgn(vec(veclen-1)-vec(0))
c     -- see if val is in the range of the table.  if it is not, set up 
c     -- for linear extrapolation.
      inrange=0
      if (dsgn(val - vec(lastel)) .eq. s) then
         inrange=1
         first=lastel-1
         pdeg=1
         goto 100
      endif
      if (dsgn(vec(0) - val) .eq. s) then
         inrange=-1
         first=0
         pdeg=1
         goto 100
      endif
c     -- get here if val is in the range of the table...
      pdeg  = 0
      first = 0
      if (vec(first) .eq. val) goto 100
      i = 1
10    continue
      first = i
      if (vec(first) .eq. val) goto 100
      if (dsgn(vec(i) - val) .eq. s) goto 11
      if (i .eq. lastel) goto 11
      i = i+1
      goto 10
11    i = i-1
      first = i
      pdeg=1
      if (first .ne. 0) then
         first = first-1
         pdeg  = pdeg+1
      endif
      if (first+pdeg .ne. lastel) pdeg=pdeg+1
100   continue
      do i = 0, pdeg
         coef(i) = 1.
         do k=0, pdeg
            if (k .ne. i) coef(i)=coef(i)*
     &       (val- vec(first+k))/(vec(first+i)-vec(first+k))
          enddo
      enddo
      return 
      end