C
C     ..................................................................
C
C        SUBROUTINE ACFI
C
C        PURPOSE
C           TO INTERPOLATE FUNCTION VALUE Y FOR A GIVEN ARGUMENT VALUE
C           X USING A GIVEN TABLE (ARG,VAL) OF ARGUMENT AND FUNCTION
C           VALUES.
C
C        USAGE
C           CALL ACFI (X,ARG,VAL,Y,NDIM,EPS,IER)
C
C        DESCRIPTION OF PARAMETERS
C           X      - THE ARGUMENT VALUE SPECIFIED BY INPUT.
C           ARG    - THE INPUT VECTOR (DIMENSION NDIM) OF ARGUMENT
C                    VALUES OF THE TABLE (POSSIBLY DESTROYED).
C           VAL    - THE INPUT VECTOR (DIMENSION NDIM) OF FUNCTION
C                    VALUES OF THE TABLE (DESTROYED).
C           Y      - THE RESULTING INTERPOLATED FUNCTION VALUE.
C           NDIM   - AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF
C                    POINTS IN TABLE (ARG,VAL).
C           EPS    - AN INPUT CONSTANT WHICH IS USED AS UPPER BOUND
C                    FOR THE ABSOLUTE ERROR.
C           IER    - A RESULTING ERROR PARAMETER.
C
C        REMARKS
C           (1) TABLE (ARG,VAL) SHOULD REPRESENT A SINGLE-VALUED
C               FUNCTION AND SHOULD BE STORED IN SUCH A WAY, THAT THE
C               DISTANCES ABS(ARG(I)-X) INCREASE WITH INCREASING
C               SUBSCRIPT I. TO GENERATE THIS ORDER IN TABLE (ARG,VAL),
C               SUBROUTINES ATSG, ATSM OR ATSE COULD BE USED IN A
C               PREVIOUS STAGE.
C           (2) NO ACTION BESIDES ERROR MESSAGE IN CASE NDIM LESS
C               THAN 1.
C           (3) INTERPOLATION IS TERMINATED EITHER IF THE DIFFERENCE
C               BETWEEN TWO SUCCESSIVE INTERPOLATED VALUES IS
C               ABSOLUTELY LESS THAN TOLERANCE EPS, OR IF THE ABSOLUTE
C               VALUE OF THIS DIFFERENCE STOPS DIMINISHING, OR AFTER
C               (NDIM-1) STEPS (THE NUMBER OF POSSIBLE STEPS IS
C               DIMINISHED IF AT ANY STAGE INFINITY ELEMENT APPEARS IN
C               THE DOWNWARD DIAGONAL OF INVERTED-DIFFERENCES-SCHEME
C               AND IF IT IS IMPOSSIBLE TO ELIMINATE THIS INFINITY
C               ELEMENT BY INTERCHANGING OF TABLE POINTS).
C               FURTHER IT IS TERMINATED IF THE PROCEDURE DISCOVERS TWO
C               ARGUMENT VALUES IN VECTOR ARG WHICH ARE IDENTICAL.
C               DEPENDENT ON THESE FOUR CASES, ERROR PARAMETER IER IS
C               CODED IN THE FOLLOWING FORM
C                IER=0 - IT WAS POSSIBLE TO REACH THE REQUIRED
C                        ACCURACY (NO ERROR).
C                IER=1 - IT WAS IMPOSSIBLE TO REACH THE REQUIRED
C                        ACCURACY BECAUSE OF ROUNDING ERRORS.
C                IER=2 - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE
C                        NDIM IS LESS THAN 2, OR THE REQUIRED ACCURACY
C                        COULD NOT BE REACHED BY MEANS OF THE GIVEN
C                        TABLE. NDIM SHOULD BE INCREASED.
C                IER=3 - THE PROCEDURE DISCOVERED TWO ARGUMENT VALUES
C                        IN VECTOR ARG WHICH ARE IDENTICAL.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           INTERPOLATION IS DONE BY CONTINUED FRACTIONS AND INVERTED-
C           DIFFERENCES-SCHEME. ON RETURN Y CONTAINS AN INTERPOLATED
C           FUNCTION VALUE AT POINT X, WHICH IS IN THE SENSE OF REMARK
C           (3) OPTIMAL WITH RESPECT TO GIVEN TABLE. FOR REFERENCE, SEE
C           F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,
C           MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.395-406.
C
C     ..................................................................
C
      SUBROUTINE ACFI(X,ARG,VAL,Y,NDIM,EPS,IER)
C
C
      DIMENSION ARG(1),VAL(1)
      IER=2
      IF(NDIM)20,20,1
    1 Y=VAL(1)
      DELT2=0.
      IF(NDIM-1)20,20,2
C
C     PREPARATIONS FOR INTERPOLATION LOOP
    2 P2=1.
      P3=Y
      Q2=0.
      Q3=1.
C
C
C     START INTERPOLATION LOOP
      DO 16 I=2,NDIM
      II=0
      P1=P2
      P2=P3
      Q1=Q2
      Q2=Q3
      Z=Y
      DELT1=DELT2
      JEND=I-1
C
C     COMPUTATION OF INVERTED DIFFERENCES
    3 AUX=VAL(I)
      DO 10 J=1,JEND
      H=VAL(I)-VAL(J)
      IF(ABS(H)-1.E-6*ABS(VAL(I)))4,4,9
    4 IF(ARG(I)-ARG(J))5,17,5
    5 IF(J-JEND)8,6,6
C
C     INTERCHANGE ROW I WITH ROW I+II
    6 II=II+1
      III=I+II
      IF(III-NDIM)7,7,19
    7 VAL(I)=VAL(III)
      VAL(III)=AUX
      AUX=ARG(I)
      ARG(I)=ARG(III)
      ARG(III)=AUX
      GOTO 3
C
C     COMPUTATION OF VAL(I) IN CASE VAL(I)=VAL(J) AND J LESS THAN I-1
    8 VAL(I)=1.7E38                                                             0
      GOTO 10
C
C     COMPUTATION OF VAL(I) IN CASE VAL(I) NOT EQUAL TO VAL(J)
    9 VAL(I)=(ARG(I)-ARG(J))/H
   10 CONTINUE
C     INVERTED DIFFERENCES ARE COMPUTED
C
C     COMPUTATION OF NEW Y
      P3=VAL(I)*P2+(X-ARG(I-1))*P1
      Q3=VAL(I)*Q2+(X-ARG(I-1))*Q1
      IF(Q3)11,12,11
   11 Y=P3/Q3
      GOTO 13
   12 Y=1.7E38                                                                  0
   13 DELT2=ABS(Z-Y)
      IF(DELT2-EPS)19,19,14
   14 IF(I-8)16,15,15
   15 IF(DELT2-DELT1)16,18,18
   16 CONTINUE
C     END OF INTERPOLATION LOOP
C
C
      RETURN
C
C     THERE ARE TWO IDENTICAL ARGUMENT VALUES IN VECTOR ARG
   17 IER=3
      RETURN
C
C     TEST VALUE DELT2 STARTS OSCILLATING
   18 Y=Z
      IER=1
      RETURN
C
C     THERE IS SATISFACTORY ACCURACY WITHIN NDIM-1 STEPS
   19 IER=0
   20 RETURN
      END