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
 
      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

d     write(6,*)'> lagrangeco'

      lastel= veclen-1
c------------------------------------------------------------------
c s indicates direction of independent variable in vec.
c
      s= dsgn(vec(veclen-1)-vec(0))

d     write(6,91)s
d  91 format(/':::s=',f4.0)
d     write(6,92)val, vec(0), dsgn(vec(0) - val)
d  92 format('val=',f6.2,'  vec(0)=',f6.2,'  dsgn=',f4.0)
d     write(6,93)val, vec(lastel), dsgn(val - vec(lastel))
d  93 format('val=',f6.2,'  vec(n)=',f6.2,'  dsgn=',f4.0)

c------------------------------------------------------------------
c See if val is in the range of the table.  If it is not, set up 
c for linear extrapolation.
c      
      inrange=0
      IF (dsgn(val - vec(lastel)) .EQ. s) THEN
        inrange=1
        first=lastel-1
        pdeg=1
        GOTO 100
      END IF

      IF (dsgn(vec(0) - val) .EQ. s) THEN
        inrange=-1
        first=0
        pdeg=1
        GOTO 100
      ENDIF

c-------------------------------------------------------------------
c Get here if val is in the range of the table...
c
      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
      END IF
      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))
        END DO
      END DO

d     write(6,*)'< lagrangeco'
      RETURN 
      END
      FUNCTION dsgn(x)
      REAL*8 dsgn, x

      IF (x) 1,2,3

    1 dsgn= -1.d0
      RETURN
    2 dsgn= 0.d0
      RETURN
    3 dsgn= 1.d0
      RETURN
      END