SUBROUTINE PREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
     1   F, JAC, IERR)
CLLL. OPTIMIZE
      EXTERNAL F, JAC
      INTEGER NEQ, NYH, IWM
      INTEGER IOWND, IOWNS,
     1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
     2   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
      INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP,
     1   MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1
      DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WM 
      DOUBLE PRECISION ROWNS, 
     1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
      DOUBLE PRECISION CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ,
     1   VNORM
      DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*),
     1   WM(*), IWM(*)
      COMMON /LS0001/ ROWNS(209),
     2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
     3   IOWND(14), IOWNS(6), 
     4   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
     5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
C-----------------------------------------------------------------------
C PREPJ IS CALLED BY STODE TO COMPUTE AND PROCESS THE MATRIX
C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF
C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
C J IS STORED IN WM AND REPLACED BY P.  IF MITER .NE. 3, P IS THEN
C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
C BY DGETRF IF MITER = 1 OR 2, AND BY DGBTRF IF MITER = 4 OR 5.
C
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
C WITH PREPJ USES THE FOLLOWING..
C Y     = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
C FTEM  = WORK ARRAY OF LENGTH N (ACOR IN STODE). 
C SAVF  = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
C WM    = REAL WORK SPACE FOR MATRICES.  ON OUTPUT IT CONTAINS THE
C         INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
C         OF P IF MITER IS 1, 2 , 4, OR 5.
C         STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
C         WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
C         WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
C         WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
C IWM   = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
C         IWM(21), IF MITER IS 1, 2, 4, OR 5.  IWM ALSO CONTAINS BAND 
C         PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
C EL0   = EL(1) (INPUT).
C IERPJ = OUTPUT ERROR FLAG,  = 0 IF NO TROUBLE, .GT. 0 IF
C         P MATRIX FOUND TO BE SINGULAR.
C JCUR  = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX
C         (OR APPROXIMATION) IS NOW CURRENT.
C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
C MITER, N, NFE, AND NJE.
C-----------------------------------------------------------------------
      NJE = NJE + 1 
      IERPJ = 0
      JCUR = 1
      HL0 = H*EL0
      GO TO (100, 200, 300, 400, 500), MITER
C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
 100  LENP = N*N
      DO 110 I = 1,LENP
 110    WM(I+2) = 0.0D0
      CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N)
      CON = -HL0
      DO 120 I = 1,LENP
 120    WM(I+2) = WM(I+2)*CON 
      GO TO 240
C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
 200  FAC = VNORM (N, SAVF, EWT)
      R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC  
      IF (R0 .EQ. 0.0D0) R0 = 1.0D0
      SRUR = WM(1)
      J1 = 2
      DO 230 J = 1,N
        YJ = Y(J)
        R = DMAX1(SRUR*DABS(YJ),R0/EWT(J))
        Y(J) = Y(J) + R
        FAC = -HL0/R
        IERR = 0
        CALL F (NEQ, TN, Y, FTEM, IERR)
        IF (IERR .LT. 0) RETURN
        DO 220 I = 1,N
 220      WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
        Y(J) = YJ
        J1 = J1 + N 
 230    CONTINUE
      NFE = NFE + N 
C ADD IDENTITY MATRIX. -------------------------------------------------
 240  J = 3
      NP1 = N + 1
      DO 250 I = 1,N
        WM(J) = WM(J) + 1.0D0 
 250    J = J + NP1 
C DO LU DECOMPOSITION ON P. --------------------------------------------
      CALL DGETRF ( N, N, WM(3), N, IWM(21), IER)
      IF (IER .NE. 0) IERPJ = 1
      RETURN
C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
 300  WM(2) = HL0
      R = EL0*0.1D0 
      DO 310 I = 1,N
 310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
      IERR = 0
      CALL F (NEQ, TN, Y, WM(3), IERR)
      IF (IERR .LT. 0) RETURN
      NFE = NFE + 1 
      DO 320 I = 1,N
        R0 = H*SAVF(I) - YH(I,2)
        DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I))
        WM(I+2) = 1.0D0
        IF (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320
        IF (DABS(DI) .EQ. 0.0D0) GO TO 330
        WM(I+2) = 0.1D0*R0/DI 
 320    CONTINUE
      RETURN
 330  IERPJ = 1
      RETURN
C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
 400  ML = IWM(1)
      MU = IWM(2)
      ML3 = ML + 3
      MBAND = ML + MU + 1
      MEBAND = MBAND + ML
      LENP = MEBAND*N
      DO 410 I = 1,LENP
 410    WM(I+2) = 0.0D0
      CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND)
      CON = -HL0
      DO 420 I = 1,LENP
 420    WM(I+2) = WM(I+2)*CON 
      GO TO 570
C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
 500  ML = IWM(1)
      MU = IWM(2)
      MBAND = ML + MU + 1
      MBA = MIN0(MBAND,N)
      MEBAND = MBAND + ML
      MEB1 = MEBAND - 1
      SRUR = WM(1)
      FAC = VNORM (N, SAVF, EWT)
      R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC  
      IF (R0 .EQ. 0.0D0) R0 = 1.0D0
      DO 560 J = 1,MBA
        DO 530 I = J,N,MBAND
          YI = Y(I) 
          R = DMAX1(SRUR*DABS(YI),R0/EWT(I))
 530      Y(I) = Y(I) + R
        IERR = 0
        CALL F (NEQ, TN, Y, FTEM, IERR)
        IF (IERR .LT. 0) RETURN
        DO 550 JJ = J,N,MBAND 
          Y(JJ) = YH(JJ,1)
          YJJ = Y(JJ)
          R = DMAX1(SRUR*DABS(YJJ),R0/EWT(JJ))
          FAC = -HL0/R
          I1 = MAX0(JJ-MU,1)
          I2 = MIN0(JJ+ML,N)
          II = JJ*MEB1 - ML + 2
          DO 540 I = I1,I2
 540        WM(II+I) = (FTEM(I) - SAVF(I))*FAC
 550      CONTINUE
 560    CONTINUE
      NFE = NFE + MBA
C ADD IDENTITY MATRIX. -------------------------------------------------
 570  II = MBAND + 2
      DO 580 I = 1,N
        WM(II) = WM(II) + 1.0D0
 580    II = II + MEBAND
C DO LU DECOMPOSITION OF P. --------------------------------------------
      CALL DGBTRF ( N, N, ML, MU, WM(3), MEBAND, IWM(21), IER)
      IF (IER .NE. 0) IERPJ = 1
      RETURN
C----------------------- END OF SUBROUTINE PREPJ -----------------------
      END