        SUBROUTINE XPENTP(JDIM,KDIM,A,B,C,D,E,U1,U2,U3,U4,F,
     *                     JL,JU,KL,KU)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAXJ=165, MAXK=48, MAXJK=MAXJ*MAXK  )
        DOUBLE PRECISION LD,LD1,LD2,LDI,LDD,LU
        DIMENSION A(JDIM,KDIM),B(JDIM,KDIM),C(JDIM,KDIM),D(JDIM,KDIM)
        DIMENSION E(JDIM,KDIM),F(JDIM,KDIM),U1(JDIM,KDIM)
        DIMENSION U2(JDIM,KDIM),U3(JDIM,KDIM),U4(JDIM,KDIM)
C
        COMMON/WORKSP/LDA(MAXJ),LDB(MAXJ),LUS(MAXJ),
     *         LDS(MAXJ),WORK(MAXJ,88)
        DOUBLE PRECISION LDA,LDB,LUS,LDS
C
C        ! START FORWARD GENERATION PROCESS AND SWEEP
C
        I = JL
        DO 10 K = KL,KU
        LD = C(I,K)
        LDI = 1./LD
        U1(I,K) = D(I,K)*LDI
        U2(I,K) = E(I,K)*LDI
        U3(I,K) = A(I,K)*LDI
        U4(I,K) = B(I,K)*LDI
        F(I,K) = F(I,K)*LDI
10      CONTINUE
C
        I = JL+1
        DO 20 K = KL,KU
        LD1 = B(I,K)
        LD = C(I,K) - LD1*U1(I-1,K)
        LDI = 1./LD
        U1(I,K) = (D(I,K) - LD1*U2(I-1,K))*LDI
        U2(I,K) = E(I,K)*LDI
        U3(I,K) = -LD1*U3(I-1,K)*LDI
        U4(I,K) = (A(I,K) -LD1*U4(I-1,K))*LDI
        F(I,K) = (F(I,K) - LD1*F(I-1,K))*LDI
20      CONTINUE
C
          DO 1 I = JL+2,JU-4
                DO 11 K = KL,KU
                LD2 = A(I,K)
                LD1 = B(I,K) - LD2*U1(I-2,K)
                LD = C(I,K) - (LD2*U2(I-2,K) + LD1*U1(I-1,K))
                LDI = 1./LD
                U1(I,K) = (D(I,K) - LD1*U2(I-1,K))*LDI
                U2(I,K) = E(I,K)*LDI
                U3(I,K) = (-LD2*U3(I-2,K) - LD1*U3(I-1,K))*LDI
                U4(I,K) = (-LD2*U4(I-2,K) - LD1*U4(I-1,K))*LDI
                F(I,K) = (F(I,K) - LD2*F(I-2,K) - LD1*F(I-1,K))*LDI
11              CONTINUE
1           CONTINUE
C
                I = JU-3
                DO 12 K = KL,KU
                LD2 = A(I,K)
                LD1 = B(I,K) - LD2*U1(I-2,K)
                LD = C(I,K) - (LD2*U2(I-2,K) + LD1*U1(I-1,K))
                LDI = 1./LD
                U1(I,K) = (D(I,K) - LD1*U2(I-1,K))*LDI
                U3(I,K) = (E(I,K) - LD2*U3(I-2,K) - LD1*U3(I-1,K))*LDI
                U4(I,K) = (-LD2*U4(I-2,K) - LD1*U4(I-1,K))*LDI
                F(I,K) = (F(I,K) - LD2*F(I-2,K) - LD1*F(I-1,K))*LDI
12              CONTINUE
C
        I = JU-2
        DO 22 K = KL,KU
        LD2 = A(I,K)
        LD1 = B(I,K) - LD2*U1(I-2,K)
        LD = C(I,K) - (LD2*U2(I-2,K) + LD1*U1(I-1,K))
        LDI = 1./LD
        U3(I,K) = (D(I,K) - LD2*U3(I-2,K) - LD1*U3(I-1,K))*LDI
        U4(I,K) = (E(I,K) - LD2*U4(I-2,K) - LD1*U4(I-1,K))*LDI
        F(I,K) = (F(I,K) - LD2*F(I-2,K) - LD1*F(I-1,K))*LDI
22      CONTINUE
C
        I = JU-1
        DO 32 K = KL,KU
        LDA(K) = E(I,K)
        LDB(K) = -LDA(K)*U1(JL,K)
        F(I,K) = F(I,K) - LDA(K)*F(JL,K) - LDB(K)*F(JL+1,K)
        LDS(K) = LDA(K)*U3(JL,K) + LDB(K)*U3(JL+1,K)
        U4(I,K) = LDA(K)*U4(JL,K) + LDB(K)*U4(JL+1,K)
32      CONTINUE

        DO 3 M = JL+2,JU-4
           DO 14 K = KL,KU
           LDD = -LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)
           LDA(K) = LDB(K)
           LDB(K) = LDD
           LDS(K) = LDS(K) + LDD*U3(M,K)
           U4(I,K) = U4(I,K) + LDD*U4(M,K)
           F(I,K) = F(I,K) - LDD*F(M,K)
14         CONTINUE
3       CONTINUE
C
        M = JU-3
        DO 15 K = KL,KU
        LDD = A(I,K) - LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)
        LDA(K) = LDB(K)
        LDB(K) = LDD
        LDS(K) = LDS(K) + LDD*U3(M,K)
        U4(I,K) = U4(I,K) + LDD*U4(M,K)
        F(I,K) = F(I,K) - LDD*F(M,K)
15      CONTINUE
        M = JU-2
        DO 25 K = KL,KU
        LDD = B(I,K) - LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)
        LDA(K) = LDB(K)
        LDB(K) = LDD
        LDS(K) = LDS(K)+ LDD*U3(M,K)
        U4(I,K) = U4(I,K) + LDD*U4(M,K)
        F(I,K) = F(I,K) - LDD*F(M,K)
25      CONTINUE
C
        DO 35 K = KL,KU
        LDS(K) = C(I,K) - LDS(K)
        LDI = 1./LDS(K)
        U4(I,K) = (D(I,K) - U4(I,K))*LDI
        F(I,K) = F(I,K)*LDI
35      CONTINUE
C
C    THIS BREAKUP OF K LOOPS IS VERY IMPORTANT
C
        I = JU
        DO 16 K = KL,KU
        LDA(K) = D(I,K)
        LDB(K) = E(I,K) - LDA(K)*U1(JL,K)
        LUS(K) = LDA(K)*U3(JL,K) + LDB(K)*U3(JL+1,K)
        LDS(K) = LDA(K)*U4(JL,K) + LDB(K)*U4(JL+1,K)
        F(I,K) = F(I,K) - LDA(K)*F(JL,K) - LDB(K)*F(JL+1,K)
C
16      CONTINUE
C
        DO 4 M = JL+2,JU-3
           DO 17 K = KL,KU
           LDD = -LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)
           LDA(K) = LDB(K)
           LDB(K) = LDD
           LUS(K) = LUS(K) + LDD*U3(M,K)
           LDS(K) = LDS(K) + LDD*U4(M,K)
           F(I,K) = F(I,K) - LDD*F(M,K)
17         CONTINUE
4       CONTINUE
C
        M = JU-2
        DO 18 K = KL,KU
        LDD = A(I,K) - LDA(K)*U2(M-2,K) - LDB(K)*U1(M-1,K)
        LDA(K) = LDB(K)
        LDB(K) = LDD
        LUS(K) = LUS(K) + LDD*U3(M,K)
        LDS(K) = LDS(K) + LDD*U4(M,K)
        F(I,K) = F(I,K) - LDD*F(M,K)
18      CONTINUE
        M = JU-1
        DO 28 K = KL,KU
        LDD = B(I,K) - LUS(K)
        F(I,K) = F(I,K) - LDD*F(M,K)
        LDS(K) = C(I,K) - LDD*U4(M,K) - LDS(K)
        LDI = 1./LDS(K)
        F(I,K) = F(I,K)*LDI
28      CONTINUE
C
C        !  BACK SWEEP SOLUTION
C
      DO 19 K = KL,KU
      F(JU,K) = F(JU,K)
      F(JU-1,K) = F(JU-1,K) - U4(JU-1,K)*F(JU,K)
      F(JU-2,K) = F(JU-2,K) - U4(JU-2,K)*F(JU,K) - U3(JU-2,K)*F(JU-1,K)
      F(JU-3,K) = F(JU-3,K) - U4(JU-3,K)*F(JU,K) - U3(JU-3,K)*F(JU-1,K)
     *                - U1(JU-3,K)*F(JU-2,K)
19      CONTINUE
C
        DO 2 I = 4,JU-JL
        IX = JU-I
          DO 50 K = KL,KU
          F(IX,K) = F(IX,K) - U4(IX,K)*F(JU,K) - U3(IX,K)*F(JU-1,K)
     *              - U1(IX,K)*F(IX+1,K) - U2(IX,K)*F(IX+2,K)
50        CONTINUE
2       CONTINUE
C
        RETURN
        END
