      FUNCTION TAB(XT,YT,N,K,X)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XT(100),YT(100),X1(6),Y1(6)
C     THIS ROUTINE USES AITKEN S METHOD OF INTERPOLATION.  FOR A GIVEN
C     INDEPENDENT VARIABLE X, AN INTERPOLATION OF ORDER K IS PERFORMED.
C     THE CALLING SEQUENCE IS Y = TAB(XT,YT,N,K,X)   WHERE:
C        XT = ARRAY OF INDEPENDENT VALUES IN A MONOTONIC SEQUENCE.
C        YT = ARRAY OF DEPENDENT VALUES.
C        N = NUMBER OF INDEPENDENT VALUES.
C        K = ORDER OF INTERPOLATION (MUST BE LESS THAN 6).
C        X = INDEPENDENT VALUE FOR WHICH INTERPOLATION IS DESIRED.
      NPTS=K+1
C     CHECK FOR ASCENDING OR DESCENDING SEQUENCE
      IF (XT(1)-XT(N))5,5,1
C     SEQUENCE IS DESCENDING, CHECK IF EXTRAPOLATION IS NEEDED
 1    IF(X-XT(1))2,2,10
 2    IF(X-XT(N))30,3,3
 3    DO 4 I=1,N
      IF(X-XT(I))4,4,50
 4    CONTINUE
C     SEQUENCE IS ASCENDING, CHECK IF EXTRAPOLATION IS NEEDED
 5    IF(X-XT(1))10,10,20
C     SET UP FOR EXTRAPOLATION AT BEGINNING OF TABLE
 10   DO 15 I=1,NPTS
      X1(I)=XT(I)
 15   Y1(I)=YT(I)
      GO TO 100
 20   IF(X-XT(N))40,30,30
C     SET UP FOR EXTRAPOLATION AT END OF TABLE
 30   DO 35 I=1,NPTS
      NN=N-NPTS+I
      X1(I)=XT(NN)
 35   Y1(I)=YT(NN)
      GO TO 100
C     INTERPOLATION NEEDED, FIND THE TWO TABLE VALUES THAT THE
C     INDEPENDENT VARIABLE,X,LIES BETWEEN
 40   DO 45 I=1,N
      IF(X-XT(I))50,45,45
 45   CONTINUE
      I=N
C     CHOOSE K+1 INDEPENDENT TABLE VALUES CLOSEST TO X
 50   IF(MOD(K,2))51,51,60
   51 IF(DABS(X-XT(I-1))-DABS(XT(I)-X)) 52,56,56
 52   IF(I-2-K/2)10,53,53
 53   IF(I-1+K/2-N)54,54,30
 54   DO 55 J=1,NPTS
      NN=I-2-K/2+J
      X1(J)=XT(NN)
 55   Y1(J)=YT(NN)
      GO TO 100
 56   IF(I-1-K/2)10,57,57
 57   IF(I+K/2-N)58,58,30
 58   DO 59 J=1,NPTS
      NN=I-K/2+J-1
      X1(J)=XT(NN)
 59   Y1(J)=YT(NN)
      GO TO 100
 60   IF(I-1+(K+1)/2-N)61,61,30
 61   IF(I-2-(K-1)/2)10,62,62
 62   DO 63 J=1,NPTS
      NN=I-1-(K+1-2*J)/2
      X1(J)=XT(NN)
 63   Y1(J)=YT(NN)
 100  DO 110 II=1,K
      SAVEY1=Y1(1)
      L=II+1
      DO 110 JJ=L,NPTS
      NN=JJ-II
      LL=NN+1
 110  Y1(NN)=(SAVEY1*(X1(JJ)-X)-Y1(LL)*(X1(II)-X))/(X1(JJ)-X1(II))
      TAB=Y1(1)
      RETURN
      END
