      SUBROUTINE LTRP (M,X,Y,N,XI,YI,NI,T,IERR)
C     ******************************************************************
C     LAGRANGE INTERPOLATION
C     ******************************************************************
      DIMENSION X(N),Y(N),XI(NI),YI(NI),T(M)
C
C     CHECK INPUT
C
      IF (M.LT.2) GO TO 130
      IF (M.GT.N) GO TO 131
      IF (NI.LT.1) GO TO 132
      IERR = 0
C
C     INITIALIZATION
C
      MM1 = M-1
      K  = 1
      XX = XI(1)
      ILOLD = 0
C
C     FIND THE SUBINTERVAL WHICH CONTAINS XX. I = 1 IF XX.LT.X(2)
C     AND I = N IF XX.GE.X(N). OTHERWISE X(I).LE.XX.LT.X(I+1).
C
      IF (XX-X(1)) 10,11,20
   10 I = 1
      IL = 1
      IR = M
      GO TO 80
   11 I = 1
      YI(K) = Y(1)
      GO TO 110
C
   20 IF (X(N)-XX) 21,22,23
   21 I = N
      IL = N-M+1
      IR = N
      GO TO 80
   22 I = N
      YI(K) = Y(N)
      GO TO 110
   23 IL = 1
      IR = N
C
C     BISECTION SEARCH
C
   30 I = (IL+IR)/2
      IF (I.EQ.IL) GO TO 50
      IF (XX-X(I)) 31,32,33
   31 IR = I
      GO TO 30
   32 YI(K) = Y(I)
      GO TO 110
   33 IL = I
      GO TO 30
C
C     LINEAR FORWARD SEARCH
C
   40 IF (XX-X(I+1)) 50,41,42
   41 I = I+1
      YI(K) = Y(I)
      GO TO 110
   42 I = I+1
      GO TO 40
C
C     POINT XX LIES IN THE OPEN SUBINTERVAL (X(I),X(I+1)).
C     FIND THE M CLOSEST POINTS X(IL),...,X(IR) TO XX.
C
   50 IF (I.GT.M) GO TO 51
      IL = 1
      NUM = I
      GO TO 60
   51 IL = I-M+1
      NUM = M
C
   60 IPM = I+M
      IF (IPM.LE.N) GO TO 61
      IR = N
      NUM = NUM+N-I
      GO TO 70
   61 IR = IPM
      NUM = NUM+M
C
   70 NUM = NUM-M
      IF (NUM.EQ.0) GO TO 80
      DL = XX-X(IL)
      DR = X(IR)-XX
      DO 72 L=1,NUM
      IF (DL.LE.DR) GO TO 71
      IL = IL+1
      DL = XX-X(IL)
      GO TO 72
   71 IR = IR-1
      DR = X(IR)-XX
   72 CONTINUE
C
C     COMPUTE THE COEFFICIENTS T(1),...,T(M) OF THE BACKWARD
C     NEWTON FORM OF THE INTERPOLATING POLYNOMIAL.
C
   80 IF (IL.EQ.ILOLD) GO TO 100
      ILOLD = IL
      ILM1 = IL-1
      DO 81 J=1,M
      L = ILM1+J
   81 T(J) = Y(L)
C
      DO 91 ISTEP=1,MM1
      JMAX = M-ISTEP
      DO 90 J=1,JMAX
      II = ILM1+J
      L = II+ISTEP
   90 T(J) = (T(J)-T(J+1))/(X(II)-X(L))
   91 CONTINUE
C
C     EVALUATION OF THE INTERPOLATING POLYNOMIAL
C
  100 YI(K) = T(1)
      DO 101 J=2,M
      L = ILM1+J
  101 YI(K) = T(J)+YI(K)*(XX-X(L))
C
C     NEXT POINT
C
  110 IF (K.GE.NI) RETURN
      K = K+1
      XX = XI(K)
      IF (XX-X(1)) 10,11,111
  111 IF (X(N)-XX) 21,22,120
C
  120 IF (XX-XI(K-1)) 121,122,40
  121 IL = 1
      IR = MIN0(I+1,N)
      GO TO 30
  122 YI(K) = YI(K-1)
      GO TO 110
C
C     ERROR RETURN
C
  130 IERR = 1
      RETURN
  131 IERR = 2
      RETURN
  132 IERR = 3
      RETURN
      END
