added README_changes.txt

[wrffire.git] / wrfv2_fire / chem / KPP / kpp / kpp-2.1 / int / oldies / kpp_ros4.f

      SUBROUTINE INTEGRATE( TIN, TOUT )

      IMPLICIT KPP_REAL (A-H,O-Z)        
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'

C TIN - Start Time
      KPP_REAL TIN
C TOUT - End Time
      KPP_REAL TOUT
      INTEGER i

      PARAMETER (LWORK=2*NVAR*NVAR+14*NVAR+20,LIWORK=NVAR+20)
      KPP_REAL WORK(LWORK)
      INTEGER IWORK(LIWORK)
      EXTERNAL FUNC_CHEM,JAC_CHEM,SOLOUT

      ITOL=1     ! --- VECTOR TOLERANCES
      IJAC=1     ! --- COMPUTE THE JACOBIAN ANALYTICALLY
      MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
      MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
      IMAS=0     ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
      IOUT=0     ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION
      IDFX=0     ! --- INTERNAL TIME DERIVATIVE

      DO i=1,20
        IWORK(i) = 0
        WORK(i) = 0.D0
      ENDDO
      IWORK(2) = 6

      CALL KPP_ROS4(NVAR,FUNC_CHEM,Autonomous,TIN,VAR,TOUT,
     &                  STEPMIN,RTOL,ATOL,ITOL,
     &                  JAC_CHEM,IJAC,MLJAC,MUJAC,FUNC_CHEM,IDFX,
     &                  FUNC_CHEM,IMAS,MLJAC,MUJAC,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,IDID)

      IF (IDID.LT.0) THEN
        print *,'KPP_ROS4: Unsucessful step at T=',
     &          TIN,' (IDID=',IDID,')'
      ENDIF

      RETURN
      END


      SUBROUTINE KPP_ROS4(N,FCN,IFCN,X,Y,XEND,H,
     &                  RelTol,AbsTol,ITOL,
     &                  JAC ,IJAC,MLJAC,MUJAC,DFX,IDFX,
     &                  MAS ,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,IDID)
C ----------------------------------------------------------
C     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
C     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS  MY'=F(X,Y).
C     THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4  
C     (WITH STEP SIZE CONTROL).
C     C.F. SECTION IV.7
C
C     AUTHORS: E. HAIRER AND G. WANNER
C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
C              CH-1211 GENEVE 24, SWITZERLAND 
C              E-MAIL:  HAIRER@CGEUGE51.BITNET,  WANNER@CGEUGE51.BITNET
C     
C     THIS CODE IS PART OF THE BOOK:
C         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
C         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
C         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
C         SPRINGER-VERLAG (1990)               
C      
C     VERSION OF NOVEMBER 17, 1992
C
C     INPUT PARAMETERS  
C     ----------------  
C     N           DIMENSION OF THE SYSTEM 
C
C     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
C                 VALUE OF F(X,Y):
C                    SUBROUTINE FCN(N,X,Y,F)
C                    KPP_REAL X,Y(N),F(N)
C                    F(1)=...   ETC.
C
C     IFCN        GIVES INFORMATION ON FCN:
C                    IFCN=0: F(X,Y) INDEPENDENT OF X (AUTONOMOUS)
C                    IFCN=1: F(X,Y) MAY DEPEND ON X (NON-AUTONOMOUS)
C
C     X           INITIAL X-VALUE
C
C     Y(N)        INITIAL VALUES FOR Y
C
C     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
C
C     H           INITIAL STEP SIZE GUESS;
C                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, 
C                 H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD.
C                 THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
C                 ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW
C                 STEPS IN SUBROUTINE "SOLOUT", WHEN YOU ARE NOT SURE.
C                 (IF H=0.D0, THE CODE PUTS H=1.D-6).
C
C     RelTol,AbsTol   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
C                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
C
C     ITOL        SWITCH FOR RelTol AND AbsTol:
C                   ITOL=0: BOTH RelTol AND AbsTol ARE SCALARS.
C                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
C                     Y(I) BELOW RelTol*ABS(Y(I))+AbsTol
C                   ITOL=1: BOTH RelTol AND AbsTol ARE VECTORS.
C                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
C                     RelTol(I)*ABS(Y(I))+AbsTol(I).
C
C     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
C                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y
C                 (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY
C                 A DUMMY SUBROUTINE IN THE CASE IJAC=0).
C                 FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM
C                    SUBROUTINE JAC(N,X,Y,DFY,LDFY)
C                    KPP_REAL X,Y(N),DFY(LDFY,N)
C                    DFY(1,1)= ...
C                 LDFY, THE COLOMN-LENGTH OF THE ARRAY, IS
C                 FURNISHED BY THE CALLING PROGRAM.
C                 IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO
C                    BE FULL AND THE PARTIAL DERIVATIVES ARE
C                    STORED IN DFY AS
C                       DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J)
C                 ELSE, THE JACOBIAN IS TAKEN AS BANDED AND
C                    THE PARTIAL DERIVATIVES ARE STORED
C                    DIAGONAL-WISE AS
C                       DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
C
C     IJAC        SWITCH FOR THE COMPUTATION OF THE JACOBIAN:
C                    IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE
C                       DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED.
C                    IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.
C
C     MLJAC       SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN:
C                    MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR
C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
C                    0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN 
C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
C                       THE MAIN DIAGONAL).
C
C     MUJAC       UPPER BANDWITH OF JACOBIAN  MATRIX (>= NUMBER OF NON-
C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
C                 NEED NOT BE DEFINED IF MLJAC=N.
C
C     DFX         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
C                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO X
C                 (THIS ROUTINE IS ONLY CALLED IF IDFX=1 AND IFCN=1;
C                 SUPPLY A DUMMY SUBROUTINE IN THE CASE IDFX=0 OR IFCN=0).
C                 OTHERWISE, THIS SUBROUTINE MUST HAVE THE FORM
C                    SUBROUTINE DFX(N,X,Y,FX)
C                    KPP_REAL X,Y(N),FX(N)
C                    FX(1)= ...
C                
C     IDFX        SWITCH FOR THE COMPUTATION OF THE DF/DX:
C                    IDFX=0: DF/DX IS COMPUTED INTERNALLY BY FINITE
C                       DIFFERENCES, SUBROUTINE "DFX" IS NEVER CALLED.
C                    IDFX=1: DF/DX IS SUPPLIED BY SUBROUTINE DFX.
C
C     ----   MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS      -----
C     ----   FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): -
C
C     MAS         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS-
C                 MATRIX M.
C                 IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY
C                 MATRIX AND NEEDS NOT TO BE DEFINED;
C                 SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
C                 IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM
C                    SUBROUTINE MAS(N,AM,LMAS)
C                    KPP_REAL AM(LMAS,N)
C                    AM(1,1)= ....
C                    IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED
C                    AS FULL MATRIX LIKE
C                         AM(I,J) = M(I,J)
C                    ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED
C                    DIAGONAL-WISE AS
C                         AM(I-J+MUMAS+1,J) = M(I,J).
C
C     IMAS       GIVES INFORMATION ON THE MASS-MATRIX:
C                    IMAS=0: M IS SUPPOSED TO BE THE IDENTITY
C                       MATRIX, MAS IS NEVER CALLED.
C                    IMAS=1: MASS-MATRIX  IS SUPPLIED.
C
C     MLMAS       SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX:
C                    MLMAS=N: THE FULL MATRIX CASE. THE LINEAR
C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
C                    0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE
C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
C                       THE MAIN DIAGONAL).
C                 MLMAS IS SUPPOSED TO BE .LE. MLJAC.
C
C     MUMAS       UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON-
C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
C                 NEED NOT BE DEFINED IF MLMAS=N.
C                 MUMAS IS SUPPOSED TO BE .LE. MUJAC.
C
C     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
C                 NUMERICAL SOLUTION DURING INTEGRATION. 
C                 IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
C                 SUPPLY A DUMMY SUBROUTINE IF IOUT=0. 
C                 IT MUST HAVE THE FORM
C                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN)
C                    KPP_REAL X,Y(N)
C                    ....  
C                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
C                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
C                    THE FIRST GRID-POINT).
C                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
C                    IS SET <0, ROS4 RETURNS TO THE CALLING PROGRAM.
C           
C     IOUT        GIVES INFORMATION ON THE SUBROUTINE SOLOUT:
C                    IOUT=0: SUBROUTINE IS NEVER CALLED
C                    IOUT=1: SUBROUTINE IS USED FOR OUTPUT
C
C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
C                 SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES.
C                 "LWORK" MUST BE AT LEAST
C                             N*(LJAC+LMAS+LE1+8)+5
C                 WHERE
C                    LJAC=N              IF MLJAC=N (FULL JACOBIAN)
C                    LJAC=MLJAC+MUJAC+1  IF MLJAC<N (BANDED JAC.)
C                 AND                  
C                    LMAS=0              IF IMAS=0
C                    LMAS=N              IF IMAS=1 AND MLMAS=N (FULL)
C                    LMAS=MLMAS+MUMAS+1  IF MLMAS<N (BANDED MASS-M.)
C                 AND
C                    LE1=N               IF MLJAC=N (FULL JACOBIAN)
C                    LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.). 
c
C                 IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
C                 MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
C                 STORAGE REQUIREMENT IS 
C                             LWORK = 2*N*N+8*N+5.
C
C     LWORK       DECLARED LENGHT OF ARRAY "WORK".
C
C     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK".
C                 "LIWORK" MUST BE AT LEAST N+2.
C
C     LIWORK      DECLARED LENGHT OF ARRAY "IWORK".
C
C ----------------------------------------------------------------------
C 
C     SOPHISTICATED SETTING OF PARAMETERS
C     -----------------------------------
C              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK 
C              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(5)
C              AS WELL AS IWORK(1),IWORK(2) DIFFERENT FROM ZERO.
C              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
C
C    IWORK(1)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
C              THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000.
C
C    IWORK(2)  SWITCH FOR THE CHOICE OF THE COEFFICIENTS
C              IF IWORK(2).EQ.1  METHOD OF SHAMPINE
C              IF IWORK(2).EQ.2  METHOD GRK4T OF KAPS-RENTROP
C              IF IWORK(2).EQ.3  METHOD GRK4A OF KAPS-RENTROP
C              IF IWORK(2).EQ.4  METHOD OF VAN VELDHUIZEN (GAMMA=1/2)
C              IF IWORK(2).EQ.5  METHOD OF VAN VELDHUIZEN ("D-STABLE")
C              IF IWORK(2).EQ.6  AN L-STABLE METHOD
C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=2.
C
C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
C
C    WORK(2)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
C
C    WORK(3), WORK(4)   PARAMETERS FOR STEP SIZE SELECTION
C              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
C                 WORK(3) <= HNEW/HOLD <= WORK(4)
C              DEFAULT VALUES: WORK(3)=0.2D0, WORK(4)=6.D0
C
C    WORK(5)   AVOID THE HUMP: AFTER TWO CONSECUTIVE STEP REJECTIONS
C              THE STEP SIZE IS MULTIPLIED BY WORK(5)
C              DEFAULT VALUES: WORK(5)=0.1D0
C
C-----------------------------------------------------------------------
C
C     OUTPUT PARAMETERS 
C     ----------------- 
C     X           X-VALUE WHERE THE SOLUTION IS COMPUTED
C                 (AFTER SUCCESSFUL RETURN X=XEND)
C
C     Y(N)        SOLUTION AT X
C  
C     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
C
C     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
C                   IDID=1  COMPUTATION SUCCESSFUL,
C                   IDID=-1 COMPUTATION UNSUCCESSFUL.
C
C --------------------------------------------------------- 
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C          DECLARATIONS 
C *** *** *** *** *** *** *** *** *** *** *** *** ***
      IMPLICIT KPP_REAL (A-H,O-Z)
      DIMENSION Y(N),AbsTol(1),RelTol(1),WORK(LWORK),IWORK(LIWORK)
      LOGICAL AUTNMS,IMPLCT,JBAND,ARRET
      EXTERNAL FCN,JAC,DFX,MAS,SOLOUT
      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
C -----------------------------------------------------
C --- COMMON STAT CAN BE USED FOR STATISTICS
C ---    NFCN      NUMBER OF FUNC_CHEMCTION EVALUATIONS (THOSE FOR NUMERICAL
C                  EVALUATION OF THE JACOBIAN ARE NOT COUNTED)  
C ---    NJAC      NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
C                  OR NUMERICALLY)
C ---    NSTEP     NUMBER OF COMPUTED STEPS
C ---    NACCPT    NUMBER OF ACCEPTED STEPS
C ---    NREJCT    NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP
C                  HAS BEEN ACCEPTED)
C ---    NDEC      NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX)
C ---    NSOL      NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS
C -----------------------------------------------------
C *** *** *** *** *** *** ***
C        SETTING THE PARAMETERS 
C *** *** *** *** *** *** ***
      NFCN=0
      NJAC=0
      NSTEP=0
      NACCPT=0
      NREJCT=0
      NDEC=0
      NSOL=0
      ARRET=.FALSE.
C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
      IF(IWORK(1).EQ.0)THEN
         NMAX=100000
      ELSE
         NMAX=IWORK(1)
         IF(NMAX.LE.0)THEN
            WRITE(6,*)' WRONG INPUT IWORK(1)=',IWORK(1)
            ARRET=.TRUE.
         END IF
      END IF
C -------- METH   COEFFICIENTS OF THE METHOD
      IF(IWORK(2).EQ.0)THEN
         METH=2
      ELSE
         METH=IWORK(2)
         IF(METH.LE.0.OR.METH.GE.7)THEN
            WRITE(6,*)' CURIOUS INPUT IWORK(2)=',IWORK(2)
            ARRET=.TRUE.
         END IF
      END IF
C -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0  
      IF(WORK(1).EQ.0.D0)THEN
         UROUND=1.D-16
      ELSE
         UROUND=WORK(1)
         IF(UROUND.LE.1.D-14.OR.UROUND.GE.1.D0)THEN
            WRITE(6,*)' COEFFICIENTS HAVE 16 DIGITS, UROUND=',WORK(1)
            ARRET=.TRUE.
         END IF
      END IF
C -------- MAXIMAL STEP SIZE
      IF(WORK(2).EQ.0.D0)THEN
         HMAX=XEND-X
      ELSE
         HMAX=WORK(2)
      END IF
C -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
      IF(WORK(3).EQ.0.D0)THEN
         FAC1=5.D0
      ELSE
         FAC1=1.D0/WORK(3)
      END IF
      IF(WORK(4).EQ.0.D0)THEN
         FAC2=1.D0/6.0D0
      ELSE
         FAC2=1.D0/WORK(4)
      END IF
C -------  FACREJ    FOR THE HUMP
      IF(WORK(5).EQ.0.D0)THEN
         FACREJ=0.1D0
      ELSE
         FACREJ=WORK(5)
      END IF
C --------- CHECK IF TOLERANCES ARE O.K.
      IF (ITOL.EQ.0) THEN
          IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN
              WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
              ARRET=.TRUE.
          END IF
      ELSE
          DO 15 I=1,N
          IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN
              WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
              ARRET=.TRUE.
          END IF
  15      CONTINUE
      END IF
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C         COMPUTATION OF ARRAY ENTRIES
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ?
      AUTNMS=IFCN.EQ.0
      IMPLCT=IMAS.NE.0
      JBAND=MLJAC.NE.N
      ARRET=.FALSE.
C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
C -- JACOBIAN 
      IF(JBAND)THEN
         LDJAC=MLJAC+MUJAC+1
      ELSE
         LDJAC=N
      END IF
C -- MATRIX E FOR LINEAR ALGEBRA
      IF(JBAND)THEN
         LDE=2*MLJAC+MUJAC+1
      ELSE
         LDE=N
      END IF
C -- MASS MATRIX
      IF (IMPLCT) THEN
          IF (MLMAS.NE.N) THEN
              LDMAS=MLMAS+MUMAS+1
          ELSE
              LDMAS=N
          END IF
C ------ BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF "JAC"
          IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN
              WRITE (6,*) 'BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF
     & "JAC"'
            ARRET=.TRUE.
          END IF
      ELSE
          LDMAS=0
      END IF
      LDMAS2=MAX(1,LDMAS)
C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
      IEYNEW=6
      IEDY1=IEYNEW+N
      IEDY=IEDY1+N
      IEAK1=IEDY+N
      IEAK2=IEAK1+N
      IEAK3=IEAK2+N
      IEAK4=IEAK3+N
      IEFX =IEAK4+N
      IEJAC=IEFX +N
      IEMAS=IEJAC+N*LDJAC
      IEE  =IEMAS+N*LDMAS
C ------ TOTAL STORAGE REQUIREMENT -----------
      ISTORE=IEE+N*LDE-1
      IF(ISTORE.GT.LWORK)THEN
         WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
         ARRET=.TRUE.
      END IF
C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
      IEIP=3
C --------- TOTAL REQUIREMENT ---------------
      ISTORE=IEIP+N-1
      IF(ISTORE.GT.LIWORK)THEN
         WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
         ARRET=.TRUE.
      END IF
C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
      IF (ARRET) THEN
         IDID=-1
         RETURN
      END IF
C -------- CALL TO CORE INTEGRATOR ------------
      CALL RO4COR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,JAC,IJAC,
     &   MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
     &   NMAX,UROUND,METH,FAC1,FAC2,FACREJ,AUTNMS,IMPLCT,JBAND,
     &   LDJAC,LDE,LDMAS2,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY),
     &   WORK(IEAK1),WORK(IEAK2),WORK(IEAK3),WORK(IEAK4),
     &   WORK(IEFX),WORK(IEJAC),WORK(IEE),WORK(IEMAS),IWORK(IEIP))
C ----------- RETURN -----------
      RETURN
      END
C
C
C
C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
C
      SUBROUTINE RO4COR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,JAC,
     &  IJAC,MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
     &  NMAX,UROUND,METH,FAC1,FAC2,FACREJ,AUTNMS,IMPLCT,BANDED,
     &  LFJAC,LE,LDMAS,YNEW,DY1,DY,AK1,AK2,AK3,AK4,FX,FJAC,E,FMAS,IP)
C ----------------------------------------------------------
C     CORE INTEGRATOR FOR ROS4
C     PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED 
C ---------------------------------------------------------- 
C         DECLARATIONS 
C ---------------------------------------------------------- 
      IMPLICIT KPP_REAL (A-H,O-Z)        
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_sparse.h'
      KPP_REAL Y(N),YNEW(N),DY1(N),DY(N),AK1(N),
     *  AK2(N),AK3(N),AK4(N),FX(N),
     *  FJAC(LU_NONZERO),E(LU_NONZERO),
     *  FMAS(LDMAS,N),AbsTol(1),RelTol(1)
      INTEGER IP(N)
      LOGICAL REJECT,RJECT2,AUTNMS,IMPLCT,BANDED
      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
      EXTERNAL FCN,JAC,MAS,SOLOUT,DFX


C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
      IF (IMPLCT) CALL MAS(N,FMAS,LDMAS)
C ---- PREPARE BANDWIDTHS -----
      IF (BANDED) THEN
          MLE=MLJAC
          MUE=MUJAC
          MBJAC=MLJAC+MUJAC+1
          MBB=MLMAS+MUMAS+1
          MDIAG=MLE+MUE+1
          MBDIAG=MUMAS+1
          MDIFF=MLE+MUE-MUMAS
      END IF
C *** *** *** *** *** *** ***
C  INITIALISATIONS
C *** *** *** *** *** *** ***
      IF (METH.EQ.1) CALL SHAMP (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      IF (METH.EQ.2) CALL GRK4T (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      IF (METH.EQ.3) CALL GRK4A (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      IF (METH.EQ.4) CALL VELDS (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      IF (METH.EQ.5) CALL VELDD (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      IF (METH.EQ.6) CALL LSTAB (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
C --- INITIAL PREPARATIONS
      POSNEG=SIGN(1.D0,XEND-X)
      HMAXN=MIN(ABS(HMAX),ABS(XEND-X))
      IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
      H=MIN(ABS(H),HMAXN) 
      H=SIGN(H,POSNEG) 
      REJECT=.FALSE.
      NSING=0
      IRTRN=1 
      XXOLD=X
      IF (IRTRN.LT.0) GOTO 79
C --- BASIC INTEGRATION STEP  
   1  IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79
      IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN
          H=HOPT
          IDID=1
          RETURN
      END IF
      HOPT=H
      IF ((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X
      CALL FCN(N,X,Y,DY1)
      NFCN=NFCN+1

C *** *** *** *** *** *** ***
C  COMPUTATION OF THE JACOBIAN
C *** *** *** *** *** *** ***
      NJAC=NJAC+1
      CALL JAC(N,X,Y,FJAC)

      IF (.NOT.AUTNMS) THEN
          IF (IDFX.EQ.0) THEN
C --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X
              DELT=DSQRT(UROUND*MAX(1.D-5,ABS(X)))
              XDELT=X+DELT
              CALL FCN(N,XDELT,Y,AK1)
              DO 19 J=1,N
  19          FX(J)=(AK1(J)-DY1(J))/DELT
          ELSE
C --- COMPUTE ANALYTICALLY THE DERIVATIVE WITH RESPECT TO X
              CALL DFX(N,X,Y,FX)
          END IF
      END IF
   2  CONTINUE

C *** *** *** *** *** *** ***
C  COMPUTE THE STAGES
C *** *** *** *** *** *** ***
      NDEC=NDEC+1
      HC21=C21/H
      HC31=C31/H
      HC32=C32/H
      HC41=C41/H
      HC42=C42/H
      HC43=C43/H
      FAC=1.D0/(H*GAMMA)
C --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX)
              DO 800 I=1,LU_NONZERO
  800         E(I)=-FJAC(I)
              DO 801 J=1,N
  801         E(LU_DIAG(J))=E(LU_DIAG(J))+FAC
              CALL KppDecomp (E,IER)
              IF (IER.NE.0) GOTO 80
              IF (AUTNMS) THEN
C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
C ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
C ---   2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX
C ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS
                  DO 803 I=1,N
  803             AK1(I)=DY1(I)
                  CALL KppSolve(E,AK1)
                  DO 810 I=1,N
  810             YNEW(I)=Y(I)+A21*AK1(I) 
                  CALL FCN(N,X,YNEW,DY)
                  DO 811 I=1,N
  811             AK2(I)=DY(I)+HC21*AK1(I)
                  CALL KppSolve(E,AK2)
                  DO 820 I=1,N
  820             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                  CALL FCN(N,X,YNEW,DY)
                  DO 821 I=1,N
  821             AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I)
                  CALL KppSolve(E,AK3)
                  DO 831 I=1,N
  831             AK4(I)=DY(I)+HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 
                  CALL KppSolve(E,AK4)
              ELSE
C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
C ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
C ---   2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX
C ---   3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS
                  HD1=H*D1
                  HD2=H*D2
                  HD3=H*D3
                  HD4=H*D4
                  DO 903 I=1,N
  903             AK1(I)=DY1(I)+HD1*FX(I)
                  CALL KppSolve(E,AK1)
                  DO 910 I=1,N
  910             YNEW(I)=Y(I)+A21*AK1(I) 
                  CALL FCN(N,X+C2*H,YNEW,DY)
                  DO 911 I=1,N
  911             AK2(I)=DY(I)+HD2*FX(I)+HC21*AK1(I)
                  CALL KppSolve(E,AK2)
                  DO 920 I=1,N
  920             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                  CALL FCN(N,X+C3*H,YNEW,DY)
                  DO 921 I=1,N
  921             AK3(I)=DY(I)+HD3*FX(I)+HC31*AK1(I)+HC32*AK2(I)
                  CALL KppSolve(E,AK3)
                  DO 931 I=1,N
  931             AK4(I)=DY(I)+HD4*FX(I)+HC41*AK1(I)+HC42*AK2(I)
     &                  +HC43*AK3(I) 
                  CALL KppSolve(E,AK4)
              END IF
      NSOL=NSOL+4
      NFCN=NFCN+2
C *** *** *** *** *** *** ***
C  ERROR ESTIMATION  
C *** *** *** *** *** *** ***
      NSTEP=NSTEP+1
C ------------ NEW SOLUTION ---------------
      DO 240 I=1,N
  240 YNEW(I)=Y(I)+B1*AK1(I)+B2*AK2(I)+B3*AK3(I)+B4*AK4(I)  
C ------------ COMPUTE ERROR ESTIMATION ----------------
      ERR=0.D0
      DO 300 I=1,N
      S=E1*AK1(I)+E2*AK2(I)+E3*AK3(I)+E4*AK4(I) 
      IF (ITOL.EQ.0) THEN
         SK=AbsTol(1)+RelTol(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
      ELSE
         SK=AbsTol(I)+RelTol(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
      END IF
  300 ERR=ERR+(S/SK)**2
      ERR=DSQRT(ERR/N)
C --- COMPUTATION OF HNEW
C --- WE REQUIRE .2<=HNEW/H<=6.
      FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.25D0/.9D0))
      HNEW=H/FAC  
C *** *** *** *** *** *** ***
C  IS THE ERROR SMALL ENOUGH ?
C *** *** *** *** *** *** ***
      IF (ERR.LE.1.D0) THEN
C --- STEP IS ACCEPTED  
         NACCPT=NACCPT+1
         DO 44 I=1,N 
  44     Y(I)=YNEW(I) 
         XXOLD=X 
         X=X+H
         IF (IRTRN.LT.0) GOTO 79
         IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN
         IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H)) 
         REJECT=.FALSE.
         RJECT2=.FALSE.
         H=HNEW
         GOTO 1
      ELSE
C --- STEP IS REJECTED  
         IF (RJECT2) HNEW=H*FACREJ
         IF (REJECT) RJECT2=.TRUE.
         REJECT=.TRUE.
         H=HNEW
         IF (NACCPT.GE.1) NREJCT=NREJCT+1
         GOTO 2
      END IF
C --- EXIT
  80  WRITE (6,*) ' MATRIX E IS SINGULAR, INFO = ',INFO
      NSING=NSING+1
      IF (NSING.GE.5) GOTO 79
      H=H*0.5D0
      GOTO 2
  79  WRITE(6,979)X,H
 979  FORMAT(' EXIT OF ROS4 AT X=',D16.7,'   H=',D16.7)
      IDID=-1
      RETURN
      END
C
      SUBROUTINE SHAMP (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      IMPLICIT KPP_REAL (A-H,O-Z)
         A21=2.D0
         A31=48.D0/25.D0
         A32=6.D0/25.D0
         C21=-8.D0
         C31=372.D0/25.D0
         C32=12.D0/5.D0
         C41=-112.D0/125.D0
         C42=-54.D0/125.D0
         C43=-2.D0/5.D0
         B1=19.D0/9.D0
         B2=1.D0/2.D0
         B3=25.D0/108.D0
         B4=125.D0/108.D0
         E1=17.D0/54.D0
         E2=7.D0/36.D0
         E3=0.D0
         E4=125.D0/108.D0
         GAMMA=.5D0
         C2= 0.1000000000000000D+01
         C3= 0.6000000000000000D+00
         D1= 0.5000000000000000D+00
         D2=-0.1500000000000000D+01
         D3= 0.2420000000000000D+01
         D4= 0.1160000000000000D+00
      RETURN
      END
C
      SUBROUTINE GRK4A (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      IMPLICIT KPP_REAL (A-H,O-Z)
       A21= 0.1108860759493671D+01
       A31= 0.2377085261983360D+01
       A32= 0.1850114988899692D+00
       C21=-0.4920188402397641D+01
       C31= 0.1055588686048583D+01
       C32= 0.3351817267668938D+01
       C41= 0.3846869007049313D+01
       C42= 0.3427109241268180D+01
       C43=-0.2162408848753263D+01
       B1= 0.1845683240405840D+01
       B2= 0.1369796894360503D+00
       B3= 0.7129097783291559D+00
       B4= 0.6329113924050632D+00
       E1= 0.4831870177201765D-01
       E2=-0.6471108651049505D+00
       E3= 0.2186876660500240D+00
       E4=-0.6329113924050632D+00
       GAMMA= 0.3950000000000000D+00
       C2= 0.4380000000000000D+00
       C3= 0.8700000000000000D+00
       D1= 0.3950000000000000D+00
       D2=-0.3726723954840920D+00
       D3= 0.6629196544571492D-01
       D4= 0.4340946962568634D+00
      RETURN
      END
C
      SUBROUTINE GRK4T (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
      IMPLICIT KPP_REAL (A-H,O-Z)
       A21= 0.2000000000000000D+01
       A31= 0.4524708207373116D+01
       A32= 0.4163528788597648D+01
       C21=-0.5071675338776316D+01
       C31= 0.6020152728650786D+01
       C32= 0.1597506846727117D+00
       C41=-0.1856343618686113D+01
       C42=-0.8505380858179826D+01
       C43=-0.2084075136023187D+01
       B1= 0.3957503746640777D+01
       B2= 0.4624892388363313D+01
       B3= 0.6174772638750108D+00
       B4= 0.1282612945269037D+01
       E1= 0.2302155402932996D+01
       E2= 0.3073634485392623D+01
       E3=-0.8732808018045032D+00
       E4=-0.1282612945269037D+01
       GAMMA= 0.2310000000000000D+00
       C2= 0.4620000000000000D+00
       C3= 0.8802083333333334D+00
       D1= 0.2310000000000000D+00
       D2=-0.3962966775244303D-01
       D3= 0.5507789395789127D+00
       D4=-0.5535098457052764D-01
      RETURN
      END
C
      SUBROUTINE VELDS (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
C --- METHOD GIVEN BY VAN VELDHUIZEN
      IMPLICIT KPP_REAL (A-H,O-Z)
       A21= 0.2000000000000000D+01
       A31= 0.1750000000000000D+01
       A32= 0.2500000000000000D+00
       C21=-0.8000000000000000D+01
       C31=-0.8000000000000000D+01
       C32=-0.1000000000000000D+01
       C41= 0.5000000000000000D+00
       C42=-0.5000000000000000D+00
       C43= 0.2000000000000000D+01
       B1= 0.1333333333333333D+01
       B2= 0.6666666666666667D+00
       B3=-0.1333333333333333D+01
       B4= 0.1333333333333333D+01
       E1=-0.3333333333333333D+00
       E2=-0.3333333333333333D+00
       E3=-0.0000000000000000D+00
       E4=-0.1333333333333333D+01
       GAMMA= 0.5000000000000000D+00
       C2= 0.1000000000000000D+01
       C3= 0.5000000000000000D+00
       D1= 0.5000000000000000D+00
       D2=-0.1500000000000000D+01
       D3=-0.7500000000000000D+00
       D4= 0.2500000000000000D+00
      RETURN
      END
C
      SUBROUTINE VELDD (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
C --- METHOD GIVEN BY VAN VELDHUIZEN
      IMPLICIT KPP_REAL (A-H,O-Z)
       A21= 0.2000000000000000D+01
       A31= 0.4812234362695436D+01
       A32= 0.4578146956747842D+01
       C21=-0.5333333333333331D+01
       C31= 0.6100529678848254D+01
       C32= 0.1804736797378427D+01
       C41=-0.2540515456634749D+01
       C42=-0.9443746328915205D+01
       C43=-0.1988471753215993D+01
       B1= 0.4289339254654537D+01
       B2= 0.5036098482851414D+01
       B3= 0.6085736420673917D+00
       B4= 0.1355958941201148D+01
       E1= 0.2175672787531755D+01
       E2= 0.2950911222575741D+01
       E3=-0.7859744544887430D+00
       E4=-0.1355958941201148D+01
       GAMMA= 0.2257081148225682D+00
       C2= 0.4514162296451364D+00
       C3= 0.8755928946018455D+00
       D1= 0.2257081148225682D+00
       D2=-0.4599403502680582D-01
       D3= 0.5177590504944076D+00
       D4=-0.3805623938054428D-01
      RETURN
      END 
C
      SUBROUTINE LSTAB (A21,A31,A32,C21,C31,C32,C41,C42,C43,
     &          B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4)
C --- AN L-STABLE METHOD
      IMPLICIT KPP_REAL (A-H,O-Z)
       A21= 0.2000000000000000D+01
       A31= 0.1867943637803922D+01
       A32= 0.2344449711399156D+00
       C21=-0.7137615036412310D+01
       C31= 0.2580708087951457D+01
       C32= 0.6515950076447975D+00
       C41=-0.2137148994382534D+01
       C42=-0.3214669691237626D+00
       C43=-0.6949742501781779D+00
C  B_i = M_i
       B1= 0.2255570073418735D+01
       B2= 0.2870493262186792D+00
       B3= 0.4353179431840180D+00
       B4= 0.1093502252409163D+01
C E_i = error estimator       
       E1=-0.2815431932141155D+00
       E2=-0.7276199124938920D-01
       E3=-0.1082196201495311D+00
       E4=-0.1093502252409163D+01
C gamma = gamma
       GAMMA= 0.5728200000000000D+00
C  C_i = alpha_i
       C2= 0.1145640000000000D+01
       C3= 0.6552168638155900D+00
C D_i = gamma_i       
       D1= 0.5728200000000000D+00
       D2=-0.1769193891319233D+01
       D3= 0.7592633437920482D+00
       D4=-0.1049021087100450D+00
      RETURN
      END 

        SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN)
C --- PRINTS SOLUTION
        IMPLICIT KPP_REAL (A-H,O-Z)
        DIMENSION Y(N)
        COMMON /INTERN/XOUT
        IF (NR.EQ.1) THEN
           WRITE (6,99) X,Y(1),Y(2),NR-1
           XOUT=0.1D0
        ELSE
           IF (X.GE.XOUT) THEN
              WRITE (6,99) X,Y(1),Y(2),NR-1
              XOUT=DMAX1(XOUT+0.1D0,X)
           END IF
        END IF
 99     FORMAT(1X,'X =',F5.2,'    Y =',2E18.10,'    NSTEP =',I4)
        RETURN
        END

      SUBROUTINE DEC (N, NDIM, A, IP, IER)
C VERSION REAL KPP_REAL
      INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J
      KPP_REAL A,T
      DIMENSION A(NDIM,N), IP(N)
C-----------------------------------------------------------------------
C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
C  INPUT..
C     N = ORDER OF MATRIX.
C     NDIM = DECLARED DIMENSION OF ARRAY  A .
C     A = MATRIX TO BE TRIANGULARIZED.
C  OUTPUT..
C     A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
C     A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
C           SINGULAR AT STAGE K.
C  USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
C  DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
C
C  REFERENCE..
C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
C     C.A.C.M. 15 (1972), P. 274.
C-----------------------------------------------------------------------
      IER = 0
      IP(N) = 1
      IF (N .EQ. 1) GO TO 70
      NM1 = N - 1
      DO 60 K = 1,NM1
        KP1 = K + 1
        M = K
        DO 10 I = KP1,N
          IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I  
 10     CONTINUE
        IP(K) = M
        T = A(M,K)
        IF (M .EQ. K) GO TO 20
        IP(N) = -IP(N)
        A(M,K) = A(K,K)
        A(K,K) = T
 20     CONTINUE
        IF (T .EQ. 0.D0) GO TO 80
        T = 1.D0/T
        DO 30 I = KP1,N
 30       A(I,K) = -A(I,K)*T
        DO 50 J = KP1,N
          T = A(M,J)
          A(M,J) = A(K,J)
          A(K,J) = T
          IF (T .EQ. 0.D0) GO TO 45
          DO 40 I = KP1,N
 40         A(I,J) = A(I,J) + A(I,K)*T
 45       CONTINUE
 50       CONTINUE
 60     CONTINUE
 70   K = N
      IF (A(N,N) .EQ. 0.D0) GO TO 80
      RETURN
 80   IER = K
      IP(N) = 0
      RETURN
C----------------------- END OF SUBROUTINE DEC -------------------------
      END
C
C
      SUBROUTINE SOL (N, NDIM, A, B, IP)
C VERSION REAL KPP_REAL
      INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1
      KPP_REAL A,B,T
      DIMENSION A(NDIM,N), B(N), IP(N)
C-----------------------------------------------------------------------
C  SOLUTION OF LINEAR SYSTEM, A*X = B .
C  INPUT..
C    N = ORDER OF MATRIX.
C    NDIM = DECLARED DIMENSION OF ARRAY  A .
C    A = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
C    B = RIGHT HAND SIDE VECTOR.
C    IP = PIVOT VECTOR OBTAINED FROM DEC.
C  DO NOT USE IF DEC HAS SET IER .NE. 0.
C  OUTPUT..
C    B = SOLUTION VECTOR, X .
C-----------------------------------------------------------------------
      IF (N .EQ. 1) GO TO 50
      NM1 = N - 1
      DO 20 K = 1,NM1
        KP1 = K + 1
        M = IP(K)
        T = B(M)
        B(M) = B(K)
        B(K) = T
        DO 10 I = KP1,N
 10       B(I) = B(I) + A(I,K)*T
 20     CONTINUE
      DO 40 KB = 1,NM1
        KM1 = N - KB
        K = KM1 + 1
        B(K) = B(K)/A(K,K)
        T = -B(K)
        DO 30 I = 1,KM1
 30       B(I) = B(I) + A(I,K)*T
 40     CONTINUE
 50   B(1) = B(1)/A(1,1)
      RETURN
C----------------------- END OF SUBROUTINE SOL -------------------------
      END


 
 
      SUBROUTINE FUNC_CHEM(N, T, Y, P)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'
      INTEGER N
      KPP_REAL   T, Told
      KPP_REAL   Y(NVAR), P(NVAR)
      Told = TIME
      TIME = T
      CALL Update_SUN()
      CALL Update_RCONST()
      CALL Fun( Y,  FIX, RCONST, P )
      TIME = Told
      RETURN
      END

 
      SUBROUTINE JAC_CHEM(N, T, Y, J)
      INCLUDE 'KPP_ROOT_params.h'
      INCLUDE 'KPP_ROOT_global.h'
      INTEGER N
      KPP_REAL   Told, T
      KPP_REAL   Y(NVAR), J(LU_NONZERO)
      Told = TIME
      TIME = T
      CALL Update_SUN()
      CALL Update_RCONST()
      CALL Jac_SP( Y,  FIX, RCONST, J )
      TIME = Told
      RETURN
      END