      SUBROUTINE ADISET(IM,IN,IMU,LM,LX,RM,RX,KK,IFN,IFUT,INFO,PARM,ARR)
 
C***********************************************************************
C  SUBROUTINE SETS UP COEFFICIENT ARRAYS FOR 2 STATE ADI ALGORITHM
C  ARGUMENTS: IM       NUMBER OF GRID INTERVALS IN L DIRECTION
C             IN       NUMBER OF GRID INTERVALS IN R DIRECTION
C             IMU      FIRST DIMENSION - 1 OF U SOLN. MATRIX IF
C                      DIFFERENT FROM IM (DEFAULTS TO IM)
C             LM       MIN. VALUE OF STATE VARIABLE L
C             LX       MAX. VALUE OF STATE VARIABLE L
C             RM       MIN. VALUE OF STATE VARIABLE R
C             RX       MAX. VALUE OF STATE VARIABLE R
C             KK       STEP SIZE IN TIME DIRECTION
C             IFN      MODEL NUMBER (IFNSET) SELECTING COEFFICIENT
C             IFUT     FLAG SET TO 1 FOR FUTURES PRICES
C             INFO()   FLAGS FOR BOUNDARY CONDITIONS (IULMN,IULMX,
C                         IURMN,IURMX,ILMNRN,ILMNRX,ILMXRN,ILMXRX)
C             PARM()   VECTOR OF MODEL PARAMETERS FOR CALCULATING COEF
C             ARRAY()  OUTPUT ARRAY OF COEFFICIENTS SET FOR ADI STEP
C                      ARRAY() CORRESPONDS TO BK1001 AS FOLLOWS:
C                      ARR(1,.,.) = D1MTX(.,.)   ARR(4,.,.) = DP1MTX()
C                      ARR(2,.,.) = D2MTX(.,.)   ARR(5,.,.) = DP2MTX()
C                      ARR(3,.,.) = D3MTX(.,.)   ARR(6,.,.) = DP3MTX()
C                      ARR(7,.,.) = C1MTX(.,.)   ARR(8,.,.) = C2MTX ()
C
C  AUTHOR:             ROBERT JONES    3 OCTOBER 1988
C
C***********************************************************************
 
      IMPLICIT   INTEGER (I-J,M-N)
      IMPLICIT   DOUBLE PRECISION (A-H,K-L,O-Z)
      DIMENSION  INFO(8), PARM(15)
      DIMENSION  ARR( 8, 0:IM, 0:IN )
      COMMON    /ADICOM/ A,B,C,D,E,F,K,        X1,X2,X3,X4,X5,
     1                   L,LMIN,LMAX,LSCALE,   R,RMIN,RMAX,RSCALE,
     2                   IULMN,IULMX,IURMN,IURMX,ILMNRN,ILMNRX,
     3                   ILMXRN,ILMXRX,IFNSET,IFUTUR,   M,N,MU
      SAVE       ADICOM
 
      M      = IM
      N      = IN
      K      = KK
      MU     = MAX( IM, IMU )
      IFNSET = IFN
      IFUTUR = IFUT
      IULMN  = INFO(1)
      IULMX  = INFO(2)
      IURMN  = INFO(3)
      IURMX  = INFO(4)
      ILMNRN = INFO(5)
      ILMNRX = INFO(6)
      ILMXRN = INFO(7)
      ILMXRX = INFO(8)
      LMIN   = LM
      LMAX   = LX
      RMIN   = RM
      RMAX   = RX
 
      DO 100 I3 = 0, N
         DO 100 I2 = 0, M
         DO 100 I1 = 1, 8
 100     ARR(I1,I2,I3) = 0.0
 
      LSCALE = DBLE(M) / ( LMAX - LMIN )
      RSCALE = DBLE(N) / ( RMAX - RMIN )
      X1     = LSCALE * LSCALE * K
      X2     = RSCALE * RSCALE * K
      X3     = LSCALE * RSCALE * K
      X4     = LSCALE * K
      X5     = RSCALE * K
 
C**** COMPUTE INTERIOR COEFFICIENTS
      DO 200 I = 1, M-1
         DO 200 J = 1, N-1
            CALL ADI1(I,J,PARM)
            ARR(1,I,J)  =  D / 4.0 - A / 2.0
            ARR(2,I,J)  =  1.0 + A
            ARR(3,I,J)  = -D / 4.0 - A / 2.0
            ARR(7,I,J)  =  C / 4.0
            ARR(8,I,J)  =  F - 2.0 * B - A
            ARR(4,I,J)  =  E / 4.0 - B / 2.0
            ARR(5,I,J)  =  1.0 + B
            ARR(6,I,J)  = -E / 4.0 - B / 2.0
 200  CONTINUE
 
C**** COMPUTE COEFFICIENTS AT LMIN
      IF     ( IULMN .EQ. 1 ) THEN
C        KNOWN SOLUTION AT LMIN
         DO 220 J = 1, N-1
 220        ARR(2,0,J) = 1.0
 
      ELSE IF ( IULMN .EQ. 2 ) THEN
C        SOLN AT LMIN SATISFIES 1 STATE VARIABLE PDE
         DO 240 J = 1, N-1
 240        ARR(2,0,J) = 1.0
         CALL ADI2(PARM,ARR)
 
      ELSE
C        QUADRATIC EXTRAPOL.
         DO 260 J = 1, N-1
            G    =  ( 3.0 + ARR(2,2,J) / ARR(3,2,J) ) / ARR(3,1,J)
            ARR(3,0,J) = -3.0 + ARR(1,2,J)/ARR(3,2,J)-G*ARR(2,1,J)
 260        ARR(2,0,J) =  1.0 - G * ARR(1,1,J)
         CALL ADI2(PARM,ARR)
      END IF
 
C**** COMPUTE COEFFICIENTS AT LMAX
      IF ( IULMX .EQ. 1 ) THEN
C        KNOWN SOLUTION AT L = MAX
         DO 310 J= 1,N-1
 310        ARR(2,M,J)= 1.0
 
      ELSE
C        FOR QUADRATIC EXTRAPOLATION
         DO 320 J = 1, N-1
            G = (3.0+ARR(2,M-2,J)/ARR(1,M-2,J)) / ARR(1,M-1,J)
            ARR(1,M,J)= -3.0+ARR(3,M-2,J)/ARR(1,M-2,J)-G*ARR(2,M-1,J)
 320        ARR(2,M,J)=  1.0 - G * ARR(3,M-1,J)
         DO 330  J= 1,N-1
            CALL ADI1(M,J,PARM)
            ARR(4,M,J) =  E/4.0 - B/2.0
            ARR(5,M,J) =  1.0 + B
            ARR(6,M,J) = -E/4.0 - B/2.0
C           USED RHS OF BOUNDARY PDE
 330        ARR(8,M,J) =  B - F
 
         IF ( ILMXRN.EQ.1) THEN
C           SOLN. KNOWN AT L=LMAX, R=0:
            ARR(5,M,0) = 1.0
         ELSE
            G = ( 3.0 + ARR(5,M,2)/ARR(6,M,2) ) / ARR(6,M,1)
            ARR(6,M,0) = -3.0 + ARR(4,M,2)/ARR(6,M,2) -G * ARR(5,M,1)
            ARR(5,M,0) =  1.0 - G * ARR(4,M,1)
         END IF
 
         IF ( ILMXRX .EQ. 1 ) THEN
C           KNOWN SOLUTION AT L=LMAX, R=RMAX:
            ARR(5,M,N) = 1.0
         ELSE
C           QUADRATIC EXTRAPOLATION:
            G = ( 3.0 + ARR(5,M,N-2)/ARR(4,M,N-2) ) / ARR(4,M,N-1)
            ARR(4,M,N) =-3.0+ARR(6,M,N-2)/ARR(4,M,N-2)-G*ARR(5,M,N-1)
            ARR(5,M,N) = 1.0 - G * ARR(6,M,N-1)
         END IF
      END IF
 
C**** COMPUTE COEFFICIENTS AT RMIN
      IF     ( IURMN .EQ. 1 ) THEN
C        KNOWN SOLUTION AT R=RMIN:
         DO 410 I = 1, M-1
 410        ARR(5,I,0) = 1.0
 
      ELSE IF ( IURMN .EQ. 2 ) THEN
C        FOR BOUNDARY PDE SOLUTION:
         DO 420 I = 1, M-1
 420        ARR(5,I,0) = 1.0
 
         DO 430 I = 1, M-1
            CALL ADI1(I,0,PARM)
            ARR(1,I,0) =  D/4.0 - A/2.0
            ARR(2,I,0) =  1.0 + A
            ARR(3,I,0) = -D/4.0 - A/2.0
C           USED RHS OF BOUNDARY PDE
 430        ARR(8,I,0) =  A - F
 
         IF (ILMNRN.EQ.1) THEN
C           KNOWN SOLUTION AT L=LMIN, R=RMIN:
             ARR(2,0,0)= 1.0
          ELSE
             G= (3.0+ARR(2,2,0)/ARR(3,2,0))/ARR(3,1,0)
             ARR(3,0,0)= -3.0+ARR(1,2,0)/ARR(3,2,0)-G*ARR(2,1,0)
             ARR(2,0,0)= 1.0-G*ARR(1,1,0)
          END IF
 
          IF ( ILMXRN .EQ. 1 ) THEN
C            KNOWN SOLUTION AT L=LMAX, R=RMIN:
             ARR(2,M,0) = 1.0
          ELSE
             G = (3.0+ARR(2,M-2,0)/ARR(1,M-2,0))/ARR(1,M-1,0)
             ARR(1,M,0) = -3.0+ARR(3,M-2,0)/ARR(1,M-2,0)-G
     /                       *ARR(2,M-1,0)
             ARR(2,M,0) = 1.0 - G * ARR(3,M-1,0)
          END IF
 
       ELSE
C         QUADRATIC EXTRAPOLATION:
          DO 440 I = 1, M-1
              G = (3.0+ARR(5,I,2)/ARR(6,I,2))/ARR(6,I,1)
              ARR(6,I,0) = -3.0 + ARR(4,I,2)/ARR(6,I,2)- G*ARR(5,I,1)
 440          ARR(5,I,0) = 1.0 - G * ARR(4,I,1)
      END IF
 
C**** COMPUTE COEFFICIENTS AT RMAX
      IF ( IURMX .EQ. 1 ) THEN
C        KNOWN SOLUTION AT R=MAX:
         DO 500 I = 1, M-1
 500        ARR(5,I,N) = 1.0
      ELSE
C        QUADRATIC EXTRAPOLATION:
         DO 510 I = 1, M-1
            G = ( 3.0 + ARR(5,I,N-2)/ARR(4,I,N-2) ) / ARR(4,I,N-1)
            ARR(4,I,N) = -3.0 + ARR(6,I,N-2)
     /                   / ARR(4,I,N-2) - G * ARR(5,I,N-1)
 510        ARR(5,I,N) = 1.0 - G * ARR(6,I,N-1)
      END IF
      RETURN
      END
 
C**** SUBROUTINES CALLED BY ADISET() **********************************
 
C**** SUBROUTINE TO CALCULATE & SCALE PDE COEFFICIENTS AT (I,J):
      SUBROUTINE ADI1(I,J,PARM)
      IMPLICIT   INTEGER (I-J,M-N)
      IMPLICIT   DOUBLE PRECISION (A-H,K-L,O-Z)
      DIMENSION  PARM(15)
      COMMON    /ADICOM/ A,B,C,D,E,F,K,        X1,X2,X3,X4,X5,
     1                   L,LMIN,LMAX,LSCALE,   R,RMIN,RMAX,RSCALE,
     2                   IULMN,IULMX,IURMN,IURMX,ILMNRN,ILMNRX,
     3                   ILMXRN,ILMXRX,IFNSET,IFUTUR,   M,N,MU
      SAVE       ADICOM

         L  = LMIN + DBLE(I) / LSCALE
         R  = RMIN + DBLE(J) / RSCALE
         IF ( IFUTUR .EQ. 1 )  THEN
            F = 0.0
         ELSE
            F = FNAKF(IFNSET,L,R,PARM)
         END IF
         K1 = 1.0 - F * K / 2.0
         A  = X1 * FNAKA(IFNSET,L,R,PARM) / K1
         B  = X2 * FNAKB(IFNSET,L,R,PARM) / K1
         C  = X3 * FNAKC(IFNSET,L,R,PARM) / K1
         D  = X4 * FNAKD(IFNSET,L,R,PARM) / K1
         E  = X5 * FNAKE(IFNSET,L,R,PARM) / K1
         F  = (2.0 - K1) / K1
      RETURN
      END
 
C**** FORMERLY SUBROUTINE AK2014(IFUTUR,PARM)
      SUBROUTINE ADI2(PARM,ARR)
      IMPLICIT   INTEGER (I-J,M-N)
      IMPLICIT   DOUBLE PRECISION (A-H,K-L,O-Z)
      COMMON    /ADICOM/ A,B,C,D,E,F,K,        X1,X2,X3,X4,X5,
     1                   L,LMIN,LMAX,LSCALE,   R,RMIN,RMAX,RSCALE,
     2                   IULMN,IULMX,IURMN,IURMX,ILMNRN,ILMNRX,
     3                   ILMXRN,ILMXRX,IFNSET,IFUTUR,   M,N,MU
      SAVE       ADICOM
      DIMENSION  PARM(15), ARR( 8 , 0:M , 0:* )
 
      DO 100 J = 1, N-1
         CALL ADI1(0,J,PARM)
         ARR(4,0,J) =  E/4.0 - B/2.0
         ARR(5,0,J) =  1.0 + B
         ARR(6,0,J) = -E/4.0 - B/2.0
C        USED RHS OF BOUNDARY
 100     ARR(8,0,J) = B-F
 
      IF ( ILMNRN .EQ. 1 ) THEN
C        IF GIVEN VALUE AT L=0, R=0:
         ARR(5,0,0) = 1.0
      ELSE
C        QUADRATIC EXTRAPLATION:
         G = ( 3.0 + ARR(5,0,2)/ARR(6,0,2) ) / ARR(6,0,1)
         ARR(6,0,0) = -3.0 + ARR(4,0,2)/ARR(6,0,2) - G * ARR(5,0,1)
         ARR(5,0,0) =  1.0 - G * ARR(4,0,1)
      END IF
 
      IF ( ILMNRX .EQ. 1 ) THEN
C        IF GIVEN VALUE AT L=0, R=RMAX:
         ARR(5,0,N) = 1.0
      ELSE
C        QUADRATIC EXTRAPLATION:
         G = ( 3.0 + ARR(5,0,N-2)/ARR(4,0,N-2) ) / ARR(4,0,N-1)
         ARR(4,0,N) = -3.0 + ARR(6,0,N-2)/ARR(4,0,N-2)-G*ARR(5,0,N-1)
         ARR(5,0,N) = 1.0 - G * ARR(6,0,N-1)
      END IF
      RETURN
      END
 
      SUBROUTINE ADSTEP(U,V,ARR)
 
C***********************************************************************
C  SUBROUTINE TAKES ONE TIME STEP USING 2 STATE ADI ALGORITHM.
C  CALL ADISET() FIRST TO SET UP NECESSARY COEFFICIENT ARRAYS ARR().
C  ARGUMENTS: U()      CURRENT SOLUTION ARRAY -- OVERWRITTEN BY ADSTEP
C             V()      INTERMEDIATE SOLUTION ARRAY
C             ARR()    INPUT ARRAY OF COEFFICIENTS SET BY ADISET
C                      ARRAY() CORRESPONDS TO BK1001 AS FOLLOWS:
C                      ARR(1,.,.) = D1MTX(.,.)   ARR(4,.,.) = DP1MTX()
C                      ARR(2,.,.) = D2MTX(.,.)   ARR(5,.,.) = DP2MTX()
C                      ARR(3,.,.) = D3MTX(.,.)   ARR(6,.,.) = DP3MTX()
C                      ARR(7,.,.) = C1MTX(.,.)   ARR(8,.,.) = C2MTX ()
C  AUTHOR:    ROBERT JONES       3 OCTOBER 1988
C***********************************************************************
 
      IMPLICIT INTEGER (I-J,M-N)
      IMPLICIT DOUBLE PRECISION (A-H,K-L,O-Z)
 
C     NOTE: INFORMATION IN ADICOM IS SET BY SUBROUTINE ADISET()
C           VARIABLES IN FIRST ROW OF ADICOM NOT USED IN ADSTEP()
      COMMON /ADICOM/ ZA,ZB,ZC,ZD,ZE,ZF,ZK, Z1,Z2,Z3,Z4,Z5,
     1                L,LMIN,LMAX,LSCALE,   R,RMIN,RMAX,RSCALE,
     2                IULMN,IULMX,IURMN,IURMX,ILMNRN,ILMNRX,
     3                ILMXRN,ILMXRX,IFNSET,IFUTUR,   M,N,MU
      SAVE    ADICOM
      
C     NOTE: PARAMETER NT MUST BE .GE. MAX(M,N) FOR TRIDAG ALGORITHM
      PARAMETER ( NT = 200 )
      COMMON /TRICOM/ A(0:NT),B(0:NT),C(0:NT),D(0:NT),GAM(0:NT),X(0:NT)
 
C     NOTE: ARR,U,V ARE THE ADJUSTABLE DIMENSION ARRAYS
      DIMENSION ARR(8, 0:M, 0:N), U(0:MU, 0:N), V(0:M, 0:N)
 
 
C**** FIRST SOLVE IN R DIRECTION ***************************************
 
      IF ( IULMN .EQ. 1 ) THEN
C        PUT LMIN BOUNDARY VALUES IN V(0,J):
         DO 100 J = 0, N
 100        V(0,J) = FNULMN( RMIN + DBLE(J)/RSCALE )
      ENDIF
 
      IF ( IULMN .EQ. 2 ) THEN
C        SOLVE BOUNDARY PDE AT L=0. PUT TRANSFORMED SOLN. IN V(0,J)
         DO 110 J = 0, N
            A(J) = ARR(4,0,J)
            B(J) = ARR(5,0,J)
 110        C(J) = ARR(6,0,J)
         DO 120 J=1,N-1
 120        D(J) = -A(J)*U(0,J-1) - ARR(8,0,J)*U(0,J) - C(J)*U(0,J+1)
 
         IF ( ILMNRN .EQ. 1 ) THEN
            D(0) = FNLNRN
         ELSE
            D(0) = D(2) / C(2) - D(1) * ( 3.0 + B(2) / C(2) ) / C(1)
         ENDIF
 
         IF ( ILMNRX .EQ. 1 ) THEN
            D(N) = FNLNRX
         ELSE
            D(N) = D(N-2)/A(N-2)-D(N-1)*(3.0 + B(N-2)/A(N-2))/A(N-1)
         ENDIF
 
         CALL ADI4(A,B,C,D,GAM,X,N)
 
C        TRANSFORM RESULT BACK TO INTERMEDIATE SOLN VALUES
         DO 130 J = 1, N-1
 130        V(0,J) = U(0,J) + ARR(4,0,J) * ( X(J-1) - U(0,J-1) )
     /               + ARR(5,0,J) * ( X(J) - U(0,J) )
     /               + ARR(6,0,J) * (X(J+1) - U(0,J+1))
         V(0,0) = U(0,0) + ARR(5,0,0) * (X(0)-U(0,0)) +
     /            ARR(6,0,0) * (X(1)-U(0,1))
         V(0,N) = U(0,N) + ARR(4,0,N) * (X(N-1) - U(0,N-1)) +
     /            ARR(5,0,N) * (X(N) - U(0,N))
      ENDIF
 
      IF ( IULMX .EQ. 1 ) THEN
C        PUT UMAX BOUNDARY VALUES IN V(M,J)
         DO 140 J = 0, N
 140        V(M,J) = FNULMX( RMIN + DBLE(J)/RSCALE )
      ENDIF
 
C     SOLVE FOR INTERMEDIATE SOLUTION VALUES
      DO 180  J = 1, N-1
         DO 150 I = 0, M
            A(I) = ARR(1,I,J)
            B(I) = ARR(2,I,J)
            C(I) = ARR(3,I,J)
 150     CONTINUE
 
         DO 160 I = 1, M-1
            D(I) = ( U(I+1,J+1) - U(I+1,J-1) - U(I-1,J+1) + U(I-1,J-1) )
     /                        * ARR(7,I,J)
     /             - U(I+1,J) * ARR(3,I,J)
     /             - U(I-1,J) * ARR(1,I,J)
     /             + U(I,J)   * ARR(8,I,J)
     /             - U(I,J+1) * ARR(6,I,J) * 2.0
     /             - U(I,J-1) * ARR(4,I,J) * 2.0
 160     CONTINUE
 
         IF ( IULMN .EQ. 1 .OR. IULMN .EQ. 2 ) THEN
            D(0) = V(0,J)
         ELSE
            D(0) = D(2) / C(2) - D(1) * ( 3.0 + B(2)/C(2) ) / C(1)
         ENDIF
 
         IF ( IULMX .EQ. 1 ) THEN
            D(M) = V(M,J)
         ELSE
            D(M) = D(M-2)/A(M-2)-D(M-1)*(3.0 + B(M-2)/A(M-2))/A(M-1)
         ENDIF
 
         CALL ADI4(A,B,C,D,GAM,X,M)
 
         DO 170  I = 0, M
 170        V(I,J) = X(I)
 180  CONTINUE
 
C**** ALTERNATE DIRECTION TO GET FINAL SOLUTION VALUES *****************
 
      IF ( IURMN .EQ. 1 ) THEN
C        PUT RMIN BOUNDARY VALUES IN V(I,0)
         DO 200 I = 0, M
 200        V(I,0) = FNURMN( LMIN + DBLE(I)/LSCALE )
      ENDIF
 
      IF ( IURMN .EQ. 2 ) THEN
C        SOLVE BOUNDARY PDE AT R=0. STORE SOLUTION IN V(I,0)
         DO 210 I = 0, M
            A(I) = ARR(1,I,0)
            B(I) = ARR(2,I,0)
 210        C(I) = ARR(3,I,0)
 
         DO 220  I=1,M-1
 220        D(I)= - A(I)*U(I-1,0) - ARR(8,I,0)*U(I,0) - C(I)*U(I+1,0)
 
         IF ( ILMNRN .EQ. 1 ) THEN
            D(0) = FNLNRN
         ELSE
            D(0) = D(2)/C(2) - D(1) * ( 3.0+B(2) / C(2) ) / C(1)
         ENDIF
 
         IF ( ILMXRN .EQ. 1 ) THEN
            D(M) = FNLXRN
         ELSE
            D(M) = D(M-2)/A(M-2) - D(M-1)*(3.0+B(M-2)/A(M-2)) / A(M-1)
         ENDIF
 
         CALL ADI4(A,B,C,D,GAM,X,M)
 
         DO 230  I = 0, M
 230        V(I,0) = X(I)
      ENDIF
 
      IF ( IURMX .EQ. 1 ) THEN
C        PUT RMAX BOUNDARY VALUES IN V(I,N)
         DO 240  I = 0, M
 240        V(I,N) = FNURMX( LMIN + DBLE(I)/LSCALE )
      ENDIF
 
C     SOLVE FOR FINAL U VALUES
      DO 250 I=1,M-1
 250     CALL ADI3(I,U,V,ARR)
 
      IF     ( IULMN .EQ. 1 ) THEN
C        COPY V TO U:
         DO 260 J = 0, N
 260     U(0,J) = V(0,J)
      ELSEIF ( IULMN .EQ. 2 ) THEN
         CALL ADI3(0,U,V,ARR)
      ELSE
C        QUADRATIC EXTRAPOLATION AT LMIN:
         DO 270  J = 0, N
 270     U(0,J) = 3.0 * U(1,J) - 3.0 * U(2,J) + U(3,J)
      ENDIF
 
      IF     ( IULMX .EQ. 1 ) THEN
C        COPY V VALUES TO U FOR I=M:
         DO 280 J = 0, N
 280     U(M,J) = V(M,J)
      ELSEIF ( IULMX .EQ. 2 ) THEN
         CALL ADI3(M,U,V,ARR)
      ELSE
C        QUADRATIC EXTRAPOLATION AT LMAX:
         DO 290 J = 0, N
 290     U(M,J) = 3.0 * U(M-1,J) - 3.0 * U(M-2,J) + U(M-3,J)
      ENDIF
 
      RETURN
      END
 
C**** SUBROUTINES CALLED BY ADSTEP (PREVIOUSLY BK4301 & BK4500 *********
 
      SUBROUTINE ADI3(I,U,V,ARR)
 
      IMPLICIT INTEGER (I-J,M-N)
      IMPLICIT DOUBLE PRECISION (A-H,K-L,O-Z)
      COMMON    /ADICOM/ ZA,ZB,ZC,ZD,ZE,ZF,ZK, Z1,Z2,Z3,Z4,Z5,
     1                   L,LMIN,LMAX,LSCALE,   R,RMIN,RMAX,RSCALE,
     2                   IULMN,IULMX,IURMN,IURMX,ILMNRN,ILMNRX,
     3                   ILMXRN,ILMXRX,IFNSET,IFUTUR,   M,N,MU
      SAVE      ADICOM
      DIMENSION ARR(8, 0:M, 0:N), U(0:MU, 0:N), V(0:M, 0:N)
 
C     NOTE: PARAMETER NT MUST BE .GE. MAX(M,N) FOR TRIDAG ALGORITHM
      PARAMETER ( NT = 200 )
      COMMON /TRICOM/ A(0:NT),B(0:NT),C(0:NT),D(0:NT),GAM(0:NT),X(0:NT)
 
      DO 100  J = 0, N
         A(J) = ARR(4,I,J)
         B(J) = ARR(5,I,J)
 100     C(J) = ARR(6,I,J)
 
      DO 110  J = 1, N-1
         D(J) = V(I,J) + ARR(4,I,J) * U(I,J-1)
     /        + (ARR(5,I,J) - 1.0) * U(I,J) + ARR(6,I,J) * U(I,J+1)
 110  CONTINUE
 
      IF ( IURMN .EQ. 1 .OR. IURMN .EQ. 2 ) THEN
         D(0) = V(I,0)
      ELSE
         D(0) = D(2)/C(2) - D(1) * ( 3.0 + B(2)/C(2) ) / C(1)
      ENDIF
 
      IF ( IURMX .EQ. 1 ) THEN
         D(N) = V(I,N)
      ELSE
         D(N) = D(N-2)/A(N-2) - D(N-1) * (3.0 + B(N-2)/A(N-2)) / A(N-1)
      ENDIF
 
      CALL ADI4(A,B,C,D,GAM,X,N)
      DO 120 J = 0, N
 120     U(I,J) = X(J)
 
      RETURN
      END
 
C**** TRIDIAGONAL SOLN. ALGORITHM FROM "NUMERICAL RECIPES", P. 40 ******
 
C     SOLVES: (ABC)X = D FOR X.  N=DIMENSION.  A,B,C,D, NOT ALTERED
C     NOTE: SUBSCRIPTS RUN FROM 0 AND SCRATCH VECTOR GAM VARIABLE DIMEN.
 
      SUBROUTINE ADI4(A,B,C,D,GAM,X,N)
      IMPLICIT INTEGER (I-J,M-N)
      IMPLICIT DOUBLE PRECISION (A-H,K-L,O-Z)
      DIMENSION A(0:*),B(0:*),C(0:*),D(0:*),GAM(0:*),X(0:*)
      BET = B(0)
      X(0) = D(0) / BET
      DO 11 J = 1, N
         GAM(J) = C(J-1) / BET
         BET = B(J) - A(J) * GAM(J)
         X(J) = (D(J) - A(J) * X(J-1)) / BET
  11  CONTINUE
      DO 12 J = N-1, 0, -1
  12     X(J) = X(J) - GAM(J+1) * X(J+1)
      RETURN
      END

