      SUBROUTINE PRJS (NEQ,Y,YH,NYH,EWT,FTEM,SAVF,WK,IWK,F,JAC)
CLLL. OPTIMIZE
      EXTERNAL F,JAC
      INTEGER NEQ, NYH, IWK
      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 IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
     1   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
     2   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
     3   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
      INTEGER I, IMUL, J, JJ, JOK, JMAX, JMIN, K, KMAX, KMIN, NG
      DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WK 
      DOUBLE PRECISION ROWNS, 
     1   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
      DOUBLE PRECISION CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH
      DOUBLE PRECISION CON, DI, FAC, HL0, PIJ, R, R0, RCON, RCONT,
     1   SRUR, VNORM
      DIMENSION NEQ(1), Y(1), YH(NYH,1), EWT(1), FTEM(1), SAVF(1),
     1   WK(1), IWK(1)
      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
      COMMON /LSS001/ CON0, CONMIN, CCMXJ, PSMALL, RBIG, SETH,
     1   IPLOST, IESP, ISTATC, IYS, IBA, IBIAN, IBJAN, IBJGP,
     2   IPIAN, IPJAN, IPJGP, IPIGP, IPR, IPC, IPIC, IPISP, IPRSP, IPA,
     3   LENYH, LENYHM, LENWK, LREQ, LRAT, LREST, LWMIN, MOSS, MSBJ,
     4   NSLJ, NGP, NLU, NNZ, NSP, NZL, NZU
C-----------------------------------------------------------------------
C PRJS IS CALLED TO COMPUTE AND PROCESS THE MATRIX
C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
C J IS COMPUTED BY COLUMNS, EITHER BY THE USER-SUPPLIED ROUTINE JAC
C IF MITER = 1, OR BY FINITE DIFFERENCING IF MITER = 2.
C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
C IF MITER = 1 OR 2, AND IF THE EXISTING VALUE OF THE JACOBIAN
C (AS CONTAINED IN P) IS CONSIDERED ACCEPTABLE, THEN A NEW VALUE OF
C P IS RECONSTRUCTED FROM THE OLD VALUE.  IN ANY CASE, WHEN MITER
C IS 1 OR 2, THE P MATRIX IS SUBJECTED TO LU DECOMPOSITION IN CDRV.
C P AND ITS LU DECOMPOSITION ARE STORED (SEPARATELY) IN WK. 
C
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
C WITH PRJS 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 WK    = REAL WORK SPACE FOR MATRICES.  ON OUTPUT IT CONTAINS THE
C         INVERSE DIAGONAL MATRIX IF MITER = 3, AND P AND ITS SPARSE
C         LU DECOMPOSITION IF MITER IS 1 OR 2.
C         STORAGE OF MATRIX ELEMENTS STARTS AT WK(3).
C         WK ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
C         WK(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
C         WK(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
C IWK   = INTEGER WORK SPACE FOR MATRIX-RELATED DATA, ASSUMED TO
C         BE EQUIVALENCED TO WK.  IN ADDITION, WK(IPRSP) AND IWK(IPISP)
C         ARE ASSUMED TO HAVE IDENTICAL LOCATIONS.
C EL0   = EL(1) (INPUT).
C IERPJ = OUTPUT ERROR FLAG (IN COMMON).
C       = 0 IF NO ERROR.
C       = 1  IF ZERO PIVOT FOUND IN CDRV.
C       = 2  IF A SINGULAR MATRIX AROSE WITH MITER = 3.
C       = -1 IF INSUFFICIENT STORAGE FOR CDRV (SHOULD NOT OCCUR HERE).
C       = -2 IF OTHER ERROR FOUND IN CDRV (SHOULD NOT OCCUR HERE).
C JCUR  = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX
C         (OR APPROXIMATION) IS NOW CURRENT.
C THIS ROUTINE ALSO USES OTHER VARIABLES IN COMMON.
C-----------------------------------------------------------------------
      HL0 = H*EL0
      CON = -HL0
      IF (MITER .EQ. 3) GO TO 300
C SEE WHETHER J SHOULD BE REEVALUATED (JOK = 0) OR NOT (JOK = 1). ------
      JOK = 1
      IF (NST .EQ. 0 .OR. NST .GE. NSLJ+MSBJ) JOK = 0
      IF (ICF .EQ. 1 .AND. DABS(RC - 1.0D0) .LT. CCMXJ) JOK = 0
      IF (ICF .EQ. 2) JOK = 0 
      IF (JOK .EQ. 1) GO TO 250
C
C MITER = 1 OR 2, AND THE JACOBIAN IS TO BE REEVALUATED. ---------------
 20   JCUR = 1
      NJE = NJE + 1 
      NSLJ = NST
      IPLOST = 0
      CONMIN = DABS(CON)
      GO TO (100, 200), MITER 
C
C IF MITER = 1, CALL JAC, MULTIPLY BY SCALAR, AND ADD IDENTITY. --------
 100  CONTINUE
      KMIN = IWK(IPIAN)
      DO 130 J = 1, N
        KMAX = IWK(IPIAN+J) - 1
        DO 110 I = 1,N
 110      FTEM(I) = 0.0D0
        CALL JAC (NEQ, TN, Y, J, IWK(IPIAN), IWK(IPJAN), FTEM)
        DO 120 K = KMIN, KMAX 
          I = IWK(IBJAN+K)
          WK(IBA+K) = FTEM(I)*CON
          IF (I .EQ. J) WK(IBA+K) = WK(IBA+K) + 1.0D0
 120      CONTINUE
        KMIN = KMAX + 1
 130    CONTINUE
      GO TO 290
C
C IF MITER = 2, MAKE NGP CALLS TO F TO APPROXIMATE J AND P. ------------
 200  CONTINUE
      FAC = VNORM(N, SAVF, EWT)
      R0 = 1000.0D0 * DABS(H) * UROUND * DBLE(N) * FAC    
      IF (R0 .EQ. 0.0D0) R0 = 1.0D0
      SRUR = WK(1)
      JMIN = IWK(IPIGP)
      DO 240 NG = 1,NGP
        JMAX = IWK(IPIGP+NG) - 1
        DO 210 J = JMIN,JMAX
          JJ = IWK(IBJGP+J)
          R = DMAX1(SRUR*DABS(Y(JJ)),R0/EWT(JJ))
 210      Y(JJ) = Y(JJ) + R
        CALL F (NEQ, TN, Y, FTEM)
        DO 230 J = JMIN,JMAX
          JJ = IWK(IBJGP+J)
          Y(JJ) = YH(JJ,1)
          R = DMAX1(SRUR*DABS(Y(JJ)),R0/EWT(JJ))
          FAC = -HL0/R
          KMIN =IWK(IBIAN+JJ) 
          KMAX =IWK(IBIAN+JJ+1) - 1
          DO 220 K = KMIN,KMAX
            I = IWK(IBJAN+K)
            WK(IBA+K) = (FTEM(I) - SAVF(I))*FAC
            IF (I .EQ. JJ) WK(IBA+K) = WK(IBA+K) + 1.0D0
 220        CONTINUE
 230      CONTINUE
        JMIN = JMAX + 1
 240    CONTINUE
      NFE = NFE + NGP
      GO TO 290
C
C IF JOK = 1, RECONSTRUCT NEW P FROM OLD P. ----------------------------
 250  JCUR = 0
      RCON = CON/CON0
      RCONT = DABS(CON)/CONMIN
      IF (RCONT .GT. RBIG .AND. IPLOST .EQ. 1) GO TO 20
      KMIN = IWK(IPIAN)
      DO 275 J = 1,N
        KMAX = IWK(IPIAN+J) - 1
        DO 270 K = KMIN,KMAX
          I = IWK(IBJAN+K)
          PIJ = WK(IBA+K)
          IF (I .NE. J) GO TO 260
          PIJ = PIJ - 1.0D0
          IF (DABS(PIJ) .GE. PSMALL) GO TO 260
            IPLOST = 1
            CONMIN = DMIN1(DABS(CON0),CONMIN)
 260      PIJ = PIJ*RCON
          IF (I .EQ. J) PIJ = PIJ + 1.0D0
          WK(IBA+K) = PIJ
 270      CONTINUE
        KMIN = KMAX + 1
 275    CONTINUE
C
C DO NUMERICAL FACTORIZATION OF P MATRIX. ------------------------------
 290  NLU = NLU + 1 
      CON0 = CON
      IERPJ = 0
      DO 295 I = 1,N
 295    FTEM(I) = 0.0D0
      CALL CDRV (N,IWK(IPR),IWK(IPC),IWK(IPIC),IWK(IPIAN),IWK(IPJAN), 
     1   WK(IPA),FTEM,FTEM,NSP,IWK(IPISP),WK(IPRSP),IESP,2,IYS)
      IF (IYS .EQ. 0) RETURN
      IMUL = (IYS - 1)/N
      IERPJ = -2
      IF (IMUL .EQ. 8) IERPJ = 1
      IF (IMUL .EQ. 10) IERPJ = -1
      RETURN
C
C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
 300  CONTINUE
      JCUR = 1
      NJE = NJE + 1 
      WK(2) = HL0
      IERPJ = 0
      R = EL0*0.1D0 
      DO 310 I = 1,N
 310    Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
      CALL F (NEQ, TN, Y, WK(3))
      NFE = NFE + 1 
      DO 320 I = 1,N
        R0 = H*SAVF(I) - YH(I,2)
        DI = 0.1D0*R0 - H*(WK(I+2) - SAVF(I))
        WK(I+2) = 1.0D0
        IF (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320
        IF (DABS(DI) .EQ. 0.0D0) GO TO 330
        WK(I+2) = 0.1D0*R0/DI 
 320    CONTINUE
      RETURN
 330  IERPJ = 2
      RETURN
C----------------------- END OF SUBROUTINE PRJS ------------------------
      END 
