      SUBROUTINE LAGRNG (X,ARG,Y,VAL,NDIM,NFS,NPTX,MAXARG,MAXFS)
C *** SEE R H LANDAU, PROGRAM LPOTT, 1981
C     LAGRANGE INTERPOLATION,UNEQUALLY SPACED POINTS
C     NPTS=2,3,4,5,6, NFS FUNCTIONS(Y*S) SIMULTANEOUSLY INTERPLTD
C     X= VALUE OF ARGUMENT, ARG IS THE TABULATED X*S
C     Y= A VECTOR OF INTERPOLAD FUNCTIONS, FROM TABULTED VAL
C     NDIM= DIMENSION OF TABLE
C     NFS= = OF FUNCTIONS SIMULT INTERPOLATED
C     MAXARG AND MAXFS ARE MAXIMUM VALUES OF SUBSCRIPTS,IE DIMENSIONS
      IMPLICIT  REAL*8   (A-H,O-Z)
      DIMENSION ARG(MAXARG), VAL(MAXFS,MAXARG), Y(MAXFS)
C     -----FIND X0 THE CLOSEST POINT TO X
      NPTS = NPTX
C                                                                   2=ST
      NI = 1
      NF = NDIM
 10   IF ((X.LE.ARG(NI)).OR.(X.GE.ARG(NF))) GO TO 30
      IF ((NF-NI+1).EQ.2) GO TO 70
      NMID = (NF+NI)/2
      IF (X.GT.ARG(NMID)) GO TO 20
      NF = NMID
      GO TO 10
 20   NI = NMID
      GO TO 10
C     ------ X IS ONE OF THE TABLTD VALUES
 30   IF (X.LE.ARG(NI)) GO TO 60
      NN = NF
 40   NUSED = 0
      DO 50 N=1,NFS
         Y(N) = VAL(N,NN)
 50   CONTINUE
      RETURN
 60   NN = NI
      GO TO 40
C     ------- 2 PTS LEFT, CHPOSE SMALLER ONE
 70   N0 = NI
      NN = NPTS-1
      GO TO (140,110,100,90,80), NN
 80   CONTINUE
      IF (((N0+3).GT.NDIM).OR.((N0-2).LT.1)) GO TO 90
      NUSED = 6
      GO TO 150
 90   CONTINUE
      IF ((N0+2).GT.NDIM) GO TO 110
      IF ((N0-2).LT.1) GO TO 100
      NUSED = 5
      GO TO 150
 100  CONTINUE
      IF (((N0+2).GT.NDIM).OR.((N0-1).LT.1)) GO TO 110
      NUSED = 4
      GO TO 150
 110  IF ((N0+1).LT.NDIM) GO TO 130
C     N0=NDIM,SPECIAL CASE
 120  NN = NDIM
      GO TO 40
 130  NUSED = 3
      IF ((N0-1).LT.1) NUSED = 2
      GO TO 150
 140  IF ((N0+1).GT.NDIM) GO TO 120
      NUSED = 2
 150  CONTINUE
C     AT LEAST 2 PTS LEFT
      Y0 = X-ARG(N0)
      Y1 = X-ARG(N0+1)
      Y01 = Y1-Y0
      C0 = Y1/Y01
      C1 = -Y0/Y01
      IF (NUSED.EQ.2) GO TO 160
C     AT LEAST 3 PTS
      YM1 = X-ARG(N0-1)
      Y0M1 = YM1-Y0
      YM11 = Y1-YM1
      CM1 = -Y0*Y1/Y0M1/YM11
      C0 = C0*YM1/Y0M1
      C1 = -C1*YM1/YM11
      IF (NUSED.EQ.3) GO TO 180
C     ------AT LEAST 4 PTS
      Y2 = X-ARG(N0+2)
      YM12 = Y2-YM1
      Y02 = Y2-Y0
      Y12 = Y2-Y1
      CM1 = CM1*Y2/YM12
      C0 = C0*Y2/Y02
      C1 = C1*Y2/Y12
      C2 = -YM1*Y0*Y1/YM12/Y02/Y12
      IF (NUSED.EQ.4) GO TO 200
C     AT LEAST 5 PTS
      YM2 = X-ARG(N0-2)
      YM2M1 = YM1-YM2
      YM20 = Y0-YM2
      YM21 = Y1-YM2
      YM22 = Y2-YM2
      CM2 = YM1*Y0*Y1*Y2/YM2M1/YM20/YM21/YM22
      CM1 = -CM1*YM2/YM2M1
      C0 = -C0*YM2/YM20
      C1 = -C1*YM2/YM21
      C2 = -C2*YM2/YM22
      IF (NUSED.EQ.5) GO TO 220
C     AT LEAST 6 PTS
      Y3 = X-ARG(N0+3)
      YM23 = Y3-YM2
      YM13 = Y3-YM1
      Y03 = Y3-Y0
      Y13 = Y3-Y1
      Y23 = Y3-Y2
      CM2 = CM2*Y3/YM23
      CM1 = CM1*Y3/YM13
      C0 = C0*Y3/Y03
      C1 = C1*Y3/Y13
      C2 = C2*Y3/Y23
      C3 = YM2*YM1*Y0*Y1*Y2/YM23/YM13/Y03/Y13/Y23
      GO TO 240
 160  CONTINUE
      DO 170 N=1,NFS
         Y(N) = C0*VAL(N,N0)+C1*VAL(N,N0+1)
 170  CONTINUE
      GO TO 260
 180  CONTINUE
      DO 190 N=1,NFS
         Y(N) = CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)
 190  CONTINUE
      GO TO 260
 200  CONTINUE
      DO 210 N=1,NFS
         Y(N) = CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C2*VAL(N,N0+
     1   2)
 210  CONTINUE
      GO TO 260
 220  CONTINUE
      DO 230 N=1,NFS
         Y(N) = CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0
     1   +1)+C2*VAL(N,N0+2)
 230  CONTINUE
      GO TO 260
 240  CONTINUE
      DO 250 N=1,NFS
         Y(N) = CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0
     1   +1)+C2*VAL(N,N0+2)+C3*VAL(N,N0+3)
 250  CONTINUE
 260  RETURN
      END
