chem/KPP/kpp/kpp-2.1/int/kpp_radau5.f90

   1 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
   2 !  RADAU5 - Runge-Kutta method based on Radau-2A quadrature               !
   3 !                       (2 stages, order 5)                               !
   4 !  By default the code employs the KPP sparse linear algebra routines     !
   5 !  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
   6 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
   7
   8 MODULE KPP_ROOT_Integrator
   9
  10   USE KPP_ROOT_Precision, ONLY: dp
  11   USE KPP_ROOT_Jacobian, ONLY: NVAR, LU_NONZERO, LU_DIAG
  12   USE KPP_ROOT_LinearAlgebra
  13
  14   IMPLICIT NONE
  15   PUBLIC
  16   SAVE
  17
  18   ! Statistics
  19   INTEGER :: Nfun, Njac, Nstp, Nacc, Nrej, Ndec, Nsol, Nsng
  20
  21   ! Method parameters
  22   KPP_REAL :: Transf(3,3), TransfInv(3,3),      &
  23                    rkA(3,3), rkB(3), rkC(3), rkE(3), &
  24                    rkGamma, rkAlpha, rkBeta
  25
  26   ! description of the error numbers IERR
  27   CHARACTER(LEN=50), PARAMETER, DIMENSION(-11:1) :: IERR_NAMES = (/ &
  28     'Matrix is repeatedly singular                     ', & ! -11
  29     'Step size too small: T + 10*H = T or H < Roundoff ', & ! -10
  30     'No of steps exceeds maximum bound                 ', & ! -9
  31     'Tolerances are too small                          ', & ! -8
  32     'Improper values for Qmin, Qmax                    ', & ! -7
  33     'Newton stopping tolerance too small               ', & ! -6
  34     'Improper value for ThetaMin                       ', & ! -5
  35     'Improper values for FacMin/FacMax/FacSafe/FacRej  ', & ! -4
  36     'Hmin/Hmax/Hstart must be positive                 ', & ! -3
  37     'Improper value for maximal no of Newton iterations', & ! -2
  38     'Improper value for maximal no of steps            ', & ! -1
  39     '                                                  ', & !  0 (not used)
  40     'Success                                           ' /) !  1
  41
  42 CONTAINS
  43
  44   ! **************************************************************************
  45
  46   SUBROUTINE INTEGRATE( TIN, TOUT, &
  47     ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
  48
  49     USE KPP_ROOT_Parameters, ONLY: NVAR
  50     USE KPP_ROOT_Global,     ONLY: ATOL,RTOL,VAR
  51
  52     IMPLICIT NONE
  53
  54     KPP_REAL :: TIN  ! TIN - Start Time
  55     KPP_REAL :: TOUT ! TOUT - End Time
  56     INTEGER,  INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
  57     KPP_REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
  58     INTEGER,  INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
  59     KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
  60     INTEGER,  INTENT(OUT), OPTIONAL :: IERR_U
  61
  62     KPP_REAL, SAVE :: H
  63     INTEGER :: IERR
  64
  65     KPP_REAL :: RCNTRL(20), RSTATUS(20)
  66     INTEGER :: ICNTRL(20), ISTATUS(20)
  67     INTEGER, SAVE :: Ntotal = 0
  68
  69     H =0.0_dp
  70
  71     !~~~> fine-tune the integrator:
  72     ICNTRL(:)  = 0
  73     ICNTRL(2)  = 0 ! 0=vector tolerances, 1=scalar tolerances
  74     ICNTRL(5)  = 8 ! Max no. of Newton iterations
  75     ICNTRL(6)  = 1 ! Starting values for Newton are interpolated (0) or zero (1)
  76     ICNTRL(11) = 1 ! Gustaffson (1) or classic(2) controller
  77     RCNTRL(1:20) = 0._dp
  78
  79     !~~~> if optional parameters are given, and if they are >0,
  80     !     then use them to overwrite default settings
  81     IF (PRESENT(ICNTRL_U)) ICNTRL(1:20) = ICNTRL_U(1:20)
  82     IF (PRESENT(RCNTRL_U)) RCNTRL(1:20) = RCNTRL_U(1:20)
  83
  84
  85     CALL RADAU5(  NVAR,TIN,TOUT,VAR,H,                  &
  86                   RTOL,ATOL,                            &
  87                   RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR  )
  88
  89 !!$    Ntotal = Ntotal + Nstp
  90 !!$    PRINT*,'NSTEPS=',Nstp,' (',Ntotal,')'
  91
  92     Nfun = Nfun + ISTATUS(1)
  93     Njac = Njac + ISTATUS(2)
  94     Nstp = Nstp + ISTATUS(3)
  95     Nacc = Nacc + ISTATUS(4)
  96     Nrej = Nrej + ISTATUS(5)
  97     Ndec = Ndec + ISTATUS(6)
  98     Nsol = Nsol + ISTATUS(7)
  99     Nsng = Nsng + ISTATUS(8)
 100
 101     ! if optional parameters are given for output
 102     ! use them to store information in them
 103     IF (PRESENT(ISTATUS_U)) THEN
 104       ISTATUS_U(:) = 0
 105       ISTATUS_U(1) = Nfun ! function calls
 106       ISTATUS_U(2) = Njac ! jacobian calls
 107       ISTATUS_U(3) = Nstp ! steps
 108       ISTATUS_U(4) = Nacc ! accepted steps
 109       ISTATUS_U(5) = Nrej ! rejected steps (except at the beginning)
 110       ISTATUS_U(6) = Ndec ! LU decompositions
 111       ISTATUS_U(7) = Nsol ! forward/backward substitutions
 112     ENDIF
 113     IF (PRESENT(RSTATUS_U)) THEN
 114       RSTATUS_U(:) = 0.
 115       RSTATUS_U(1) = TOUT ! final time
 116     ENDIF
 117     IF (PRESENT(IERR_U)) IERR_U = IERR
 118
 119 ! mz_rs_20050716: IERR is returned to the user who then decides what to do
 120 ! about it, i.e. either stop the run or ignore it.
 121 !!$    IF (IERR < 0) THEN
 122 !!$      PRINT *,'RADAU: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')'
 123 !!$      STOP
 124 !!$    ENDIF
 125
 126   END SUBROUTINE INTEGRATE
 127
 128    SUBROUTINE RADAU5(N,T,Tend,Y,H,RelTol,AbsTol,    &
 129                        RCNTRL,ICNTRL,RSTATUS,ISTATUS,IDID)
 130
 131 !~~~>-----------------------------------------------
 132 !     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
 133 !     SYSTEM OF FirstStep 0RDER ORDINARY DIFFERENTIAL EQUATIONS
 134 !                     M*Y'=F(T,Y).
 135 !     THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M  /=  I)
 136 !     OR EXPLICIT (M=I).
 137 !     THE METHOD USED IS AN IMPLICIT RUNGE-KUTTA METHOD (RADAU IIA)
 138 !     OF ORDER 5 WITH STEP SIZE CONTROL AND CONTINUOUS OUTPUT.
 139 !     C.F. SECTION IV.8
 140 !
 141 !     AUTHORS: E. HAIRER AND G. WANNER
 142 !              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
 143 !              CH-1211 GENEVE 24, SWITZERLAND
 144 !              E-MAIL:  HAIRER@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
 145 !
 146 !     THIS CODE IS PART OF THE BOOK:
 147 !         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
 148 !         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
 149 !         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
 150 !         SPRINGER-VERLAG (1991)
 151 !
 152 !     VERSION OF SEPTEMBER 30, 1995
 153 !
 154 !     INPUT PARAMETERS
 155 !     ----------------
 156 !     N           DIMENSION OF THE SYSTEM
 157 !
 158 !     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
 159 !                 VALUE OF F(T,Y):
 160 !                 SUBROUTINE FCN(N,T,Y,F)
 161 !                    KPP_REAL T,Y(N),F(N)
 162 !                    F(1)=...   ETC.
 163 !                 RPAR, IPAR (SEE BELOW)
 164 !
 165 !     T           INITIAL TIME VALUE
 166 !
 167 !     Tend        FINAL T-VALUE (Tend-T MAY BE POSITIVE OR NEGATIVE)
 168 !
 169 !     Y(N)        INITIAL VALUES FOR Y
 170 !
 171 !     H           INITIALL STEP SIZE GUESS;
 172 !                 FOR STIFF EQUATIONS WITH INITIALL TRANSIENT,
 173 !                 H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD.
 174 !                 THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS
 175 !                 QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6).
 176 !
 177 !     RelTol,AbsTol   RELATIVE AND ABSOLUTE ERROR TOLERANCES.
 178 !          for ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
 179 !                        = 1: AbsTol, RelTol are scalars
 180 !
 181 !          -----  CONTINUOUS OUTPUT: -----
 182 !                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
 183 !                 FOR THE INTERVAL [Told,T] IS AVAILABLE THROUGH
 184 !                 THE FUNCTION
 185 !                        >>>   CONTR5(I,S,CONT,LRC)   <<<
 186 !                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH
 187 !                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
 188 !                 S SHOULD LIE IN THE INTERVAL [Told,T].
 189 !                 DO NOT CHANGE THE ENTRIES OF CONT(LRC), IF THE
 190 !                 DENSE OUTPUT FUNCTION IS USED.
 191 !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 192 !
 193 !~~~>     INPUT PARAMETERS:
 194 !
 195 !    Note: For input parameters equal to zero the default values of the
 196 !          corresponding variables are used.
 197 !
 198 !~~~>  Integer input parameters:
 199 !
 200 !    ICNTRL(1) = not used
 201 !
 202 !    ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
 203 !              = 1: AbsTol, RelTol are scalars
 204 !
 205 !    ICNTRL(3) = not used
 206 !
 207 !    ICNTRL(4)  -> maximum number of integration steps
 208 !        For ICNTRL(4)=0 the default value of 10000 is used
 209 !
 210 !    ICNTRL(5)  -> maximum number of Newton iterations
 211 !        For ICNTRL(5)=0 the default value of 8 is used
 212 !
 213 !    ICNTRL(6)  -> starting values of Newton iterations:
 214 !        ICNTRL(6)=0 : starting values are obtained from
 215 !                      the extrapolated collocation solution
 216 !                      (the default)
 217 !        ICNTRL(6)=1 : starting values are zero
 218 !
 219 !    ICNTRL(11) -> switch for step size strategy
 220 !              ICNTRL(8) == 1:  mod. predictive controller (Gustafsson)
 221 !              ICNTRL(8) == 2:  classical step size control
 222 !              the default value (for iwork(8)=0) is iwork(8)=1.
 223 !              the choice iwork(8) == 1 seems to produce safer results;
 224 !              for simple problems, the choice iwork(8) == 2 produces
 225 !              often slightly faster runs
 226 !              ( currently unused )
 227 !
 228 !~~~>  Real input parameters:
 229 !
 230 !    RCNTRL(1)  -> not used
 231 !
 232 !    RCNTRL(2)  -> Hmax, upper bound for the integration step size
 233 !
 234 !    RCNTRL(3)  -> not used
 235 !
 236 !    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
 237 !
 238 !    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
 239 !
 240 !    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
 241 !                 (default=0.1)
 242 !
 243 !    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
 244 !                  than the predicted value  (default=0.9)
 245 !
 246 !    RCNTRL(8)  -> ThetaMin. If Newton convergence rate smaller
 247 !                  than ThetaMin the Jacobian is not recomputed;
 248 !                  (default=0.001)
 249 !
 250 !    RCNTRL(9)  -> NewtonTol, stopping criterion for Newton's method
 251 !                  (default=0.03)
 252 !
 253 !    RCNTRL(10) -> Qmin
 254 !
 255 !    RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
 256 !                  step size is kept constant and the LU factorization
 257 !                  reused (default Qmin=1, Qmax=1.2)
 258 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 259 !
 260 !
 261 !    OUTPUT PARAMETERS
 262 !    -----------------
 263 !    T           T-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED
 264 !                     (AFTER SUCCESSFUL RETURN T=Tend).
 265 !
 266 !    Y(N)        NUMERICAL SOLUTION AT T
 267 !
 268 !    H           PREDICTED STEP SIZE OF THE LastStep ACCEPTED STEP
 269 !
 270 !    IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
 271 !
 272 !    ISTATUS(1) = No. of function calls
 273 !    ISTATUS(2) = No. of jacobian calls
 274 !    ISTATUS(3) = No. of steps
 275 !    ISTATUS(4) = No. of accepted steps
 276 !    ISTATUS(5) = No. of rejected steps (except at the beginning)
 277 !    ISTATUS(6) = No. of LU decompositions
 278 !    ISTATUS(7) = No. of forward/backward substitutions
 279 !    ISTATUS(8) = No. of singular matrix decompositions
 280 !
 281 !    RSTATUS(1)  -> Texit, the time corresponding to the
 282 !                   computed Y upon return
 283 !    RSTATUS(2)  -> Hexit, last accepted step before exit
 284 !                   For multiple restarts, use Hexit as Hstart
 285 !                   in the subsequent run
 286 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 287
 288       IMPLICIT NONE
 289
 290       INTEGER :: N
 291       KPP_REAL :: Y(N),AbsTol(N),RelTol(N),RCNTRL(20),RSTATUS(20)
 292       INTEGER :: ICNTRL(20), ISTATUS(20)
 293       LOGICAL :: StartNewton, Gustafsson
 294       INTEGER :: IDID, ITOL
 295       KPP_REAL :: H,Tend,T
 296
 297       !~~~> Control arguments
 298       INTEGER :: Max_no_steps, NewtonMaxit
 299       KPP_REAL :: Hstart,Hmin,Hmax,Qmin,Qmax
 300       KPP_REAL :: Roundoff, ThetaMin,TolNewton
 301       KPP_REAL :: FacSafe,FacMin,FacMax,FacRej
 302       !~~~> Local variables
 303       INTEGER :: i
 304       KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
 305
 306       !~~~> variables from the former COMMON block /CONRA5/
 307       ! KPP_REAL :: Tsol, Hsol
 308
 309 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 310 !        SETTING THE PARAMETERS
 311 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 312        Nfun=0
 313        Njac=0
 314        Nstp=0
 315        Nacc=0
 316        Nrej=0
 317        Ndec=0
 318        Nsol=0
 319        IDID = 0
 320
 321 !~~~> ICNTRL(1) - autonomous system - not used
 322 !~~~> ITOL: 1 for vector and 0 for scalar AbsTol/RelTol
 323       IF (ICNTRL(2) == 0) THEN
 324          ITOL = 1
 325       ELSE
 326          ITOL = 0
 327       END IF
 328 !~~~> ICNTRL(3) - method selection - not used
 329 !~~~> Max_no_steps: the maximal number of time steps
 330       IF (ICNTRL(4) == 0) THEN
 331          Max_no_steps = 10000
 332       ELSE
 333          Max_no_steps=ICNTRL(4)
 334          IF (Max_no_steps <= 0) THEN
 335             WRITE(6,*) 'ICNTRL(4)=',ICNTRL(4)
 336             CALL RAD_ErrorMsg(-1,T,ZERO,IDID)
 337          END IF
 338       END IF
 339 !~~~> NewtonMaxit    MAXIMAL NUMBER OF NEWTON ITERATIONS
 340       IF (ICNTRL(5) == 0) THEN
 341          NewtonMaxit = 8
 342       ELSE
 343          NewtonMaxit=ICNTRL(5)
 344          IF (NewtonMaxit <= 0) THEN
 345             WRITE(6,*) 'ICNTRL(5)=',ICNTRL(5)
 346             CALL RAD_ErrorMsg(-2,T,ZERO,IDID)
 347           END IF
 348       END IF
 349 !~~~> StartNewton:  Use extrapolation for starting values of Newton iterations
 350       IF (ICNTRL(6) == 0) THEN
 351          StartNewton = .TRUE.
 352       ELSE
 353          StartNewton = .FALSE.
 354       END IF
 355 !~~~> Gustafsson: step size controller
 356       IF(ICNTRL(11) == 0)THEN
 357          Gustafsson=.TRUE.
 358       ELSE
 359          Gustafsson=.FALSE.
 360       END IF
 361
 362
 363 !~~~> Roundoff   SMALLEST NUMBER SATISFYING 1.0d0+Roundoff>1.0d0
 364       Roundoff=WLAMCH('E');
 365
 366 !~~~> RCNTRL(1) = Hmin - not used
 367       Hmin = ZERO
 368 !~~~> Hmax = maximal step size
 369       IF (RCNTRL(2) == ZERO) THEN
 370          Hmax=Tend-T
 371       ELSE
 372          Hmax=MAX(ABS(RCNTRL(7)),ABS(Tend-T))
 373       END IF
 374 !~~~> RCNTRL(3) = Hstart - not used
 375       Hstart = ZERO
 376 !~~~> FacMin: lower bound on step decrease factor
 377       IF(RCNTRL(4) == ZERO)THEN
 378          FacMin = 0.2d0
 379       ELSE
 380          FacMin = RCNTRL(4)
 381       END IF
 382 !~~~> FacMax: upper bound on step increase factor
 383       IF(RCNTRL(5) == ZERO)THEN
 384          FacMax = 8.D0
 385       ELSE
 386          FacMax = RCNTRL(5)
 387       END IF
 388 !~~~> FacRej: step decrease factor after 2 consecutive rejections
 389       IF(RCNTRL(6) == ZERO)THEN
 390          FacRej = 0.1d0
 391       ELSE
 392          FacRej = RCNTRL(6)
 393       END IF
 394 !~~~> FacSafe:  by which the new step is slightly smaller
 395 !               than the predicted value
 396       IF (RCNTRL(7) == ZERO) THEN
 397          FacSafe=0.9d0
 398       ELSE
 399          FacSafe=RCNTRL(7)
 400       END IF
 401       IF ( (FacMax < ONE) .OR. (FacMin > ONE) .OR. &
 402            (FacSafe <= 0.001D0) .OR. (FacSafe >= 1.0d0) ) THEN
 403             WRITE(6,*)'RCNTRL(4:7)=',RCNTRL(4:7)
 404             CALL RAD_ErrorMsg(-4,T,ZERO,IDID)
 405       END IF
 406
 407 !~~~> ThetaMin:  decides whether the Jacobian should be recomputed
 408       IF (RCNTRL(8) == ZERO) THEN
 409          ThetaMin = 1.0d-3
 410       ELSE
 411          ThetaMin=RCNTRL(8)
 412          IF (ThetaMin <= 0.0d0 .OR. ThetaMin >= 1.0d0) THEN
 413             WRITE(6,*) 'RCNTRL(8)=', RCNTRL(8)
 414             CALL RAD_ErrorMsg(-5,T,ZERO,IDID)
 415          END IF
 416       END IF
 417 !~~~> TolNewton:  stopping crierion for Newton's method
 418       IF (RCNTRL(9) == ZERO) THEN
 419          TolNewton = 3.0d-2
 420       ELSE
 421          TolNewton = RCNTRL(9)
 422          IF (TolNewton <= Roundoff) THEN
 423             WRITE(6,*) 'RCNTRL(9)=',RCNTRL(9)
 424             CALL RAD_ErrorMsg(-6,T,ZERO,IDID)
 425          END IF
 426       END IF
 427 !~~~> Qmin AND Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST.
 428       IF (RCNTRL(10) == ZERO) THEN
 429          Qmin=1.D0
 430       ELSE
 431          Qmin=RCNTRL(10)
 432       END IF
 433       IF (RCNTRL(11) == ZERO) THEN
 434          Qmax=1.2D0
 435       ELSE
 436          Qmax=RCNTRL(11)
 437       END IF
 438       IF (Qmin > ONE .OR. Qmax < ONE) THEN
 439          WRITE(6,*) 'RCNTRL(10:11)=',Qmin,Qmax
 440          CALL RAD_ErrorMsg(-7,T,ZERO,IDID)
 441       END IF
 442 !~~~> Check if tolerances are reasonable
 443       IF (ITOL == 0) THEN
 444           IF (AbsTol(1) <= ZERO.OR.RelTol(1) <= 10.d0*Roundoff) THEN
 445               WRITE (6,*) 'AbsTol/RelTol=',AbsTol,RelTol
 446               CALL RAD_ErrorMsg(-8,T,ZERO,IDID)
 447           END IF
 448       ELSE
 449           DO i=1,N
 450           IF (AbsTol(i) <= ZERO.OR.RelTol(i) <= 10.d0*Roundoff) THEN
 451               WRITE (6,*) 'AbsTol/RelTol(',i,')=',AbsTol(i),RelTol(i)
 452               CALL RAD_ErrorMsg(-8,T,ZERO,IDID)
 453           END IF
 454           END DO
 455       END IF
 456
 457 !~~~> WHEN A FAIL HAS OCCURED, RETURN
 458       IF (IDID < 0) RETURN
 459
 460
 461 !~~~>  CALL TO CORE INTEGRATOR ------------
 462       CALL RAD_Integrator( N,T,Y,Tend,Hmax,H,RelTol,AbsTol,ITOL,IDID, &
 463         Max_no_steps,Roundoff,FacSafe,ThetaMin,TolNewton,Qmin,Qmax,   &
 464         NewtonMaxit,StartNewton,Gustafsson,FacMin,FacMax,FacRej )
 465
 466       ISTATUS(1)=Nfun
 467       ISTATUS(2)=Njac
 468       ISTATUS(3)=Nstp
 469       ISTATUS(4)=Nacc
 470       ISTATUS(5)=Nrej
 471       ISTATUS(6)=Ndec
 472       ISTATUS(7)=Nsol
 473       ISTATUS(8)=Nsng
 474
 475    ! END SUBROUTINE RADAU5
 476    CONTAINS ! INTERNAL PROCEDURES TO RADAU5
 477
 478
 479 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 480    SUBROUTINE RAD_Integrator( N,T,Y,Tend,Hmax,H,RelTol,AbsTol,ITOL,IDID, &
 481         Max_no_steps,Roundoff,FacSafe,ThetaMin,TolNewton,Qmin,Qmax,      &
 482         NewtonMaxit,StartNewton,Gustafsson,FacMin,FacMax,FacRej )
 483 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 484 !     CORE INTEGRATOR FOR RADAU5
 485 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 486
 487       IMPLICIT NONE
 488       INTEGER :: N
 489       KPP_REAL  Y(NVAR),Z1(NVAR),Z2(NVAR),Z3(NVAR),Y0(NVAR),&
 490                      SCAL(NVAR),F1(NVAR),F2(NVAR),F3(NVAR),      &
 491                      CONT(N,4),AbsTol(NVAR),RelTol(NVAR)
 492
 493 #ifdef FULL_ALGEBRA
 494       KPP_REAL  :: FJAC(NVAR,NVAR), E1(NVAR,NVAR)
 495       DOUBLE COMPLEX :: E2(NVAR,NVAR)
 496 #else
 497       KPP_REAL  :: FJAC(LU_NONZERO), E1(LU_NONZERO)
 498       DOUBLE COMPLEX :: E2(LU_NONZERO)
 499 #endif
 500
 501       !~~~> Local variables
 502       KPP_REAL  :: TMP(NVAR), T, Tend, Tdirection, &
 503                  H, Hmax, HmaxN, Hacc, Hnew, Hopt, Hold, &
 504                  Fac, FacMin, Facmax, FacSafe, FacRej, FacGus, FacConv, &
 505                  Theta, ThetaMin, TolNewton, ERR, ERRACC,   &
 506                  Qmin, Qmax, DYNO, Roundoff, &
 507                  AK, AK1, AK2, AK3, C3Q, &
 508                  Qnewton, DYTH, THQ, THQOLD, DYNOLD, &
 509                  DENOM, C1Q, C2Q, ALPHA, BETA, GAMMA, CFAC, ACONT3, QT
 510       INTEGER :: IP1(NVAR),IP2(NVAR), ITOL, IDID, Max_no_steps, &
 511                  NewtonIter, NewtonMaxit, ISING
 512       LOGICAL :: REJECT, FirstStep, FreshJac, LastStep, &
 513                  Gustafsson, StartNewton, NewtonDone
 514       ! KPP_REAL, PARAMETER  :: ONE = 1.0d0, ZERO = 0.0d0
 515
 516 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 517 !  INITIALISATIONS
 518 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 519
 520       CALL RAD_Coefficients
 521
 522       Tdirection=SIGN(1.D0,Tend-T)
 523       HmaxN=MIN(ABS(Hmax),ABS(Tend-T))
 524       H=MIN(ABS(Hmin),ABS(Hstart))
 525       H=MIN(ABS(H),HmaxN)
 526       IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6
 527       H=SIGN(H,Tdirection)
 528       Hold=H
 529       REJECT=.FALSE.
 530       FirstStep=.TRUE.
 531       LastStep=.FALSE.
 532       FreshJac=.FALSE.; Theta=1.0d0
 533       IF ((T+H*1.0001D0-Tend)*Tdirection >= 0.D0) THEN
 534          H=Tend-T
 535          LastStep=.TRUE.
 536       END IF
 537       FacConv=1.D0
 538       CFAC=FacSafe*(1+2*NewtonMaxit)
 539       Nsng=0
 540 !      Told=T
 541       CALL RAD_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
 542       CALL FUN_CHEM(T,Y,Y0)
 543
 544 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 545 !~~~>  Time loop begins
 546 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 547 Tloop: DO WHILE ( (Tend-T)*Tdirection - Roundoff > ZERO )
 548
 549       !~~~> COMPUTE JACOBIAN MATRIX ANALYTICALLY
 550       IF ( (.NOT.FreshJac)  .AND. (Theta > ThetaMin) ) THEN
 551         CALL JAC_CHEM(T,Y,FJAC)
 552         FreshJac=.TRUE.
 553       END IF
 554
 555       !~~~> Compute the matrices E1 and E2 and their decompositions
 556       GAMMA  = rkGamma/H
 557       ALPHA = rkAlpha/H
 558       BETA  = rkBeta/H
 559       CALL RAD_DecompReal(N,FJAC,GAMMA,E1,IP1,ISING)
 560       IF (ISING /= 0) THEN
 561           Nsng=Nsng+1
 562           IF (Nsng >= 5) THEN
 563             CALL RAD_ErrorMsg(-12,T,H,IDID); RETURN
 564           END IF
 565           H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE.
 566           CYCLE Tloop
 567       END IF
 568       CALL RAD_DecompCmplx(N,FJAC,ALPHA,BETA,E2,IP2,ISING)
 569       IF (ISING /= 0) THEN
 570           Nsng=Nsng+1
 571           IF (Nsng >= 5) THEN
 572             CALL RAD_ErrorMsg(-12,T,H,IDID); RETURN
 573           END IF
 574           H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE.
 575           CYCLE Tloop
 576       END IF
 577
 578    30 CONTINUE
 579       Nstp=Nstp+1
 580       IF (Nstp > Max_no_steps) THEN
 581         PRINT*,'Max number of time steps is ',Max_no_steps
 582         CALL RAD_ErrorMsg(-9,T,H,IDID); RETURN
 583       END IF
 584       IF (0.1D0*ABS(H) <= ABS(T)*Roundoff)  THEN
 585         CALL RAD_ErrorMsg(-10,T,H,IDID); RETURN
 586       END IF
 587
 588 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 589 !  STARTING VALUES FOR NEWTON ITERATION
 590 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 591       IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN
 592          CALL Set2zero(N,Z1)
 593          CALL Set2zero(N,Z2)
 594          CALL Set2zero(N,Z3)
 595          CALL Set2zero(N,F1)
 596          CALL Set2zero(N,F2)
 597          CALL Set2zero(N,F3)
 598       ELSE
 599          C3Q=H/Hold
 600          C1Q=rkC(1)*C3Q
 601          C2Q=rkC(2)*C3Q
 602          DO i=1,N
 603             AK1=CONT(i,2)
 604             AK2=CONT(i,3)
 605             AK3=CONT(i,4)
 606             Z1(i)=C1Q*(AK1+(C1Q-rkC(2)+ONE)*(AK2+(C1Q-rkC(1)+ONE)*AK3))
 607             Z2(i)=C2Q*(AK1+(C2Q-rkC(2)+ONE)*(AK2+(C2Q-rkC(1)+ONE)*AK3))
 608             Z3(i)=C3Q*(AK1+(C3Q-rkC(2)+ONE)*(AK2+(C3Q-rkC(1)+ONE)*AK3))
 609          END DO
 610          !  F(1,2,3) = TransfInv x Z(1,2,3)
 611          CALL RAD_Transform(N,TransfInv,Z1,Z2,Z3,F1,F2,F3)
 612       END IF
 613 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 614 !  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
 615 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 616
 617             FacConv = MAX(FacConv,Roundoff)**0.8D0
 618             Theta=ABS(ThetaMin)
 619
 620 NewtonLoop:DO  NewtonIter = 1, NewtonMaxit
 621
 622             !~~~>  The right-hand side
 623             DO i=1,N
 624                TMP(i)=Y(i)+Z1(i)
 625             END DO
 626             CALL FUN_CHEM(T+rkC(1)*H,TMP,Z1)
 627             DO i=1,N
 628                TMP(i)=Y(i)+Z2(i)
 629             END DO
 630             CALL FUN_CHEM(T+rkC(2)*H,TMP,Z2)
 631             DO i=1,N
 632                TMP(i)=Y(i)+Z3(i)
 633             END DO
 634             CALL FUN_CHEM(T+rkC(3)*H,TMP,Z3)
 635
 636             !~~~> Solve the linear systems
 637             !  Z(1,2,3) = TransfInv x Z(1,2,3)
 638             CALL RAD_Transform(N,TransfInv,Z1,Z2,Z3,Z1,Z2,Z3)
 639             CALL RAD_Solve( N,FJAC,GAMMA,ALPHA,BETA,E1,E2, &
 640                         Z1,Z2,Z3,F1,F2,F3,CONT,IP1,IP2,ISING )
 641             Nsol=Nsol+1
 642
 643             DYNO=0.0d0
 644             DO i=1,N
 645                DENOM=SCAL(i)
 646                DYNO=DYNO+(Z1(i)/DENOM)**2+(Z2(i)/DENOM)**2+(Z3(i)/DENOM)**2
 647             END DO
 648             DYNO=SQRT(DYNO/(3*N))
 649
 650             !~~~> Bad convergence or number of iterations too large
 651             IF ( (NewtonIter > 1) .AND. (NewtonIter < NewtonMaxit) ) THEN
 652                 THQ=DYNO/DYNOLD
 653                 IF (NewtonIter == 2) THEN
 654                    Theta=THQ
 655                 ELSE
 656                    Theta=SQRT(THQ*THQOLD)
 657                 END IF
 658                 THQOLD=THQ
 659                 IF (Theta < 0.99d0) THEN
 660                     FacConv=Theta/(1.0d0-Theta)
 661                     DYTH=FacConv*DYNO*Theta**(NewtonMaxit-1-NewtonIter)/TolNewton
 662                     IF (DYTH >= 1.0d0) THEN
 663                          Qnewton=MAX(1.0D-4,MIN(20.0d0,DYTH))
 664                          FAC=.8D0*Qnewton**(-1.0d0/(4.0d0+NewtonMaxit-1-NewtonIter))
 665                          H=FAC*H
 666                          REJECT=.TRUE.
 667                          LastStep=.FALSE.
 668                          CYCLE Tloop
 669                     END IF
 670                 ELSE ! Non-convergence of Newton
 671                    H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE.
 672                    CYCLE Tloop
 673                 END IF
 674             END IF
 675             DYNOLD=MAX(DYNO,Roundoff)
 676             CALL WAXPY(N,ONE,Z1,1,F1,1) ! F1 <- F1 + Z1
 677             CALL WAXPY(N,ONE,Z2,1,F2,1) ! F2 <- F2 + Z2
 678             CALL WAXPY(N,ONE,Z3,1,F3,1) ! F3 <- F3 + Z3
 679             !  Z(1,2,3) = Transf x F(1,2,3)
 680             CALL RAD_Transform(N,Transf,F1,F2,F3,Z1,Z2,Z3)
 681             NewtonDone = (FacConv*DYNO <= TolNewton)
 682             IF (NewtonDone) EXIT NewtonLoop
 683
 684        END DO NewtonLoop
 685
 686        IF (.NOT.NewtonDone) THEN
 687            CALL RAD_ErrorMsg(-8,T,H,IDID);
 688            H=H*0.5D0; REJECT=.TRUE.; LastStep=.FALSE.
 689            CYCLE Tloop
 690        END IF
 691
 692
 693 !~~~> ERROR ESTIMATION
 694       CALL  RAD_ErrorEstimate(N,FJAC,H,Y0,Y,T, &
 695                E1,Z1,Z2,Z3,IP1,SCAL,ERR,       &
 696                FirstStep,REJECT,GAMMA)
 697 !~~~> COMPUTATION OF Hnew
 698       Fac  = ERR**(-0.25d0)*   &
 699              MIN(FacSafe,(NewtonIter+2*NewtonMaxit)/CFAC)
 700       Fac  = MIN(FacMax,MAX(FacMin,Fac))
 701       Hnew = Fac*H
 702
 703 !~~~> IS THE ERROR SMALL ENOUGH ?
 704 accept:IF (ERR < ONE) THEN !~~~> STEP IS ACCEPTED
 705          FirstStep=.FALSE.
 706          Nacc=Nacc+1
 707          IF (Gustafsson) THEN
 708             !~~~> Predictive controller of Gustafsson
 709             !~~~> Currently not implemented
 710             IF (Nacc > 1) THEN
 711                FacGus=FacSafe*(H/Hacc)*(ERR**2/ERRACC)**(-0.25d0)
 712                FacGus=MIN(FacMax,MAX(FacMin,FacGus))
 713                Fac=MIN(Fac,FacGus)
 714                Hnew=H*Fac
 715             END IF
 716             Hacc=H
 717             ERRACC=MAX(1.0D-2,ERR)
 718          END IF
 719          ! Told = T
 720          Hold = H
 721          T=T+H
 722          DO i=1,N
 723             Y(i)=Y(i)+Z3(i)
 724             CONT(i,2)=(Z2(i)-Z3(i))/(rkC(2)-ONE)
 725             AK=(Z1(i)-Z2(i))/(rkC(1)-rkC(2))
 726             ACONT3=Z1(i)/rkC(1)
 727             ACONT3=(AK-ACONT3)/rkC(2)
 728             CONT(i,3)=(AK-CONT(i,2))/(rkC(1)-ONE)
 729             CONT(i,4)=CONT(i,3)-ACONT3
 730          END DO
 731          CALL RAD_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
 732          FreshJac=.FALSE.
 733          IF (LastStep) THEN
 734             H=Hopt
 735             IDID=1
 736             RETURN
 737          END IF
 738          CALL FUN_CHEM(T,Y,Y0)
 739          Hnew=Tdirection*MIN(ABS(Hnew),HmaxN)
 740          Hopt=Hnew
 741          Hopt=MIN(H,Hnew)
 742          IF (REJECT) Hnew=Tdirection*MIN(ABS(Hnew),ABS(H))
 743          REJECT=.FALSE.
 744          IF ((T+Hnew/Qmin-Tend)*Tdirection >= 0.D0) THEN
 745             H=Tend-T
 746             LastStep=.TRUE.
 747          ELSE
 748             QT=Hnew/H
 749             IF ( (Theta<=ThetaMin) .AND. (QT>=Qmin) &
 750                             .AND. (QT<=Qmax) ) GOTO 30
 751             H=Hnew
 752          END IF
 753          CYCLE Tloop
 754       ELSE accept !~~~> STEP IS REJECTED
 755          REJECT=.TRUE.
 756          LastStep=.FALSE.
 757          IF (FirstStep) THEN
 758              H=H*FacRej
 759          ELSE
 760              H=Hnew
 761          END IF
 762          IF (Nacc >= 1) Nrej=Nrej+1
 763          CYCLE Tloop
 764       END IF accept
 765
 766
 767       END DO Tloop
 768
 769
 770    END SUBROUTINE RAD_Integrator
 771
 772 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 773  SUBROUTINE RAD_ErrorMsg(Code,T,H,IERR)
 774 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 775 !    Handles all error messages
 776 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 777
 778    IMPLICIT NONE
 779    KPP_REAL, INTENT(IN) :: T, H
 780    INTEGER, INTENT(IN)  :: Code
 781    INTEGER, INTENT(OUT) :: IERR
 782
 783    IERR = Code
 784    PRINT * , &
 785      'Forced exit from RADAU5 due to the following error:'
 786    IF ((Code>=-11).AND.(Code<=-1)) THEN
 787      PRINT *, IERR_NAMES(Code)
 788    ELSE
 789      PRINT *, 'Unknown Error code: ', Code
 790    ENDIF
 791
 792    PRINT *, "T=", T, "and H=", H
 793
 794  END SUBROUTINE RAD_ErrorMsg
 795
 796
 797 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 798  SUBROUTINE RAD_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
 799 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 800 !    Handles all error messages
 801 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 802    IMPLICIT NONE
 803    INTEGER, INTENT(IN)  :: N, ITOL
 804    KPP_REAL, INTENT(IN) :: AbsTol(*), RelTol(*), Y(N)
 805    KPP_REAL, INTENT(OUT) :: SCAL(N)
 806    INTEGER :: i
 807
 808    IF (ITOL==0) THEN
 809        DO i=1,N
 810           SCAL(i)=AbsTol(1)+RelTol(1)*ABS(Y(i))
 811        END DO
 812    ELSE
 813        DO i=1,N
 814           SCAL(i)=AbsTol(i)+RelTol(i)*ABS(Y(i))
 815        END DO
 816    END IF
 817
 818  END SUBROUTINE RAD_ErrorScale
 819
 820 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 821   SUBROUTINE RAD_Transform(N,Tr,Z1,Z2,Z3,F1,F2,F3)
 822 !~~~>                 F = Tr x Z
 823 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 824       IMPLICIT NONE
 825       INTEGER :: N, i
 826       KPP_REAL :: Tr(3,3),Z1(N),Z2(N),Z3(N),F1(N),F2(N),F3(N)
 827       KPP_REAL :: x1, x2, x3
 828       DO i=1,N
 829           x1 = Z1(i); x2 = Z2(i); x3 = Z3(i)
 830           F1(i) = Tr(1,1)*x1 + Tr(1,2)*x2 + Tr(1,3)*x3
 831           F2(i) = Tr(2,1)*x1 + Tr(2,2)*x2 + Tr(2,3)*x3
 832           F3(i) = Tr(3,1)*x1 + Tr(3,2)*x2 + Tr(3,3)*x3
 833       END DO
 834   END SUBROUTINE RAD_Transform
 835
 836 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 837   SUBROUTINE RAD_Coefficients
 838 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 839       IMPLICIT NONE
 840       KPP_REAL :: s2,s3,s6,x1,x2,x3,x4,y1,y2,y3,y4,y5
 841
 842       s2 = SQRT(2.0d0);
 843       s3 = SQRT(3.0d0);
 844       s6 = SQRT(6.0d0);
 845       x1 = 3.d0**(1.d0/3.d0);
 846       x2 = 3.d0**(2.d0/3.d0);
 847       x3 = 3.d0**(1.d0/6.d0);
 848       x4 = 3.d0**(5.d0/6.d0);
 849
 850       rkA(1,1) =  11.d0/45.d0-7.d0/360.d0*s6
 851       rkA(1,2) =  37.d0/225.d0-169.d0/1800.d0*s6
 852       rkA(1,3) =  -2.d0/225.d0+s6/75
 853       rkA(2,1) =  37.d0/225.d0+169.d0/1800.d0*s6
 854       rkA(2,2) =  11.d0/45.d0+7.d0/360.d0*s6
 855       rkA(2,3) =  -2.d0/225.d0-s6/75
 856       rkA(3,1) =  4.d0/9.d0-s6/36
 857       rkA(3,2) =  4.d0/9.d0+s6/36
 858       rkA(3,3) =  1.d0/9.d0
 859
 860       rkB(1)   =  4.d0/9.d0-s6/36
 861       rkB(2)   =  4.d0/9.d0+s6/36
 862       rkB(3)   =  1.d0/9.d0
 863
 864       rkC(1)   =  2.d0/5.d0-s6/10
 865       rkC(2)   =  2.d0/5.d0+s6/10
 866       rkC(3)   =  1.d0
 867
 868       ! Error estimation
 869       rkE(1) = -(13.d0+7.d0*s6)/3.d0
 870       rkE(2) =  (-13.d0+7.d0*s6)/3.d0
 871       rkE(3) = -1.d0/3.d0
 872
 873       !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 874       !~~~> Diagonalize the RK matrix:
 875       ! TransfInv * inv(rkA) * Transf =
 876       !       |  rkGamma      0           0       |
 877       !       |      0      rkAlpha   -rkBeta   |
 878       !       |      0      rkBeta     rkAlpha  |
 879       !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 880       rkGamma =  3-x1+x2
 881       rkAlpha =  x1/2-x2/2+3
 882       rkBeta  =  -x4/2-3.d0/2.d0*x3
 883
 884       y1 = 36.d0/625.d0*s6
 885       y2 = 129.d0/2500.d0*x1
 886       y3 = 111.d0/2500.d0*x3*s2
 887       Transf(1,1) = -31.d0/1250.d0*s6*x1+37.d0/1250.d0*s6*x2-y1 &
 888                     +129.d0/1250.d0*x1-33.d0/1250.d0*x2+49.d0/625.d0
 889       Transf(1,2) =  -y1-y2-y3 &
 890                     +31.d0/2500.d0*x4*s2+33.d0/2500.d0*x2+49.d0/625.d0
 891       Transf(1,3) =  3.d0/2500.d0*x3*(-33-43*x2+31*x3*s2+37*s3*s2)
 892       Transf(2,1) =  31.d0/1250.d0*s6*x1-37.d0/1250.d0*s6*x2+y1 &
 893                     +129.d0/1250.d0*x1-33.d0/1250.d0*x2+49.d0/625.d0
 894       Transf(2,2) =  y1-y2+y3&
 895                     -31.d0/2500.d0*x4*s2+33.d0/2500.d0*x2+49.d0/625.d0
 896       Transf(2,3) =  -3.d0/2500.d0*x3*(33+43*x2+31*x3*s2+37*s3*s2)
 897       Transf(3,1) =  1.d0
 898       Transf(3,2) =  1.d0
 899       Transf(3,3) =  0.d0
 900
 901       y1 = 11.d0/36.d0*x3*s2 + 43.d0/108.d0*x4*s2
 902       y2 = 11.d0/36.d0*s2*x2 - 43.d0/36.d0*s2*x1
 903       y3 = 31.d0/54.d0*x1 + 37.d0/54.d0*x2
 904       y4 = 31.d0/54.d0*x4-37.d0/18.d0*x3
 905       y5 = -x2/27+5.d0/27.d0*x1
 906       TransfInv(1,1) =  y1 + y3
 907       TransfInv(1,2) = -y1 + y3
 908       TransfInv(1,3) =  y5 + 1.d0/3.d0
 909       TransfInv(2,1) = -y1 - y3
 910       TransfInv(2,2) =  y1 - y3
 911       TransfInv(2,3) = -y5 + 2.d0/3.d0
 912       TransfInv(3,1) =  y4 - y2
 913       TransfInv(3,2) =  y4 + y2
 914       TransfInv(3,3) =  x3/9+5.d0/27.d0*x4
 915
 916    END SUBROUTINE RAD_Coefficients
 917
 918
 919 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 920    SUBROUTINE RAD_DecompReal(N,FJAC,GAMMA,E1,IP1,ISING)
 921 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 922       IMPLICIT NONE
 923
 924       INTEGER :: N, ISING
 925       KPP_REAL :: GAMMA
 926 #ifdef FULL_ALGEBRA
 927       KPP_REAL :: FJAC(NVAR,NVAR),E1(NVAR,NVAR)
 928 #else
 929       KPP_REAL :: FJAC(LU_NONZERO),E1(LU_NONZERO)
 930 #endif
 931       INTEGER :: IP1(N), i, j
 932
 933 #ifdef FULL_ALGEBRA
 934       DO j=1,N
 935          DO  i=1,N
 936             E1(i,j)=-FJAC(i,j)
 937          END DO
 938          E1(j,j)=E1(j,j)+GAMMA
 939       END DO
 940       CALL DGETRF(N,N,E1,N,IP1,ISING)
 941 #else
 942       DO i=1,LU_NONZERO
 943          E1(i)=-FJAC(i)
 944       END DO
 945       DO i=1,NVAR
 946          j = LU_DIAG(i); E1(j)=E1(j)+GAMMA
 947       END DO
 948       CALL KppDecomp(E1,ISING)
 949 #endif
 950
 951    END SUBROUTINE RAD_DecompReal
 952
 953
 954 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 955    SUBROUTINE RAD_DecompCmplx(N,FJAC,ALPHA,BETA,E2,IP2,ISING)
 956 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 957       IMPLICIT NONE
 958       INTEGER :: N, ISING
 959 #ifdef FULL_ALGEBRA
 960       KPP_REAL ::  FJAC(N,N)
 961       DOUBLE COMPLEX :: E2(N,N)
 962 #else
 963       KPP_REAL ::  FJAC(LU_NONZERO)
 964       DOUBLE COMPLEX :: E2(LU_NONZERO)
 965 #endif
 966       KPP_REAL :: ALPHA, BETA
 967       INTEGER :: IP2(N), i, j
 968
 969 #ifdef FULL_ALGEBRA
 970       DO j=1,N
 971         DO i=1,N
 972           E2(i,j) = DCMPLX( -FJAC(i,j), 0.0d0 )
 973         END DO
 974         E2(j,j) = E2(j,j) + DCMPLX( ALPHA, BETA )
 975       END DO
 976       CALL ZGETRF(N,N,E2,N,IP2,ISING)
 977 #else
 978       DO i=1,LU_NONZERO
 979          E2(i) = DCMPLX( -FJAC(i), 0.0d0 )
 980       END DO
 981       DO i=1,NVAR
 982          j = LU_DIAG(i); E2(j)=E2(j)+DCMPLX( ALPHA, BETA )
 983       END DO
 984       CALL KppDecompCmplx(E2,ISING)
 985 #endif
 986       Ndec=Ndec+1
 987
 988    END SUBROUTINE RAD_DecompCmplx
 989
 990
 991 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 992    SUBROUTINE RAD_Solve(N,FJAC,GAMMA,ALPHA,BETA,E1,E2,&
 993                Z1,Z2,Z3,F1,F2,F3,CONT,IP1,IP2,ISING)
 994 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 995       IMPLICIT NONE
 996       INTEGER :: N,IP1(NVAR),IP2(NVAR),ISING
 997 #ifdef FULL_ALGEBRA
 998       KPP_REAL  :: FJAC(NVAR,NVAR), E1(NVAR,NVAR)
 999       DOUBLE COMPLEX :: E2(NVAR,NVAR)
1000 #else
1001       KPP_REAL  :: FJAC(LU_NONZERO), E1(LU_NONZERO)
1002       DOUBLE COMPLEX :: E2(LU_NONZERO)
1003 #endif
1004       KPP_REAL :: Z1(N),Z2(N),Z3(N),   &
1005                F1(N),F2(N),F3(N),CONT(N),   &
1006                GAMMA,ALPHA,BETA
1007       DOUBLE COMPLEX :: BC(N)
1008       INTEGER :: i,j
1009       KPP_REAL :: S2, S3
1010 !
1011       DO i=1,N
1012          S2=-F2(i)
1013          S3=-F3(i)
1014          Z1(i)=Z1(i)-F1(i)*GAMMA
1015          Z2(i)=Z2(i)+S2*ALPHA-S3*BETA
1016          Z3(i)=Z3(i)+S3*ALPHA+S2*BETA
1017       END DO
1018 #ifdef FULL_ALGEBRA
1019       CALL DGETRS ('N',N,1,E1,N,IP1,Z1,N,0)
1020 #else
1021       CALL KppSolve (E1,Z1)
1022 #endif
1023
1024       DO j=1,N
1025         BC(j) = DCMPLX(Z2(j),Z3(j))
1026       END DO
1027 #ifdef FULL_ALGEBRA
1028       CALL ZGETRS ('N',N,1,E2,N,IP2,BC,N,0)
1029 #else
1030       CALL KppSolveCmplx (E2,BC)
1031 #endif
1032       DO j=1,N
1033         Z2(j) = DBLE( BC(j) )
1034         Z3(j) = AIMAG( BC(j) )
1035       END DO
1036
1037    END SUBROUTINE RAD_Solve
1038
1039
1040 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1041    SUBROUTINE RAD_ErrorEstimate(N,FJAC,H,Y0,Y,T,&
1042                E1,Z1,Z2,Z3,IP1,SCAL,ERR,        &
1043                FirstStep,REJECT,GAMMA)
1044 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1045       IMPLICIT NONE
1046
1047       INTEGER :: N
1048 #ifdef FULL_ALGEBRA
1049       KPP_REAL :: FJAC(NVAR,NVAR),  E1(NVAR,NVAR)
1050 #else
1051       KPP_REAL :: FJAC(LU_NONZERO), E1(LU_NONZERO)
1052 #endif
1053       KPP_REAL :: SCAL(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N), &
1054                Y0(N),Y(N),TMP(N),T,H,GAMMA
1055       INTEGER :: IP1(N), i
1056       LOGICAL FirstStep,REJECT
1057       KPP_REAL :: HEE1,HEE2,HEE3,ERR
1058
1059       HEE1 = rkE(1)/H
1060       HEE2 = rkE(2)/H
1061       HEE3 = rkE(3)/H
1062
1063       DO  i=1,N
1064          F2(i)=HEE1*Z1(i)+HEE2*Z2(i)+HEE3*Z3(i)
1065          TMP(i)=F2(i)+Y0(i)
1066       END DO
1067
1068 #ifdef FULL_ALGEBRA
1069       CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,0)
1070 #else
1071       CALL KppSolve (E1, TMP)
1072 #endif
1073
1074       ERR=0.D0
1075       DO  i=1,N
1076          ERR=ERR+(TMP(i)/SCAL(i))**2
1077       END DO
1078       ERR=MAX(SQRT(ERR/N),1.D-10)
1079 !
1080       IF (ERR < 1.D0) RETURN
1081 firej:IF (FirstStep.OR.REJECT) THEN
1082           DO i=1,N
1083              TMP(i)=Y(i)+TMP(i)
1084           END DO
1085           CALL FUN_CHEM(T,TMP,F1)
1086           DO i=1,N
1087              TMP(i)=F1(i)+F2(i)
1088           END DO
1089
1090 #ifdef FULL_ALGEBRA
1091           CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,0)
1092 #else
1093           CALL KppSolve (E1, TMP)
1094 #endif
1095           ERR=0.D0
1096           DO i=1,N
1097              ERR=ERR+(TMP(i)/SCAL(i))**2
1098           END DO
1099           ERR=MAX(SQRT(ERR/N),1.0d-10)
1100        END IF firej
1101
1102    END SUBROUTINE RAD_ErrorEstimate
1103
1104 !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1105 !   KPP_REAL FUNCTION CONTR5(I,N,T,CONT)
1106 !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1107 !!     THIS FUNCTION CAN BE USED FOR CONTINUOUS OUTPUT. IT PROVIDES AN
1108 !!     APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT T.
1109 !!     IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR
1110 !!     THE STEP SUCCESSFULLY COMPUTED STEP (BY RADAU5).
1111 !!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1112 !      IMPLICIT NONE
1113 !      INTEGER :: I, N
1114 !      KPP_REAL :: T, CONT(N,4)
1115 !      KPP_REAL :: S
1116 !      KPP_REAL, PARAMETER :: ONE = 1.0d0
1117 !      S=(T-Tsol)/Hsol
1118 !      CONTR5=CONT(i,1)+S* &
1119 !        (CONT(i,2)+(S-rkC(2)+ONE)*(CONT(i,3)+(S-rkC(1)+ONE)*CONT(i,4)))
1120 !   END FUNCTION CONTR5
1121
1122 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1123   END SUBROUTINE RADAU5 ! AND ALL ITS INTERNAL PROCEDURES
1124 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1125
1126 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1127   SUBROUTINE FUN_CHEM(T, V, FCT)
1128 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1129
1130     USE KPP_ROOT_Parameters
1131     USE KPP_ROOT_Global
1132     USE KPP_ROOT_Function, ONLY: Fun
1133     USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1134
1135     IMPLICIT NONE
1136
1137     KPP_REAL :: V(NVAR), FCT(NVAR)
1138     KPP_REAL :: T, Told
1139
1140     !Told = TIME
1141     !TIME = T
1142     !CALL Update_SUN()
1143     !CALL Update_RCONST()
1144     !CALL Update_PHOTO()
1145     !TIME = Told
1146
1147     CALL Fun(V, FIX, RCONST, FCT)
1148
1149     Nfun=Nfun+1
1150
1151   END SUBROUTINE FUN_CHEM
1152
1153 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1154   SUBROUTINE JAC_CHEM (T, V, JF)
1155 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1156
1157     USE KPP_ROOT_Parameters
1158     USE KPP_ROOT_Global
1159     USE KPP_ROOT_JacobianSP
1160     USE KPP_ROOT_Jacobian, ONLY: Jac_SP
1161     USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1162
1163     IMPLICIT NONE
1164
1165     KPP_REAL :: V(NVAR), T, Told
1166 #ifdef FULL_ALGEBRA
1167     KPP_REAL :: JV(LU_NONZERO), JF(NVAR,NVAR)
1168     INTEGER :: i, j
1169 #else
1170     KPP_REAL :: JF(LU_NONZERO)
1171 #endif
1172
1173     !Told = TIME
1174     !TIME = T
1175     !CALL Update_SUN()
1176     !CALL Update_RCONST()
1177     !CALL Update_PHOTO()
1178     !TIME = Told
1179
1180 #ifdef FULL_ALGEBRA
1181     CALL Jac_SP(V, FIX, RCONST, JV)
1182     DO j=1,NVAR
1183       DO i=1,NVAR
1184          JF(i,j) = 0.0d0
1185       END DO
1186     END DO
1187     DO i=1,LU_NONZERO
1188        JF(LU_IROW(i),LU_ICOL(i)) = JV(i)
1189     END DO
1190 #else
1191     CALL Jac_SP(V, FIX, RCONST, JF)
1192 #endif
1193
1194     Njac=Njac+1
1195
1196   END SUBROUTINE JAC_CHEM
1197
1198 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1199
1200 END MODULE KPP_ROOT_Integrator