chem/KPP/kpp/kpp-2.1/int/oldies/kpp_rodas.f

   1       SUBROUTINE INTEGRATE( TIN, TOUT )
   2
   3       IMPLICIT KPP_REAL (A-H,O-Z)
   4       INCLUDE 'KPP_ROOT_params.h'
   5       INCLUDE 'KPP_ROOT_global.h'
   6
   7 C TIN - Start Time
   8       KPP_REAL TIN
   9 C TOUT - End Time
  10       KPP_REAL TOUT
  11       INTEGER i
  12
  13       PARAMETER (LWORK=2*NVAR*NVAR+14*NVAR+20,LIWORK=NVAR+20)
  14       KPP_REAL WORK(LWORK)
  15       INTEGER IWORK(LIWORK)
  16       EXTERNAL FUNC_CHEM,JAC_CHEM
  17
  18
  19       ITOL=1     ! --- VECTOR TOLERANCES
  20       IJAC=1     ! --- COMPUTE THE JACOBIAN ANALYTICALLY
  21       MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
  22       MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
  23       IMAS=0     ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
  24       IOUT=0     ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION
  25       IDFX=0     ! --- INTERNAL TIME DERIVATIVE
  26
  27       DO i=1,20
  28         IWORK(i) = 0
  29         WORK(i) = 0.D0
  30       ENDDO
  31
  32       IWORK(3) = 1
  33
  34       CALL ATMRODAS(NVAR,FUNC_CHEM,Autonomous,TIN,VAR,TOUT,
  35      &                  STEPMIN,RTOL,ATOL,ITOL,
  36      &                  JAC_CHEM,IJAC,MLJAC,MUJAC,FUNC_CHEM,IDFX,
  37      &                  FUNC_CHEM,IMAS,
  38      &                  WORK,LWORK,IWORK,LIWORK,IDID)
  39
  40       IF (IDID.LT.0) THEN
  41         print *,'ATMRODAS: Unsucessfull exit at T=',
  42      &          TIN,' (IDID=',IDID,')'
  43       ENDIF
  44
  45
  46       RETURN
  47       END
  48
  49
  50       SUBROUTINE ATMRODAS(N,FCN,IFCN,X,Y,XEND,H,
  51      &                  RelTol,AbsTol,ITOL,
  52      &                  JAC ,IJAC,MLJAC,MUJAC,DFX,IDFX,
  53      &                  MAS ,IMAS,
  54      &                  WORK,LWORK,IWORK,LIWORK,IDID)
  55 C ----------------------------------------------------------
  56 C     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
  57 C     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS  MY'=F(X,Y).
  58 C     THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4
  59 C     (WITH STEP SIZE CONTROL).
  60 C     C.F. SECTIONS IV.7  AND VI.3
  61 C
  62 C     AUTHORS: E. HAIRER AND G. WANNER
  63 C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
  64 C              CH-1211 GENEVE 24, SWITZERLAND
  65 C              E-MAIL:  HAIRER@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
  66 C ---------------------------------------------------------
  67 C *** *** *** *** *** *** *** *** *** *** *** *** ***
  68 C          DECLARATIONS
  69 C *** *** *** *** *** *** *** *** *** *** *** *** ***
  70       IMPLICIT KPP_REAL (A-H,O-Z)
  71       DIMENSION Y(N),AbsTol(*),RelTol(*),WORK(LWORK),IWORK(LIWORK)
  72       LOGICAL AUTNMS,IMPLCT,JBAND,ARRET,PRED
  73       EXTERNAL FCN,JAC,DFX,MAS
  74       COMMON /STATISTICS/ NFCN,NACCPT,NREJCT,NSTEP,NJAC,NDEC,NSOL
  75 C *** *** *** *** *** *** ***
  76 C        SETTING THE PARAMETERS
  77 C *** *** *** *** *** *** ***
  78       ARRET=.FALSE.
  79       METH = 1
  80       NMAX=100000
  81 C -------- PRED   STEP SIZE CONTROL
  82       IF(IWORK(3).LE.1)THEN
  83          PRED=.TRUE.
  84       ELSE
  85          PRED=.FALSE.
  86       END IF
  87       UROUND=1.D-16
  88       NM1 = N
  89       M1  = N
  90       M2  = N
  91 C -------- MAXIMAL STEP SIZE
  92       IF(WORK(2).EQ.0.D0)THEN
  93          HMAX=XEND-X
  94       ELSE
  95          HMAX=WORK(2)
  96       END IF
  97 C -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
  98       IF(WORK(3).EQ.0.D0)THEN
  99          FAC1=5.D0
 100       ELSE
 101          FAC1=1.D0/WORK(3)
 102       END IF
 103       IF(WORK(4).EQ.0.D0)THEN
 104          FAC2=1.D0/6.0D0
 105       ELSE
 106          FAC2=1.D0/WORK(4)
 107       END IF
 108       IF (FAC1.LT.1.0D0.OR.FAC2.GT.1.0D0) THEN
 109             WRITE(6,*)' CURIOUS INPUT WORK(3,4)=',WORK(3),WORK(4)
 110             ARRET=.TRUE.
 111          END IF
 112 C --------- SAFE     SAFETY FACTOR IN STEP SIZE PREDICTION
 113       SAFE=0.9D0
 114 C --------- CHECK IF TOLERANCES ARE O.K.
 115       IF (ITOL.EQ.0) THEN
 116           IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN
 117               WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
 118               ARRET=.TRUE.
 119           END IF
 120       ELSE
 121           DO I=1,N
 122              IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN
 123                 WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
 124                 ARRET=.TRUE.
 125              END IF
 126           END DO
 127       END IF
 128
 129       IF (ARRET) STOP
 130       NM1 = N
 131 C *** *** *** *** *** *** *** *** *** *** *** *** ***
 132 C         COMPUTATION OF ARRAY ENTRIES
 133 C *** *** *** *** *** *** *** *** *** *** *** *** ***
 134 C ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ?
 135       AUTNMS=IFCN.EQ.0
 136       IMPLCT=IMAS.NE.0
 137       JBAND=MLJAC.LT.NM1
 138 C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
 139 C -- JACOBIAN AND MATRIX E
 140          MLJAC=NM1
 141          MUJAC=NM1
 142          LDJAC=NM1
 143          LDE=NM1
 144 C -- MASS MATRIX
 145       IF (IMPLCT) THEN
 146          print *, 'Implicit 1'
 147       ELSE
 148          LDMAS=0
 149          IJOB=1
 150       END IF
 151       LDMAS2=MAX(1,LDMAS)
 152 C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
 153       IEYNEW=21
 154       IEDY1=IEYNEW+N
 155       IEDY=IEDY1+N
 156       IEAK1=IEDY+N
 157       IEAK2=IEAK1+N
 158       IEAK3=IEAK2+N
 159       IEAK4=IEAK3+N
 160       IEAK5=IEAK4+N
 161       IEAK6=IEAK5+N
 162       IEFX =IEAK6+N
 163       IECON=IEFX+N
 164       IEJAC=IECON+4*N
 165       IEMAS=IEJAC+N*LDJAC
 166       IEE  =IEMAS+NM1*LDMAS
 167 C ------ TOTAL STORAGE REQUIREMENT -----------
 168       ISTORE=IEE+NM1*LDE-1
 169       IF(ISTORE.GT.LWORK)THEN
 170          WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
 171          ARRET=.TRUE.
 172       END IF
 173 C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
 174       IEIP=21
 175       ISTORE=IEIP+NM1-1
 176       IF(ISTORE.GT.LIWORK)THEN
 177          WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
 178          ARRET=.TRUE.
 179       END IF
 180 C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
 181       IF (ARRET) THEN
 182          IDID=-1
 183          RETURN
 184       END IF
 185 C -------- CALL TO CORE INTEGRATOR ------------
 186       CALL ROSCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,JAC,IJAC,
 187      &   MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,IOUT,IDID,NMAX,
 188      &   UROUND,METH,IJOB,FAC1,FAC2,SAFE,AUTNMS,IMPLCT,JBAND,PRED,LDJAC,
 189      &   LDE,LDMAS2,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY),WORK(IEAK1),
 190      &   WORK(IEAK2),WORK(IEAK3),WORK(IEAK4),WORK(IEAK5),WORK(IEAK6),
 191      &   WORK(IEFX),WORK(IEJAC),WORK(IEE),WORK(IEMAS),IWORK(IEIP),
 192      &   WORK(IECON),
 193      &   M1,M2,NM1,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL)
 194       IWORK(14)=NFCN
 195       IWORK(15)=NJAC
 196       IWORK(16)=NSTEP
 197       IWORK(17)=NACCPT
 198       IWORK(18)=NREJCT
 199       IWORK(19)=NDEC
 200       IWORK(20)=NSOL
 201 C ----------- RETURN -----------
 202       RETURN
 203       END
 204 C
 205 C
 206 C
 207 C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
 208 C
 209       SUBROUTINE ROSCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,
 210      &  ITOL,JAC,IJAC,
 211      &  MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,IOUT,IDID,NMAX,
 212      &  UROUND,METH,IJOB,FAC1,FAC2,SAFE,AUTNMS,IMPLCT,BANDED,
 213      &  PRED,LDJAC,
 214      &  LDE,LDMAS,YNEW,DY1,DY,AK1,AK2,AK3,AK4,AK5,AK6,
 215      &  FX,FJAC,E,FMAS,IP,CONT,
 216      &  M1,M2,NM1,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL)
 217 C ----------------------------------------------------------
 218 C     CORE INTEGRATOR FOR RODAS
 219 C     PARAMETERS SAME AS IN RODAS WITH WORKSPACE ADDED
 220 C ----------------------------------------------------------
 221 C         DECLARATIONS
 222 C ----------------------------------------------------------
 223       IMPLICIT KPP_REAL (A-H,O-Z)
 224       INCLUDE 'KPP_ROOT_params.h'
 225       INCLUDE 'KPP_ROOT_sparse.h'
 226       DIMENSION Y(N),YNEW(N),DY1(N),DY(N),AK1(N),
 227      *  AK2(N),AK3(N),AK4(N),AK5(N),AK6(N),FX(N),
 228      *  FJAC(LU_NONZERO),E(LDE,NM1),FMAS(LDMAS,NM1),
 229      *  AbsTol(*),RelTol(*)
 230       DIMENSION CONT(4*N)
 231       INTEGER IP(NM1)
 232       LOGICAL REJECT,AUTNMS,IMPLCT,BANDED
 233       LOGICAL ONE,LAST,PRED,SINGULAR
 234       EXTERNAL FCN, MAS, JAC, DFX
 235       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
 236       COMMON /CONROS/XOLD,HOUT,NN
 237 C *** *** *** *** *** *** ***
 238 C  INITIALISATIONS
 239 C *** *** *** *** *** *** ***
 240       NN=N
 241       NN2=2*N
 242       NN3=3*N
 243       LRC=4*N
 244 C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
 245       IF (IMPLCT) CALL MAS (NM1,FMAS,LDMAS)
 246 C ------ SET THE PARAMETERS OF THE METHOD -----
 247       CALL ROCOE(METH,A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,
 248      &    C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,
 249      &    C62,C63,C64,C65,GAMMA,C2,C3,C4,D1,D2,D3,D4,
 250      &    D21,D22,D23,D24,D25,D31,D32,D33,D34,D35)
 251 C --- INITIAL PREPARATIONS
 252       IF (M1.GT.0) IJOB=IJOB+10
 253       POSNEG=SIGN(1.D0,XEND-X)
 254       HMAXN=DMIN1(DABS(HMAX),DABS(XEND-X))
 255       IF (DABS(H).LE.10.D0*UROUND) H=1.0D-6
 256       H=DMIN1(DABS(H),HMAXN)
 257       H=SIGN(H,POSNEG)
 258       HACC = H
 259       ERRACC = 1.0d0
 260       REJECT=.FALSE.
 261       LAST=.FALSE.
 262       NSING=0
 263       IRTRN=1
 264       IF (AUTNMS) THEN
 265          HD1=0.0D0
 266          HD2=0.0D0
 267          HD3=0.0D0
 268          HD4=0.0D0
 269       END IF
 270 C -------- PREPARE BAND-WIDTHS --------
 271       MBDIAG=MUMAS+1
 272
 273 C --- BASIC INTEGRATION STEP
 274       LAST = .FALSE.
 275       DO WHILE (.NOT.LAST)
 276       IF (.NOT. REJECT) THEN
 277       IF (NSTEP.GT.NMAX) CALL FAIL_EXIT(3,X,IDID,H,NMAX)
 278       IF ( 0.1D0*DABS(H) .LE. DABS(X)*UROUND )
 279      *       CALL FAIL_EXIT(2,X,IDID,H,NMAX)
 280       HOPT=H
 281       IF ((X+H*1.0001D0-XEND)*POSNEG.GE.0.D0) THEN
 282          H=XEND-X
 283          LAST=.TRUE.
 284       END IF
 285 C *** *** *** *** *** *** ***
 286 C  COMPUTATION OF THE JACOBIAN
 287 C *** *** *** *** *** *** ***
 288       CALL FCN(N,X,Y,DY1)
 289       CALL JAC(N,X,Y,FJAC,LDJAC)
 290       NFCN=NFCN+1
 291       NJAC=NJAC+1
 292
 293       IF (.NOT.AUTNMS) THEN
 294 C --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X
 295             DELT=DSQRT(UROUND*DMAX1(1.D-5,DABS(X)))
 296             XDELT=X+DELT
 297             CALL FCN(N,XDELT,Y,AK1)
 298             DO J=1,N
 299                FX(J)=(AK1(J)-DY1(J))/DELT
 300             END DO
 301       END IF
 302       END IF
 303
 304 C *** *** *** *** *** *** ***
 305 C  COMPUTE THE STAGES
 306 C *** *** *** *** *** *** ***
 307       SINGULAR = .TRUE.
 308       DO WHILE (SINGULAR)
 309         FAC=1.D0/(H*GAMMA)
 310         CALL DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
 311      &            M1,M2,NM1,FAC,E,LDE,IP,IER,IJOB,IMPLCT,IP)
 312         SINGULAR = IER.NE.0
 313         IF (SINGULAR) THEN
 314           NSING=NSING+1
 315           IF (NSING.GE.5) CALL FAIL_EXIT(1,X,IDID,H,NMAX)
 316           H=H*0.5D0
 317           REJECT=.TRUE.
 318           LAST=.FALSE.
 319           ONE = .FALSE.
 320         END IF
 321       END DO
 322
 323       NDEC=NDEC+1
 324 C --- PREPARE FOR THE COMPUTATION OF THE 6 STAGES
 325       HC21=C21/H
 326       HC31=C31/H
 327       HC32=C32/H
 328       HC41=C41/H
 329       HC42=C42/H
 330       HC43=C43/H
 331       HC51=C51/H
 332       HC52=C52/H
 333       HC53=C53/H
 334       HC54=C54/H
 335       HC61=C61/H
 336       HC62=C62/H
 337       HC63=C63/H
 338       HC64=C64/H
 339       HC65=C65/H
 340       IF (.NOT.AUTNMS) THEN
 341          HD1=H*D1
 342          HD2=H*D2
 343          HD3=H*D3
 344          HD4=H*D4
 345       END IF
 346 C --- THE STAGES
 347       CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
 348      &    M1,M2,NM1,FAC,E,LDE,IP,DY1,AK1,FX,YNEW,HD1,IJOB,.FALSE.)
 349       DO  I=1,N
 350          YNEW(I)=Y(I)+A21*AK1(I)
 351       END DO
 352       CALL FCN(N,X+C2*H,YNEW,DY)
 353       DO I=1,N
 354          YNEW(I)=HC21*AK1(I)
 355       END DO
 356       CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
 357      &    M1,M2,NM1,FAC,E,LDE,IP,DY,AK2,FX,YNEW,HD2,IJOB,.TRUE.)
 358       DO I=1,N
 359          YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)
 360       END DO
 361       CALL FCN(N,X+C3*H,YNEW,DY)
 362       DO I=1,N
 363          YNEW(I)=HC31*AK1(I)+HC32*AK2(I)
 364       END DO
 365       CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
 366      &    M1,M2,NM1,FAC,E,LDE,IP,DY,AK3,FX,YNEW,HD3,IJOB,.TRUE.)
 367       DO I=1,N
 368          YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)
 369       END DO
 370       CALL FCN(N,X+C4*H,YNEW,DY)
 371       DO I=1,N
 372          YNEW(I)=HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I)
 373       END DO
 374       CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
 375      &    M1,M2,NM1,FAC,E,LDE,IP,DY,AK4,FX,YNEW,HD4,IJOB,.TRUE.)
 376       DO I=1,N
 377          YNEW(I)=Y(I)+A51*AK1(I)+A52*AK2(I)+A53*AK3(I)+A54*AK4(I)
 378       END DO
 379       CALL FCN(N,X+H,YNEW,DY)
 380       DO I=1,N
 381          AK6(I)=HC52*AK2(I)+HC54*AK4(I)+HC51*AK1(I)+HC53*AK3(I)
 382       END DO
 383       CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
 384      &    M1,M2,NM1,FAC,E,LDE,IP,DY,AK5,FX,AK6,0.D0,IJOB,.TRUE.)
 385 C ------------ EMBEDDED SOLUTION ---------------
 386       DO I=1,N
 387          YNEW(I)=YNEW(I)+AK5(I)
 388       END DO
 389       CALL FCN(N,X+H,YNEW,DY)
 390       DO I=1,N
 391          AK5(I)=HC61*AK1(I)+HC62*AK2(I)+HC65*AK5(I)
 392      &              +HC64*AK4(I)+HC63*AK3(I)
 393       END DO
 394       CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
 395      &    M1,M2,NM1,FAC,E,LDE,IP,DY,AK6,FX,AK5,0.D0,IJOB,.TRUE.)
 396 C ------------ NEW SOLUTION ---------------
 397       DO  I=1,N
 398          YNEW(I)=YNEW(I)+AK6(I)
 399       END DO
 400       NSOL=NSOL+6
 401       NFCN=NFCN+5
 402
 403 C *** *** *** *** *** *** ***
 404 C  ERROR ESTIMATION
 405 C *** *** *** *** *** *** ***
 406       NSTEP=NSTEP+1
 407 C ------------ COMPUTE ERROR ESTIMATION ----------------
 408       ERR=0.0D0
 409       DO I=1,N
 410          IF (ITOL.EQ.0) THEN
 411             SK=AbsTol(1)+RelTol(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
 412          ELSE
 413             SK=AbsTol(I)+RelTol(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
 414          END IF
 415          ERR=ERR+(AK6(I)/SK)**2
 416 c2           ERR = DMAX1(ERR, AK6(I)/SK)
 417       END DO
 418       ERR=DSQRT(ERR/N)
 419
 420 C --- COMPUTATION OF HNEW
 421 C --- WE REQUIRE .2<=HNEW/H<=6.
 422       FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**0.25D0/SAFE))
 423       HNEW=DMAX1(H/FAC, STEPMIN)
 424
 425 C *** *** *** *** *** *** ***
 426 C  IS THE ERROR SMALL ENOUGH ?
 427 C *** *** *** *** *** *** ***
 428
 429       IF ( (ERR.LE.1.D0).or.(H.LE.STEPMIN) ) THEN
 430 C --- STEP IS ACCEPTED
 431          NACCPT=NACCPT+1
 432          IF (PRED) THEN
 433 C       --- PREDICTIVE CONTROLLER OF GUSTAFSSON
 434             IF (NACCPT.GT.1) THEN
 435                FACGUS=(HACC/H)*(ERR**2/ERRACC)**0.25D0/SAFE
 436                FACGUS=DMAX1(FAC2,DMIN1(FAC1,FACGUS))
 437                FAC=DMAX1(FAC,FACGUS)
 438                HNEW=DMAX1(H/FAC, STEPMIN)
 439             END IF
 440             HACC=H
 441             ERRACC=DMAX1(1.0D-2,ERR)
 442          END IF
 443          DO I=1,N
 444             Y(I)=YNEW(I)
 445          END DO
 446          XOLD=X
 447          X=X+H
 448          IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN
 449          IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H))
 450          REJECT=.FALSE.
 451          H=HNEW
 452       ELSE
 453 C --- STEP IS REJECTED
 454          REJECT=.TRUE.
 455          LAST=.FALSE.
 456          H=HNEW
 457          IF (NACCPT.GE.1) NREJCT=NREJCT+1
 458       END IF
 459       END DO
 460       RETURN
 461       END
 462 C
 463       SUBROUTINE FAIL_EXIT(NERR,X,IDID,H,NMAX)
 464       INTEGER NERR, NMAX
 465       KPP_REAL X, H
 466       GO TO (1,2,3,4) NERR
 467  1    CONTINUE
 468       WRITE(6,979)X
 469       WRITE(6,*) ' MATRIX IS REPEATEDLY SINGULAR, IER=',IER
 470       IDID=-4
 471       STOP
 472  2    CONTINUE
 473       WRITE(6,979)X
 474       WRITE(6,*) ' STEP SIZE TOO SMALL, H=',H
 475       IDID=-3
 476       STOP
 477  3    CONTINUE
 478       WRITE(6,979)X
 479       WRITE(6,*) ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED'
 480       IDID=-2
 481       STOP
 482 C --- EXIT CAUSED BY solout
 483  4    CONTINUE
 484       WRITE(6,979)X
 485  979  FORMAT(' EXIT OF RODAS AT X=',E18.4)
 486       IDID=2
 487       RETURN
 488       END
 489
 490       SUBROUTINE ROCOE(METH,A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,
 491      &  C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,
 492      &  C62,C63,C64,C65,GAMMA,C2,C3,C4,D1,D2,D3,D4,
 493      &  D21,D22,D23,D24,D25,D31,D32,D33,D34,D35)
 494       IMPLICIT KPP_REAL (A-H,O-Z)
 495
 496       if (METH.ne.1) print *, 'WRONG CHOICE OF METHOD'
 497         C2=0.386D0
 498         C3=0.21D0
 499         C4=0.63D0
 500         BET2P=0.0317D0
 501         BET3P=0.0635D0
 502         BET4P=0.3438D0
 503        D1= 0.2500000000000000D+00
 504        D2=-0.1043000000000000D+00
 505        D3= 0.1035000000000000D+00
 506        D4=-0.3620000000000023D-01
 507        A21= 0.1544000000000000D+01
 508        A31= 0.9466785280815826D+00
 509        A32= 0.2557011698983284D+00
 510        A41= 0.3314825187068521D+01
 511        A42= 0.2896124015972201D+01
 512        A43= 0.9986419139977817D+00
 513        A51= 0.1221224509226641D+01
 514        A52= 0.6019134481288629D+01
 515        A53= 0.1253708332932087D+02
 516        A54=-0.6878860361058950D+00
 517        C21=-0.5668800000000000D+01
 518        C31=-0.2430093356833875D+01
 519        C32=-0.2063599157091915D+00
 520        C41=-0.1073529058151375D+00
 521        C42=-0.9594562251023355D+01
 522        C43=-0.2047028614809616D+02
 523        C51= 0.7496443313967647D+01
 524        C52=-0.1024680431464352D+02
 525        C53=-0.3399990352819905D+02
 526        C54= 0.1170890893206160D+02
 527        C61= 0.8083246795921522D+01
 528        C62=-0.7981132988064893D+01
 529        C63=-0.3152159432874371D+02
 530        C64= 0.1631930543123136D+02
 531        C65=-0.6058818238834054D+01
 532        GAMMA= 0.2500000000000000D+00
 533
 534        D21= 0.1012623508344586D+02
 535        D22=-0.7487995877610167D+01
 536        D23=-0.3480091861555747D+02
 537        D24=-0.7992771707568823D+01
 538        D25= 0.1025137723295662D+01
 539        D31=-0.6762803392801253D+00
 540        D32= 0.6087714651680015D+01
 541        D33= 0.1643084320892478D+02
 542        D34= 0.2476722511418386D+02
 543        D35=-0.6594389125716872D+01
 544       RETURN
 545       END
 546 C
 547
 548 C ******************************************
 549 C     VERSION OF SEPTEMBER 18, 1995
 550 C ******************************************
 551 C
 552       SUBROUTINE DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
 553      &            M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES)
 554       IMPLICIT KPP_REAL (A-H,O-Z)
 555       INCLUDE 'KPP_ROOT_params.h'
 556       INCLUDE 'KPP_ROOT_sparse.h'
 557       DIMENSION FJAC(LU_NONZERO),FMAS(LDMAS,NM1),E1(LU_NONZERO),
 558      &          IP1(NM1),IPHES(N)
 559       LOGICAL CALHES
 560       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
 561 C
 562
 563
 564 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
 565       DO J=1,LU_NONZERO
 566          E1(J) = -FJAC(J)
 567       END DO
 568       DO J=1,N
 569          E1(LU_DIAG(J)) = E1(LU_DIAG(J)) + FAC1
 570       END DO
 571       CALL KppDecomp (E1,IER)
 572       RETURN
 573       END
 574 C
 575 C     END OF SUBROUTINE DECOMR
 576 C
 577 C ***********************************************************
 578 C
 579 C
 580 C
 581 C
 582       SUBROUTINE SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
 583      &          M1,M2,NM1,FAC1,E,LDE,IP,DY,AK,FX,YNEW,HD,IJOB,STAGE1)
 584       IMPLICIT KPP_REAL (A-H,O-Z)
 585       INCLUDE 'KPP_ROOT_params.h'
 586       INCLUDE 'KPP_ROOT_sparse.h'
 587       DIMENSION FJAC(LU_NONZERO),FMAS(LDMAS,NM1),E(LU_NONZERO),
 588      &          IP(NM1),DY(N),AK(N),FX(N),YNEW(N)
 589       LOGICAL STAGE1
 590       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
 591 C
 592       IF (HD.EQ.0.D0) THEN
 593          DO  I=1,N
 594            AK(I)=DY(I)
 595          END DO
 596       ELSE
 597          DO I=1,N
 598             AK(I)=DY(I)+HD*FX(I)
 599          END DO
 600       END IF
 601
 602 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
 603       IF (STAGE1) THEN
 604          DO I=1,N
 605             AK(I)=AK(I)+YNEW(I)
 606          END DO
 607       END IF
 608       CALL KppSolve (E,AK)
 609       RETURN
 610       END
 611 C
 612 C     END OF SUBROUTINE SLVROD
 613 C
 614 C
 615 C ***********************************************************
 616
 617
 618       SUBROUTINE FUNC_CHEM(N, T, Y, P)
 619       INCLUDE 'KPP_ROOT_params.h'
 620       INCLUDE 'KPP_ROOT_global.h'
 621       INTEGER N
 622       KPP_REAL   T, Told
 623       KPP_REAL   Y(NVAR), P(NVAR)
 624       Told = TIME
 625       TIME = T
 626       CALL Update_SUN()
 627       CALL Update_RCONST()
 628       CALL Fun( Y,  FIX, RCONST, P )
 629       TIME = Told
 630       RETURN
 631       END
 632
 633
 634       SUBROUTINE JAC_CHEM(N, T, Y, J)
 635       INCLUDE 'KPP_ROOT_params.h'
 636       INCLUDE 'KPP_ROOT_global.h'
 637       INTEGER N
 638       KPP_REAL   Told, T
 639       KPP_REAL   Y(NVAR), J(LU_NONZERO)
 640       Told = TIME
 641       TIME = T
 642       CALL Update_SUN()
 643       CALL Update_RCONST()
 644       CALL Jac_SP( Y,  FIX, RCONST, J )
 645       TIME = Told
 646       RETURN
 647       END