chem/KPP/kpp/kpp-2.1/int/kpp_sdirk.f

   1       SUBROUTINE INTEGRATE( TIN, TOUT )
   2
   3       IMPLICIT KPP_REAL (A-H,O-Z)
   4       INCLUDE 'KPP_ROOT_Parameters.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+12*NVAR+7,LIWORK=2*NVAR+7)
  14       PARAMETER (LRCONT=5*NVAR+2)
  15
  16       KPP_REAL WORK(LWORK)
  17       INTEGER IWORK(LIWORK)
  18       COMMON /CONT/ ICONT(4),RCONT(LRCONT)
  19       EXTERNAL FUNC_CHEM,JAC_CHEM
  20
  21       ITOL=1     ! --- VECTOR TOLERANCES
  22       IJAC=1     ! --- COMPUTE THE JACOBIAN ANALYTICALLY
  23       IMAS=0     ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
  24
  25       DO i=1,20
  26         IWORK(i) = 0
  27         WORK(i) = 0.D0
  28       ENDDO
  29
  30       IWORK(3) = 8
  31
  32       CALL ATMSDIRK(NVAR,FUNC_CHEM,TIN,VAR,TOUT,STEPMIN,
  33      &                  RTOL,ATOL,ITOL,
  34      &                  JAC_CHEM ,IJAC, FUNC_CHEM ,IMAS,
  35      &                  WORK,LWORK,IWORK,LIWORK,LRCONT,IDID)
  36
  37       IF (IDID.LT.0) THEN
  38         print *,'ATMSDIRK: Unsucessfull exit at T=',
  39      &          TIN,' (IDID=',IDID,')'
  40       ENDIF
  41
  42       RETURN
  43       END
  44
  45
  46       SUBROUTINE ATMSDIRK(N,FCN,X,Y,XEND,H,
  47      &                  RelTol,AbsTol,ITOL,
  48      &                  JAC ,IJAC, MAS ,IMAS,
  49      &                  WORK,LWORK,IWORK,LIWORK,LRCONT,IDID)
  50 C ----------------------------------------------------------
  51 C *** *** *** *** *** *** *** *** *** *** *** *** ***
  52 C          DECLARATIONS
  53 C *** *** *** *** *** *** *** *** *** *** *** *** ***
  54       IMPLICIT KPP_REAL (A-H,O-Z)
  55       DIMENSION Y(N),AbsTol(1),RelTol(1),WORK(LWORK),IWORK(LIWORK)
  56       LOGICAL IMPLCT,JBAND,ARRET
  57       EXTERNAL FCN,JAC,MAS
  58       COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
  59
  60 C *** *** *** *** *** *** ***
  61 C        SETTING THE PARAMETERS
  62 C *** *** *** *** *** *** ***
  63        NFCN=0
  64        NJAC=0
  65        NSTEP=0
  66        NACCPT=0
  67        NREJCT=0
  68        NDEC=0
  69        NSOL=0
  70        ARRET=.FALSE.
  71 C -------- SWITCH FOR TRANSFORMATION OF JACOBIAN TO HESS_CHEM FORM ---
  72       NHESS1 = 0         ! ADRIAN
  73 C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
  74       IF(IWORK(2).EQ.0)THEN
  75          NMAX=100000
  76       ELSE
  77          NMAX=IWORK(2)
  78          IF(NMAX.LE.0)THEN
  79             WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2)
  80             ARRET=.TRUE.
  81          END IF
  82       END IF
  83 C -------- NIT    MAXIMAL NUMBER OF NEWTON ITERATIONS
  84       IF(IWORK(3).EQ.0)THEN
  85          NIT=8
  86       ELSE
  87          NIT=IWORK(3)
  88          IF(NIT.LE.0)THEN
  89             WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
  90             ARRET=.TRUE.
  91          END IF
  92       END IF
  93 C -------- METH    SWITCH FOR THE COEFFICIENTS OF THE METHOD
  94       METH = 2
  95 C -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0
  96       IF(WORK(1).EQ.0.D0)THEN
  97          UROUND=1.D-16
  98       ELSE
  99          UROUND=WORK(1)
 100          IF(UROUND.LE.1.D-19.OR.UROUND.GE.1.D0)THEN
 101             WRITE(6,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1)
 102             ARRET=.TRUE.
 103          END IF
 104       END IF
 105 C --------- SAFE     SAFETY FACTOR IN STEP SIZE PREDICTION
 106       IF(WORK(2).EQ.0.D0)THEN
 107          SAFE=0.9D0
 108       ELSE
 109          SAFE=WORK(2)
 110          IF(SAFE.LE..001D0.OR.SAFE.GE.1.D0)THEN
 111             WRITE(6,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2)
 112             ARRET=.TRUE.
 113          END IF
 114       END IF
 115 C ------ THET     DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
 116       IF(WORK(3).EQ.0.D0)THEN
 117          THET=0.001D0
 118       ELSE
 119          THET=WORK(3)
 120       END IF
 121 C --- FNEWT   STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1.
 122       IF(WORK(4).EQ.0.D0)THEN
 123          FNEWT=0.03D0
 124       ELSE
 125          FNEWT=WORK(4)
 126       END IF
 127 C --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST.
 128       IF(WORK(5).EQ.0.D0)THEN
 129          QUOT1=1.D0
 130       ELSE
 131          QUOT1=WORK(5)
 132       END IF
 133       IF(WORK(6).EQ.0.D0)THEN
 134          QUOT2=1.2D0
 135       ELSE
 136          QUOT2=WORK(6)
 137       END IF
 138 C -------- MAXIMAL STEP SIZE
 139       IF(WORK(7).EQ.0.D0)THEN
 140          HMAX=XEND-X
 141       ELSE
 142          HMAX=WORK(7)
 143       END IF
 144 C --------- CHECK IF TOLERANCES ARE O.K.
 145       IF (ITOL.EQ.0) THEN
 146           IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN
 147               WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
 148               ARRET=.TRUE.
 149           END IF
 150       ELSE
 151           DO 15 I=1,N
 152           IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN
 153               WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
 154               ARRET=.TRUE.
 155           END IF
 156   15      CONTINUE
 157       END IF
 158
 159 C *** *** *** *** *** *** *** *** *** *** *** *** ***
 160 C         COMPUTATION OF ARRAY ENTRIES
 161 C *** *** *** *** *** *** *** *** *** *** *** *** ***
 162 C ---- IMPLICIT, BANDED OR NOT ?
 163       IMPLCT=IMAS.NE.0
 164       ARRET=.FALSE.
 165 C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
 166 C -- JACOBIAN
 167          LDJAC=N
 168 C -- MATRIX E FOR LINEAR ALGEBRA
 169          LDE=N
 170 C -- MASS MATRIX
 171       IF (IMPLCT) THEN
 172           print *,'IMPLCT 1'
 173       ELSE
 174           LDMAS=0
 175       END IF
 176       LDMAS2=MAX(1,LDMAS)
 177
 178 C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
 179       IEYHAT=8
 180       IEZ=IEYHAT+N
 181       IEY0=IEZ+N
 182       IEZ1=IEY0+N
 183       IEZ2=IEZ1+N
 184       IEZ3=IEZ2+N
 185       IEZ4=IEZ3+N
 186       IEZ5=IEZ4+N
 187       IESCAL=IEZ5+N
 188       IEF1=IESCAL+N
 189       IEG1=IEF1+N
 190       IEH1=IEG1+N
 191       IEJAC=IEH1+N
 192       IEMAS=IEJAC+N*LDJAC
 193       IEE=IEMAS+N*LDMAS
 194
 195 C ------ TOTAL STORAGE REQUIREMENT -----------
 196       ISTORE=IEE+N*LDE-1
 197       IF(ISTORE.GT.LWORK)THEN
 198          WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
 199          ARRET=.TRUE.
 200       END IF
 201 C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
 202       IEIP=5
 203       IEHES=IEIP+N
 204 C --------- TOTAL REQUIREMENT ---------------
 205       ISTORE=IEHES+N-1
 206       IF(ISTORE.GT.LIWORK)THEN
 207          WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
 208          ARRET=.TRUE.
 209       END IF
 210 C --------- CONTROL OF LENGTH OF COMMON BLOCK "CONT" -------
 211       IF(LRCONT.LT.(5*N+2))THEN
 212          WRITE(6,*)' INSUFF. STORAGE FOR RCONT, MIN. LRCONT=',5*N+2
 213          ARRET=.TRUE.
 214       END IF
 215 C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
 216       IF (ARRET) THEN
 217          IDID=-1
 218          RETURN
 219       END IF
 220 C -------- CALL TO CORE INTEGRATOR ------------
 221       CALL SDICOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,
 222      &   JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,IOUT,IDID,
 223      &   NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,METH,NHESS1,
 224      &   IMPLCT,JBAND,LDJAC,LDE,LDMAS2,
 225      &   WORK(IEYHAT),WORK(IEZ),WORK(IEY0),WORK(IEZ1),WORK(IEZ2),
 226      &   WORK(IEZ3),WORK(IEZ4),WORK(IEZ5),WORK(IESCAL),WORK(IEF1),
 227      &   WORK(IEG1),WORK(IEH1),WORK(IEJAC),WORK(IEE),
 228      &   WORK(IEMAS),IWORK(IEIP),IWORK(IEHES))
 229 C ----------- RETURN -----------
 230       RETURN
 231       END
 232 C
 233 C
 234 C
 235 C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
 236 C
 237       SUBROUTINE SDICOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,
 238      &   JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,IOUT,IDID,
 239      &   NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,METH,NHESS1,
 240      &   IMPLCT,BANDED,LDJAC,LE,LDMAS,
 241      &   YHAT,Z,Y0,Z1,Z2,Z3,Z4,Z5,SCAL,F1,G1,H1,FJAC,E,FMAS,IP,IPHES)
 242 C ----------------------------------------------------------
 243 C     CORE INTEGRATOR FOR SDIRK4
 244 C     PARAMETERS SAME AS IN SDIRK4 WITH WORKSPACE ADDED
 245 C ----------------------------------------------------------
 246 C         DECLARATIONS
 247 C ----------------------------------------------------------
 248       IMPLICIT KPP_REAL (A-H,O-Z)
 249       INCLUDE 'KPP_ROOT_Parameters.h'
 250       INCLUDE 'KPP_ROOT_Sparse.h'
 251       KPP_REAL Y(N),YHAT(N),Z(N),Y0(N),Z1(N),Z2(N),Z3(N),Z4(N),Z5(N)
 252       KPP_REAL SCAL(N),F1(N),G1(N),H1(N)
 253       KPP_REAL FJAC(LU_NONZERO),E(LU_NONZERO),FMAS(LDMAS,N)
 254       KPP_REAL AbsTol(1),RelTol(1)
 255       INTEGER IP(N),IPHES(N)
 256       LOGICAL REJECT,FIRST,IMPLCT,BANDED,CALJAC,NEWTRE
 257       COMMON /CONT/NN,NN2,NN3,NN4,XOLD,HSOL,CONT(5*NVAR)
 258       COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
 259       EXTERNAL MAS, FCN, JAC
 260
 261 C *** *** *** *** *** *** ***
 262 C  INITIALISATIONS
 263 C *** *** *** *** *** *** ***
 264
 265 C --------- DUPLIFY N FOR COMMON BLOCK CONT -----
 266       NN=N
 267       NN2=2*N
 268       NN3=3*N
 269       NN4=4*N
 270
 271 C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
 272       IF(IMPLCT) CALL MAS(N,FMAS,LDMAS)
 273
 274 C ---------- CONSTANTS ---------
 275       MBDIAG=MUMAS+1
 276       IF (METH.EQ.2) THEN
 277 C ---------- METHOD WITH GAMMA = 4/15 ---------------
 278           GAMMA=4.0D0/15.0D0
 279           C2=23.0D0/30.0D0
 280           C3=17.0D0/30.0D0
 281           C4=2881.0D0/28965.0D0+GAMMA
 282           ALPH21=15.0D0/8.0D0
 283           ALPH31=1577061.0D0/922880.0D0
 284           ALPH32=-23427.0D0/115360.0D0
 285           ALPH41=647163682356923881.0D0/2414496535205978880.0D0
 286           ALPH42=-593512117011179.0D0/3245291041943520.0D0
 287           ALPH43=559907973726451.0D0/1886325418129671.0D0
 288           ALPH51=724545451.0D0/796538880.0D0
 289           ALPH52=-830832077.0D0/267298560.0D0
 290           ALPH53=30957577.0D0/2509272.0D0
 291           ALPH54=-69863904375173.0D0/6212571137048.0D0
 292           E1=7752107607.0D0/11393456128.0D0
 293           E2=-17881415427.0D0/11470078208.0D0
 294           E3=2433277665.0D0/179459416.0D0
 295           E4=-96203066666797.0D0/6212571137048.0D0
 296           D11= 24.74416644927758D0
 297           D12= -4.325375951824688D0
 298           D13= 41.39683763286316D0
 299           D14= -61.04144619901784D0
 300           D15= -3.391332232917013D0
 301           D21= -51.98245719616925D0
 302           D22= 10.52501981094525D0
 303           D23= -154.2067922191855D0
 304           D24= 214.3082125319825D0
 305           D25= 14.71166018088679D0
 306           D31= 33.14347947522142D0
 307           D32= -19.72986789558523D0
 308           D33= 230.4878502285804D0
 309           D34= -287.6629744338197D0
 310           D35= -18.99932366302254D0
 311           D41= -5.905188728329743D0
 312           D42= 13.53022403646467D0
 313           D43= -117.6778956422581D0
 314           D44= 134.3962081008550D0
 315           D45= 8.678995715052762D0
 316          ETA1=23.D0/8.D0
 317          ANU1= 0.9838473040915402D0
 318          ANU2= 0.3969226768377252D0
 319          AMU1= 0.6563374010466914D0
 320          AMU3= 0.3372498196189311D0
 321       ELSE
 322          PRINT *, 'WRONG  CHOICE OF <METH>'
 323       END IF
 324       POSNEG=SIGN(1.D0,XEND-X)
 325       HMAX1=MIN(ABS(HMAX),ABS(XEND-X))
 326       IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
 327       H=MIN(ABS(H),HMAX1)
 328       H=SIGN(H,POSNEG)
 329       HOLD=H
 330       CFAC=SAFE*(1+2*NIT)
 331       NEWTRE=.FALSE.
 332       REJECT=.FALSE.
 333       FIRST=.TRUE.
 334       FACCO1=1.D0
 335       FACCO2=1.D0
 336       FACCO3=1.D0
 337       FACCO4=1.D0
 338       FACCO5=1.D0
 339       NSING=0
 340       XOLD=X
 341       IF (ITOL.EQ.0) THEN
 342           DO 8 I=1,N
 343    8      SCAL(I)=1.D0 / ( AbsTol(1)+RelTol(1)*DABS(Y(I)) )
 344       ELSE
 345           DO 9 I=1,N
 346    9      SCAL(I)=1.D0 / ( AbsTol(I)+RelTol(I)*DABS(Y(I)) )
 347       END IF
 348
 349 C --- BASIC INTEGRATION STEP
 350   10  CONTINUE
 351
 352 C *** *** *** *** *** *** ***
 353 C  COMPUTATION OF THE JACOBIAN
 354 C *** *** *** *** *** *** ***
 355       NJAC=NJAC+1
 356       CALL JAC(N,X,Y,FJAC)
 357       CALJAC=.TRUE.
 358   20  CONTINUE
 359
 360 C *** *** *** *** *** *** ***
 361 C  COMPUTE THE MATRIX E AND ITS DECOMPOSITION
 362 C *** *** *** *** *** *** ***
 363       FAC1=1.D0/(H*GAMMA)
 364       IF (IMPLCT) THEN
 365          print *, 'IMPLCT 4'
 366       ELSE  ! EXPLICIT SYSTEM
 367 C --- THE MATRIX E (MAS=IDENTITY, JACOBIAN A FULL MATRIX)
 368 c         DO 526 J=1,N
 369 c           DO 525 I=1,N
 370 c 525         E(I,J)=-FJAC(I,J)
 371 c 526       E(J,J)=E(J,J)+FAC1
 372 c         CALL DEC(N,LE,E,IP,IER)
 373           DO K=1,LU_NONZERO
 374              E(K) = -FJAC(K)
 375           END DO
 376           DO I=1,N
 377              IDG = LU_DIAG(I)
 378              E(IDG) = E(IDG) + FAC1
 379           END DO
 380           CALL KppDecomp ( E, IER)
 381
 382          IF (IER.NE.0) GOTO 79
 383       END IF
 384       NDEC=NDEC+1
 385   30  CONTINUE
 386
 387       IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79
 388       XPH=X+H
 389 C --- LOOP FOR THE 5 STAGES
 390       FACCO1=DMAX1(FACCO1,UROUND)**0.8D0
 391       FACCO2=DMAX1(FACCO2,UROUND)**0.8D0
 392       FACCO3=DMAX1(FACCO3,UROUND)**0.8D0
 393       FACCO4=DMAX1(FACCO4,UROUND)**0.8D0
 394       FACCO5=DMAX1(FACCO5,UROUND)**0.8D0
 395
 396 C *** *** *** *** *** *** ***
 397 C  STARTING VALUES FOR NEWTON ITERATION
 398 C *** *** *** *** *** *** ***
 399       DO 59 ISTAGE=1,5
 400       IF (ISTAGE.EQ.1) THEN
 401           XCH=X+GAMMA*H
 402           IF (FIRST.OR.NEWTRE) THEN
 403               DO 132 I=1,N
 404  132          Z(I)=0.D0
 405           ELSE
 406               S=1.D0+GAMMA*H/HOLD
 407               DO 232 I=1,N
 408 c  232          Z(I) = 0.D0
 409  232          Z(I)=S*(CONT(I+NN)+S*(CONT(I+NN2)+S*(CONT(I+NN3)
 410      &             +S*CONT(I+NN4))))-YHAT(I)
 411
 412           END IF
 413           DO 31 I=1,N
 414   31      G1(I)=0.D0
 415           FACCON=FACCO1
 416       END IF
 417       IF (ISTAGE.EQ.2) THEN
 418           XCH=X+C2*H
 419           DO 131 I=1,N
 420           Z1I=Z1(I)
 421           Z(I)=ETA1*Z1I
 422  131      G1(I)=ALPH21*Z1I
 423           FACCON=FACCO2
 424       END IF
 425       IF (ISTAGE.EQ.3) THEN
 426           XCH=X+C3*H
 427           DO 231 I=1,N
 428           Z1I=Z1(I)
 429           Z2I=Z2(I)
 430           Z(I)=ANU1*Z1I+ANU2*Z2I
 431  231      G1(I)=ALPH31*Z1I+ALPH32*Z2I
 432           FACCON=FACCO3
 433       END IF
 434       IF (ISTAGE.EQ.4) THEN
 435           XCH=X+C4*H
 436           DO 331 I=1,N
 437           Z1I=Z1(I)
 438           Z3I=Z3(I)
 439           Z(I)=AMU1*Z1I+AMU3*Z3I
 440  331      G1(I)=ALPH41*Z1I+ALPH42*Z2(I)+ALPH43*Z3I
 441           FACCON=FACCO4
 442       END IF
 443       IF (ISTAGE.EQ.5) THEN
 444           XCH=XPH
 445           DO 431 I=1,N
 446           Z1I=Z1(I)
 447           Z2I=Z2(I)
 448           Z3I=Z3(I)
 449           Z4I=Z4(I)
 450           Z(I)=E1*Z1I+E2*Z2I+E3*Z3I+E4*Z4I
 451           YHAT(I)=Z(I)
 452  431      G1(I)=ALPH51*Z1I+ALPH52*Z2I+ALPH53*Z3I+ALPH54*Z4I
 453           FACCON=FACCO5
 454       END IF
 455
 456
 457
 458 C *** *** *** *** *** *** *** *** *** *** ***
 459 C  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
 460 C *** *** *** *** *** *** *** *** *** *** ***
 461             NEWT=0
 462             THETA=ABS(THET)
 463             IF (REJECT) THETA=2*ABS(THET)
 464   40        CONTINUE
 465             IF (NEWT.GE.NIT) THEN
 466                 H=H/2.D0
 467                 REJECT=.TRUE.
 468                 NEWTRE=.TRUE.
 469                 IF (CALJAC) GOTO 20
 470                 GOTO 10
 471             END IF
 472
 473 C ---     COMPUTE THE RIGHT-HAND SIDE
 474             DO 41 I=1,N
 475             H1(I)=G1(I)-Z(I)
 476   41        CONT(I)=Y(I)+Z(I)
 477             CALL FCN(N,XCH,CONT,F1)
 478             NFCN=NFCN+1
 479
 480 C ---     KppSolve THE LINEAR SYSTEMS
 481             IF (IMPLCT) THEN
 482                 print *, 'IMPLCT 2'
 483             ELSE
 484                 DO 345 I=1,N
 485  345            F1(I)=H1(I)*FAC1+F1(I)
 486 C                CALL SOL(N,LE,E,F1,IP)
 487                 CALL KppSolve(E, F1)
 488             END IF
 489             NEWT=NEWT+1
 490             DYNO=0.D0
 491 C --- NORM 2 ---
 492             DO 57 I=1,N
 493   57        DYNO=DYNO+(F1(I)*SCAL(I))**2
 494             DYNO=DSQRT(DYNO/N)
 495 C --- NORM INF ---
 496 C            DO 57 I=1,N
 497 C  57        DYNO=DMAX1( DYNO, DABS(F1(I)*SCAL(I)) )
 498
 499
 500 C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE
 501             IF (NEWT.GE.2.AND.NEWT.LT.NIT) THEN
 502                 THETA=DYNO/DYNOLD
 503                 IF (THETA.LT.0.99D0) THEN
 504                     FACCON=THETA/(1.0D0-THETA)
 505                     DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)
 506                     QNEWT=DMAX1(1.0D-4,DMIN1(16.0D0,DYTH/FNEWT))
 507                     IF (QNEWT.GE.1.0D0) THEN
 508                          H=.8D0*H*QNEWT**(-1.0D0/(NIT-NEWT))
 509                          REJECT=.TRUE.
 510                          NEWTRE=.TRUE.
 511                          IF (CALJAC) GOTO 20
 512                          GOTO 10
 513                     END IF
 514                 ELSE
 515                     NEWTRE=.TRUE.
 516                     GOTO 78
 517                 END IF
 518             END IF
 519             DYNOLD=DYNO
 520             DO 58 I=1,N
 521   58        Z(I)=Z(I)+F1(I)
 522             NSOL=NSOL+1
 523             IF (FACCON*DYNO.GT.FNEWT) GOTO 40
 524
 525 C --- END OF SIMPILFIED NEWTON
 526       IF (ISTAGE.EQ.1) THEN
 527           DO I=1,N
 528             Z1(I) = Z(I)
 529           END DO
 530           FACCO1=FACCON
 531       END IF
 532       IF (ISTAGE.EQ.2) THEN
 533           DO I=1,N
 534             Z2(I) = Z(I)
 535           END DO
 536           FACCO2=FACCON
 537       END IF
 538       IF (ISTAGE.EQ.3) THEN
 539           DO I=1,N
 540             Z3(I) = Z(I)
 541           END DO
 542           FACCO3=FACCON
 543       END IF
 544       IF (ISTAGE.EQ.4) THEN
 545           DO I=1,N
 546             Z4(I) = Z(I)
 547           END DO
 548           FACCO4=FACCON
 549       END IF
 550       IF (ISTAGE.EQ.5) THEN
 551           DO I=1,N
 552             Z5(I) = Z(I)
 553           END DO
 554           FACCO5=FACCON
 555       END IF
 556   59  CONTINUE
 557
 558
 559 C *** *** *** *** *** *** ***
 560 C  ERROR ESTIMATION
 561 C *** *** *** *** *** *** ***
 562       NSTEP=NSTEP+1
 563       IF (IMPLCT) THEN
 564           print *,'IMPLCT 3'
 565       ELSE
 566           DO 461 I=1,N
 567  461      CONT(I)=FAC1*(Z5(I)-YHAT(I))
 568       END IF
 569
 570       CALL KppSolve(E, CONT)
 571
 572       ERR=0.D0
 573 C ---- NORM 2 ---
 574       DO 64 I=1,N
 575   64  ERR=ERR+(CONT(I)*SCAL(I))**2
 576       ERR=DMAX1(DSQRT(ERR/N),1.D-10)
 577
 578 C ---- NORM INF ---
 579 C      DO 64 I=1,N
 580 c  64  ERR=DMAX1( ERR, DABS( CONT(I)*SCAL(I) ) )
 581
 582 C --- COMPUTATION OF HNEW
 583 C --- WE REQUIRE .25<=HNEW/H<=10.
 584       FAC=DMIN1(SAFE,CFAC/(NEWT+2*NIT))
 585       QUOT=DMAX1(.25D0,DMIN1(10.D0,(ERR)**.25D0/FAC))
 586       HNEW= H/QUOT
 587
 588 C *** *** *** *** *** *** ***
 589 C  IS THE ERROR SMALL ENOUGH ?
 590 C *** *** *** *** *** *** ***
 591       IF (ERR.LT.1.D0) THEN
 592 C --- STEP IS ACCEPTED
 593          FIRST=.FALSE.
 594          NACCPT=NACCPT+1
 595          HOLD=H
 596          XOLD=X
 597 C --- COEFFICIENTS FOR CONTINUOUS SOLUTION
 598          DO 74 I=1,N
 599           Z1I=Z1(I)
 600           Z2I=Z2(I)
 601           Z3I=Z3(I)
 602           Z4I=Z4(I)
 603           Z5I=Z5(I)
 604          CONT(I)=Y(I)
 605          Y(I)=Y(I)+Z5I
 606          CONT(I+NN) =D11*Z1I+D12*Z2I+D13*Z3I+D14*Z4I+D15*Z5I
 607          CONT(I+NN2)=D21*Z1I+D22*Z2I+D23*Z3I+D24*Z4I+D25*Z5I
 608          CONT(I+NN3)=D31*Z1I+D32*Z2I+D33*Z3I+D34*Z4I+D35*Z5I
 609          CONT(I+NN4)=D41*Z1I+D42*Z2I+D43*Z3I+D44*Z4I+D45*Z5I
 610          YHAT(I)=Z5I
 611          IF (ITOL.EQ.0) THEN
 612            SCAL(I)=1.D0/( AbsTol(1)+RelTol(1)*DABS(Y(I)) )
 613          ELSE
 614            SCAL(I)=1.D0/( AbsTol(I)+RelTol(I)*DABS(Y(I)) )
 615          END IF
 616   74     CONTINUE
 617          X=XPH
 618          CALJAC=.FALSE.
 619          IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN
 620             H=HOPT
 621             IDID=1
 622             RETURN
 623          END IF
 624          IF (IJAC.EQ.0) CALL FCN(N,X,Y,Y0)
 625          NFCN=NFCN+1
 626          HNEW=POSNEG*DMIN1(DABS(HNEW),HMAX1)
 627          HOPT=HNEW
 628          IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H))
 629          REJECT=.FALSE.
 630          NEWTRE=.FALSE.
 631          IF ((X+HNEW/QUOT1-XEND)*POSNEG.GT.0.D0) THEN
 632             H=XEND-X
 633          ELSE
 634             QT=HNEW/H
 635             IF (THETA.LE.THET.AND.QT.GE.QUOT1.AND.QT.LE.QUOT2) GOTO 30
 636             H = HNEW
 637          END IF
 638          IF (THETA.LE.THET) GOTO 20
 639          GOTO 10
 640
 641       ELSE
 642 C --- STEP IS REJECTED
 643          REJECT=.TRUE.
 644          IF (FIRST) THEN
 645              H=H/10.D0
 646          ELSE
 647              H=HNEW
 648          END IF
 649          IF (NACCPT.GE.1) NREJCT=NREJCT+1
 650          IF (CALJAC) GOTO 20
 651          GOTO 10
 652       END IF
 653
 654 C --- UNEXPECTED STEP-REJECTION
 655   78  CONTINUE
 656       IF (IER.NE.0) THEN
 657           WRITE (6,*) ' MATRIX IS SINGULAR, IER=',IER,' X=',X,' H=',H
 658           NSING=NSING+1
 659           IF (NSING.GE.6) GOTO 79
 660       END IF
 661       H=H*0.5D0
 662       REJECT=.TRUE.
 663       IF (CALJAC) GOTO 20
 664       GOTO 10
 665
 666 C --- FAIL EXIT
 667   79  WRITE (6,979) X,H,IER
 668  979  FORMAT(' EXIT OF SDIRK4 AT X=',D14.7,'   H=',D14.7,'   IER=',I4)
 669       IDID=-1
 670       RETURN
 671       END
 672 C
 673
 674
 675       SUBROUTINE FUNC_CHEM(N, T, Y, P)
 676       INCLUDE 'KPP_ROOT_Parameters.h'
 677       INCLUDE 'KPP_ROOT_Global.h'
 678       INTEGER N
 679       KPP_REAL   T, Told
 680       KPP_REAL   Y(NVAR), P(NVAR)
 681       Told = TIME
 682       TIME = T
 683       CALL Update_SUN()
 684       CALL Update_RCONST()
 685       CALL Fun( Y,  FIX, RCONST, P )
 686       TIME = Told
 687       RETURN
 688       END
 689
 690
 691       SUBROUTINE JAC_CHEM(N, T, Y, J)
 692       INCLUDE 'KPP_ROOT_Parameters.h'
 693       INCLUDE 'KPP_ROOT_Global.h'
 694       INTEGER N
 695       KPP_REAL   Told, T
 696       KPP_REAL   Y(NVAR), J(LU_NONZERO)
 697       Told = TIME
 698       TIME = T
 699       CALL Update_SUN()
 700       CALL Update_RCONST()
 701       CALL Jac_SP( Y,  FIX, RCONST, J )
 702       TIME = Told
 703       RETURN
 704       END
 705