[2696] | 1 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! |
---|
| 2 | ! LSODE - Stiff method based on backward differentiation formulas (BDF) ! |
---|
| 3 | ! By default the code employs the KPP sparse linear algebra routines ! |
---|
| 4 | ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! |
---|
| 5 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! |
---|
| 6 | ! A. Sandu - version of July 2005 |
---|
| 7 | |
---|
| 8 | MODULE KPP_ROOT_Integrator |
---|
| 9 | |
---|
| 10 | USE KPP_ROOT_Precision |
---|
| 11 | USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME, ATOL, RTOL |
---|
| 12 | USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO |
---|
| 13 | USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG |
---|
| 14 | USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve, & |
---|
| 15 | Set2zero, WLAMCH |
---|
| 16 | |
---|
| 17 | IMPLICIT NONE |
---|
| 18 | PUBLIC |
---|
| 19 | SAVE |
---|
| 20 | |
---|
| 21 | !~~~> Statistics on the work performed by the LSODE method |
---|
| 22 | INTEGER :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng |
---|
| 23 | INTEGER, PARAMETER :: ifun=1, ijac=2, istp=3, iacc=4, & |
---|
| 24 | irej=5, idec=6, isol=7, isng=8, itexit=1, ihexit=2 |
---|
| 25 | ! SDIRK method coefficients |
---|
| 26 | KPP_REAL :: rkAlpha(5,4), rkBeta(5,4), rkD(4,5), & |
---|
| 27 | rkGamma, rkA(5,5), rkB(5), rkC(5) |
---|
| 28 | |
---|
| 29 | ! mz_rs_20050717: TODO: use strings of IERR_NAMES for error messages |
---|
| 30 | ! description of the error numbers IERR |
---|
| 31 | CHARACTER(LEN=50), PARAMETER, DIMENSION(-8:1) :: IERR_NAMES = (/ & |
---|
| 32 | 'Matrix is repeatedly singular ', & ! -8 |
---|
| 33 | 'Step size too small ', & ! -7 |
---|
| 34 | 'No of steps exceeds maximum bound ', & ! -6 |
---|
| 35 | 'Improper tolerance values ', & ! -5 |
---|
| 36 | 'FacMin/FacMax/FacRej must be positive ', & ! -4 |
---|
| 37 | 'Hmin/Hmax/Hstart must be positive ', & ! -3 |
---|
| 38 | 'Improper value for maximal no of Newton iterations', & ! -2 |
---|
| 39 | 'Improper value for maximal no of steps ', & ! -1 |
---|
| 40 | ' ', & ! 0 (not used) |
---|
| 41 | 'Success ' /) ! 1 |
---|
| 42 | |
---|
| 43 | CONTAINS |
---|
| 44 | |
---|
| 45 | SUBROUTINE INTEGRATE( TIN, TOUT, & |
---|
| 46 | ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U ) |
---|
| 47 | |
---|
| 48 | USE KPP_ROOT_Parameters |
---|
| 49 | USE KPP_ROOT_Global |
---|
| 50 | IMPLICIT NONE |
---|
| 51 | |
---|
| 52 | KPP_REAL, INTENT(IN) :: TIN ! Start Time |
---|
| 53 | KPP_REAL, INTENT(IN) :: TOUT ! End Time |
---|
| 54 | ! Optional input parameters and statistics |
---|
| 55 | INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20) |
---|
| 56 | KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20) |
---|
| 57 | INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20) |
---|
| 58 | KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20) |
---|
| 59 | INTEGER, INTENT(OUT), OPTIONAL :: IERR_U |
---|
| 60 | |
---|
| 61 | KPP_REAL :: RCNTRL(20), RSTATUS(20) |
---|
| 62 | INTEGER :: ICNTRL(20), ISTATUS(20), IERR |
---|
| 63 | !!$ INTEGER, SAVE :: Ntotal = 0 |
---|
| 64 | |
---|
| 65 | ICNTRL(:) = 0 |
---|
| 66 | RCNTRL(:) = 0.0_dp |
---|
| 67 | ISTATUS(:) = 0 |
---|
| 68 | RSTATUS(:) = 0.0_dp |
---|
| 69 | |
---|
| 70 | ICNTRL(5) = 2 ! maximal order |
---|
| 71 | |
---|
| 72 | ! If optional parameters are given, and if they are >0, |
---|
| 73 | ! then they overwrite default settings. |
---|
| 74 | IF (PRESENT(ICNTRL_U)) THEN |
---|
| 75 | WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:) |
---|
| 76 | END IF |
---|
| 77 | IF (PRESENT(RCNTRL_U)) THEN |
---|
| 78 | WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:) |
---|
| 79 | END IF |
---|
| 80 | |
---|
| 81 | CALL KppLsode( TIN,TOUT,VAR,RTOL,ATOL, & |
---|
| 82 | RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR ) |
---|
| 83 | |
---|
| 84 | ! mz_rs_20050716: IERR and ISTATUS are returned to the user who then |
---|
| 85 | ! decides what to do about it, i.e. either stop the run or ignore it. |
---|
| 86 | !!$ IF (IERR < 0) THEN |
---|
| 87 | !!$ PRINT*,'LSODE: Unsuccessful exit at T=',TIN,' (IERR=',IERR,')' |
---|
| 88 | !!$ STOP |
---|
| 89 | !!$ ENDIF |
---|
| 90 | !!$ Ntotal = Ntotal + ISTATUS(3) |
---|
| 91 | !!$ PRINT*,'Nsteps = ', ISTATUS(3),' (',Ntotal,')' |
---|
| 92 | |
---|
| 93 | STEPMIN = RSTATUS(ihexit) ! Save last step |
---|
| 94 | |
---|
| 95 | ! if optional parameters are given for output they to return information |
---|
| 96 | IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(1:20) |
---|
| 97 | IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(1:20) |
---|
| 98 | IF (PRESENT(IERR_U)) THEN |
---|
| 99 | IF (IERR==2) THEN ! DLSODE returns "2" after successful completion |
---|
| 100 | IERR_U = 1 ! IERR_U will return "1" for successful completion |
---|
| 101 | ELSE |
---|
| 102 | IERR_U = IERR |
---|
| 103 | ENDIF |
---|
| 104 | ENDIF |
---|
| 105 | |
---|
| 106 | END SUBROUTINE INTEGRATE |
---|
| 107 | |
---|
| 108 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 109 | SUBROUTINE KppLsode( TIN,TOUT,Y,RelTol,AbsTol, & |
---|
| 110 | RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR ) |
---|
| 111 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 112 | ! |
---|
| 113 | !~~~> INPUT PARAMETERS: |
---|
| 114 | ! |
---|
| 115 | ! Note: For input parameters equal to zero the default values of the |
---|
| 116 | ! corresponding variables are used. |
---|
| 117 | ! |
---|
| 118 | ! Note: For input parameters equal to zero the default values of the |
---|
| 119 | ! corresponding variables are used. |
---|
| 120 | !~~~> |
---|
| 121 | ! ICNTRL(1) = not used |
---|
| 122 | ! |
---|
| 123 | ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors |
---|
| 124 | ! = 1: AbsTol, RelTol are scalars |
---|
| 125 | ! |
---|
| 126 | ! ICNTRL(3) = not used |
---|
| 127 | ! |
---|
| 128 | ! ICNTRL(4) -> maximum number of integration steps |
---|
| 129 | ! For ICNTRL(4)=0 the default value of 100000 is used |
---|
| 130 | ! |
---|
| 131 | ! ICNTRL(5) -> maximum order of the integration formula allowed |
---|
| 132 | ! |
---|
| 133 | !~~~> Real parameters |
---|
| 134 | ! |
---|
| 135 | ! RCNTRL(1) -> Hmin, lower bound for the integration step size |
---|
| 136 | ! It is strongly recommended to keep Hmin = ZERO |
---|
| 137 | ! RCNTRL(2) -> Hmax, upper bound for the integration step size |
---|
| 138 | ! RCNTRL(3) -> Hstart, starting value for the integration step size |
---|
| 139 | ! |
---|
| 140 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 141 | ! |
---|
| 142 | !~~~> OUTPUT PARAMETERS: |
---|
| 143 | ! |
---|
| 144 | ! Note: each call to Rosenbrock adds the current no. of fcn calls |
---|
| 145 | ! to previous value of ISTATUS(1), and similar for the other params. |
---|
| 146 | ! Set ISTATUS(1:10) = 0 before call to avoid this accumulation. |
---|
| 147 | ! |
---|
| 148 | ! ISTATUS(1) = No. of function calls |
---|
| 149 | ! ISTATUS(2) = No. of jacobian calls |
---|
| 150 | ! ISTATUS(3) = No. of steps |
---|
| 151 | ! |
---|
| 152 | ! RSTATUS(1) -> Texit, the time corresponding to the |
---|
| 153 | ! computed Y upon return |
---|
| 154 | ! RSTATUS(2) -> Hexit, last predicted step before exit |
---|
| 155 | ! For multiple restarts, use Hexit as Hstart in the following run |
---|
| 156 | ! |
---|
| 157 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 158 | |
---|
| 159 | IMPLICIT NONE |
---|
| 160 | KPP_REAL :: Y(NVAR), AbsTol(NVAR), RelTol(NVAR), TIN, TOUT |
---|
| 161 | KPP_REAL :: RCNTRL(20), RSTATUS(20) |
---|
| 162 | INTEGER :: ICNTRL(20), ISTATUS(20) |
---|
| 163 | INTEGER, PARAMETER :: LRW = 25 + 9*NVAR+2*NVAR*NVAR, & |
---|
| 164 | LIW = 32 + NVAR |
---|
| 165 | KPP_REAL :: RWORK(LRW), RPAR(1) |
---|
| 166 | INTEGER :: IWORK(LIW), IPAR(1), ITOL, ITASK, & |
---|
| 167 | IERR, IOPT, MF |
---|
| 168 | |
---|
| 169 | !~~~> NORMAL COMPUTATION |
---|
| 170 | ITASK = 1 |
---|
| 171 | IERR = 1 |
---|
| 172 | IOPT = 1 ! 0=no/1=use optional input |
---|
| 173 | |
---|
| 174 | RWORK(1:30) = 0.0d0 |
---|
| 175 | IWORK(1:30) = 0 |
---|
| 176 | |
---|
| 177 | IF (ICNTRL(2)==0) THEN |
---|
| 178 | ITOL = 4 ! Abs/RelTol are both vectors |
---|
| 179 | ELSE |
---|
| 180 | ITOL = 1 ! Abs/RelTol are both scalars |
---|
| 181 | END IF |
---|
| 182 | IWORK(6) = ICNTRL(4) ! max number of internal steps |
---|
| 183 | IWORK(5) = ICNTRL(5) ! maximal order |
---|
| 184 | |
---|
| 185 | MF = 21 !~~~> stiff case, analytic full Jacobian |
---|
| 186 | |
---|
| 187 | RWORK(5) = RCNTRL(3) ! Hstart |
---|
| 188 | RWORK(6) = RCNTRL(2) ! Hmax |
---|
| 189 | RWORK(7) = RCNTRL(1) ! Hmin |
---|
| 190 | |
---|
| 191 | CALL DLSODE ( FUN_CHEM, NVAR, Y, TIN, TOUT, ITOL, RelTol, AbsTol, ITASK,& |
---|
| 192 | IERR, IOPT, RWORK, LRW, IWORK, LIW, JAC_CHEM, MF) |
---|
| 193 | |
---|
| 194 | ISTATUS(1) = IWORK(12) ! Number of function evaluations |
---|
| 195 | ISTATUS(2) = IWORK(13) ! Number of Jacobian evaluations |
---|
| 196 | ISTATUS(3) = IWORK(11) ! Number of steps |
---|
| 197 | |
---|
| 198 | RSTATUS(1) = TOUT ! mz_rs_20050717 |
---|
| 199 | RSTATUS(2) = RWORK(11) ! mz_rs_20050717 |
---|
| 200 | |
---|
| 201 | END SUBROUTINE KppLsode |
---|
| 202 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 203 | |
---|
| 204 | |
---|
| 205 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 206 | !DECK DLSODE |
---|
| 207 | SUBROUTINE DLSODE (F, NEQ, Y, T, TOUT, ITOL, RelTol, AbsTol, ITASK, & |
---|
| 208 | ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF) |
---|
| 209 | EXTERNAL F, JAC |
---|
| 210 | INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, LIW, IWORK(LIW), MF |
---|
| 211 | KPP_REAL Y(*), T, TOUT, RelTol(*), AbsTol(*), RWORK(LRW) |
---|
| 212 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 213 | !***BEGIN PROLOGUE DLSODE |
---|
| 214 | !***PURPOSE Livermore Solver for Ordinary Differential Equations. |
---|
| 215 | ! DLSODE solves the initial-value problem for stiff or |
---|
| 216 | ! nonstiff systems of first-order ODE's, |
---|
| 217 | ! dy/dt = f(t,y), or, in component form, |
---|
| 218 | ! dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(N)), i=1,...,N. |
---|
| 219 | !***CATEGORY I1A |
---|
| 220 | !***TYPE KPP_REAL (SLSODE-S, DLSODE-D) |
---|
| 221 | !***KEYWORDS ORDINARY DIFFERENTIAL EQUATIONS, INITIAL VALUE PROBLEM, |
---|
| 222 | ! STIFF, NONSTIFF |
---|
| 223 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 224 | ! Center for Applied Scientific Computing, L-561 |
---|
| 225 | ! Lawrence Livermore National Laboratory |
---|
| 226 | ! Livermore, CA 94551. |
---|
| 227 | !***DESCRIPTION |
---|
| 228 | ! |
---|
| 229 | ! NOTE: The "Usage" and "Arguments" sections treat only a subset of |
---|
| 230 | ! available options, in condensed fashion. The options |
---|
| 231 | ! covered and the information supplied will support most |
---|
| 232 | ! standard uses of DLSODE. |
---|
| 233 | ! |
---|
| 234 | ! For more sophisticated uses, full details on all options are |
---|
| 235 | ! given in the concluding section, headed "Long Description." |
---|
| 236 | ! A synopsis of the DLSODE Long Description is provided at the |
---|
| 237 | ! beginning of that section; general topics covered are: |
---|
| 238 | ! - Elements of the call sequence; optional input and output |
---|
| 239 | ! - Optional supplemental routines in the DLSODE package |
---|
| 240 | ! - internal COMMON block |
---|
| 241 | ! |
---|
| 242 | ! *Usage: |
---|
| 243 | ! Communication between the user and the DLSODE package, for normal |
---|
| 244 | ! situations, is summarized here. This summary describes a subset |
---|
| 245 | ! of the available options. See "Long Description" for complete |
---|
| 246 | ! details, including optional communication, nonstandard options, |
---|
| 247 | ! and instructions for special situations. |
---|
| 248 | ! |
---|
| 249 | ! A sample program is given in the "Examples" section. |
---|
| 250 | ! |
---|
| 251 | ! Refer to the argument descriptions for the definitions of the |
---|
| 252 | ! quantities that appear in the following sample declarations. |
---|
| 253 | ! |
---|
| 254 | ! For MF = 10, |
---|
| 255 | ! PARAMETER (LRW = 20 + 16*NEQ, LIW = 20) |
---|
| 256 | ! For MF = 21 or 22, |
---|
| 257 | ! PARAMETER (LRW = 22 + 9*NEQ + NEQ**2, LIW = 20 + NEQ) |
---|
| 258 | ! For MF = 24 or 25, |
---|
| 259 | ! PARAMETER (LRW = 22 + 10*NEQ + (2*ML+MU)*NEQ, |
---|
| 260 | ! * LIW = 20 + NEQ) |
---|
| 261 | ! |
---|
| 262 | ! EXTERNAL F, JAC |
---|
| 263 | ! INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK(LIW), |
---|
| 264 | ! * LIW, MF |
---|
| 265 | ! KPP_REAL Y(NEQ), T, TOUT, RelTol, AbsTol(ntol), RWORK(LRW) |
---|
| 266 | ! |
---|
| 267 | ! CALL DLSODE (F, NEQ, Y, T, TOUT, ITOL, RelTol, AbsTol, ITASK, |
---|
| 268 | ! * ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF) |
---|
| 269 | ! |
---|
| 270 | ! *Arguments: |
---|
| 271 | ! F :EXT Name of subroutine for right-hand-side vector f. |
---|
| 272 | ! This name must be declared EXTERNAL in calling |
---|
| 273 | ! program. The form of F must be: |
---|
| 274 | ! |
---|
| 275 | ! SUBROUTINE F (NEQ, T, Y, YDOT) |
---|
| 276 | ! INTEGER NEQ |
---|
| 277 | ! KPP_REAL T, Y(*), YDOT(*) |
---|
| 278 | ! |
---|
| 279 | ! The inputs are NEQ, T, Y. F is to set |
---|
| 280 | ! |
---|
| 281 | ! YDOT(i) = f(i,T,Y(1),Y(2),...,Y(NEQ)), |
---|
| 282 | ! i = 1, ..., NEQ . |
---|
| 283 | ! |
---|
| 284 | ! NEQ :IN Number of first-order ODE's. |
---|
| 285 | ! |
---|
| 286 | ! Y :INOUT Array of values of the y(t) vector, of length NEQ. |
---|
| 287 | ! Input: For the first call, Y should contain the |
---|
| 288 | ! values of y(t) at t = T. (Y is an input |
---|
| 289 | ! variable only if ISTATE = 1.) |
---|
| 290 | ! Output: On return, Y will contain the values at the |
---|
| 291 | ! new t-value. |
---|
| 292 | ! |
---|
| 293 | ! T :INOUT Value of the independent variable. On return it |
---|
| 294 | ! will be the current value of t (normally TOUT). |
---|
| 295 | ! |
---|
| 296 | ! TOUT :IN Next point where output is desired (.NE. T). |
---|
| 297 | ! |
---|
| 298 | ! ITOL :IN 1 or 2 according as AbsTol (below) is a scalar or |
---|
| 299 | ! an array. |
---|
| 300 | ! |
---|
| 301 | ! RelTol :IN Relative tolerance parameter (scalar). |
---|
| 302 | ! |
---|
| 303 | ! AbsTol :IN Absolute tolerance parameter (scalar or array). |
---|
| 304 | ! If ITOL = 1, AbsTol need not be dimensioned. |
---|
| 305 | ! If ITOL = 2, AbsTol must be dimensioned at least NEQ. |
---|
| 306 | ! |
---|
| 307 | ! The estimated local error in Y(i) will be controlled |
---|
| 308 | ! so as to be roughly less (in magnitude) than |
---|
| 309 | ! |
---|
| 310 | ! EWT(i) = RelTol*ABS(Y(i)) + AbsTol if ITOL = 1, or |
---|
| 311 | ! EWT(i) = RelTol*ABS(Y(i)) + AbsTol(i) if ITOL = 2. |
---|
| 312 | ! |
---|
| 313 | ! Thus the local error test passes if, in each |
---|
| 314 | ! component, either the absolute error is less than |
---|
| 315 | ! AbsTol (or AbsTol(i)), or the relative error is less |
---|
| 316 | ! than RelTol. |
---|
| 317 | ! |
---|
| 318 | ! Use RelTol = 0.0 for pure absolute error control, and |
---|
| 319 | ! use AbsTol = 0.0 (or AbsTol(i) = 0.0) for pure relative |
---|
| 320 | ! error control. Caution: Actual (global) errors may |
---|
| 321 | ! exceed these local tolerances, so choose them |
---|
| 322 | ! conservatively. |
---|
| 323 | ! |
---|
| 324 | ! ITASK :IN Flag indicating the task DLSODE is to perform. |
---|
| 325 | ! Use ITASK = 1 for normal computation of output |
---|
| 326 | ! values of y at t = TOUT. |
---|
| 327 | ! |
---|
| 328 | ! ISTATE:INOUT Index used for input and output to specify the state |
---|
| 329 | ! of the calculation. |
---|
| 330 | ! Input: |
---|
| 331 | ! 1 This is the first call for a problem. |
---|
| 332 | ! 2 This is a subsequent call. |
---|
| 333 | ! Output: |
---|
| 334 | ! 1 Nothing was done, because TOUT was equal to T. |
---|
| 335 | ! 2 DLSODE was successful (otherwise, negative). |
---|
| 336 | ! Note that ISTATE need not be modified after a |
---|
| 337 | ! successful return. |
---|
| 338 | ! -1 Excess work done on this call (perhaps wrong |
---|
| 339 | ! MF). |
---|
| 340 | ! -2 Excess accuracy requested (tolerances too |
---|
| 341 | ! small). |
---|
| 342 | ! -3 Illegal input detected (see printed message). |
---|
| 343 | ! -4 Repeated error test failures (check all |
---|
| 344 | ! inputs). |
---|
| 345 | ! -5 Repeated convergence failures (perhaps bad |
---|
| 346 | ! Jacobian supplied or wrong choice of MF or |
---|
| 347 | ! tolerances). |
---|
| 348 | ! -6 Error weight became zero during problem |
---|
| 349 | ! (solution component i vanished, and AbsTol or |
---|
| 350 | ! AbsTol(i) = 0.). |
---|
| 351 | ! |
---|
| 352 | ! IOPT :IN Flag indicating whether optional inputs are used: |
---|
| 353 | ! 0 No. |
---|
| 354 | ! 1 Yes. (See "Optional inputs" under "Long |
---|
| 355 | ! Description," Part 1.) |
---|
| 356 | ! |
---|
| 357 | ! RWORK :WORK Real work array of length at least: |
---|
| 358 | ! 20 + 16*NEQ for MF = 10, |
---|
| 359 | ! 22 + 9*NEQ + NEQ**2 for MF = 21 or 22, |
---|
| 360 | ! 22 + 10*NEQ + (2*ML + MU)*NEQ for MF = 24 or 25. |
---|
| 361 | ! |
---|
| 362 | ! LRW :IN Declared length of RWORK (in user's DIMENSION |
---|
| 363 | ! statement). |
---|
| 364 | ! |
---|
| 365 | ! IWORK :WORK Integer work array of length at least: |
---|
| 366 | ! 20 for MF = 10, |
---|
| 367 | ! 20 + NEQ for MF = 21, 22, 24, or 25. |
---|
| 368 | ! |
---|
| 369 | ! If MF = 24 or 25, input in IWORK(1),IWORK(2) the |
---|
| 370 | ! lower and upper Jacobian half-bandwidths ML,MU. |
---|
| 371 | ! |
---|
| 372 | ! On return, IWORK contains information that may be |
---|
| 373 | ! of interest to the user: |
---|
| 374 | ! |
---|
| 375 | ! Name Location Meaning |
---|
| 376 | ! ----- --------- ----------------------------------------- |
---|
| 377 | ! NST IWORK(11) Number of steps taken for the problem so |
---|
| 378 | ! far. |
---|
| 379 | ! NFE IWORK(12) Number of f evaluations for the problem |
---|
| 380 | ! so far. |
---|
| 381 | ! NJE IWORK(13) Number of Jacobian evaluations (and of |
---|
| 382 | ! matrix LU decompositions) for the problem |
---|
| 383 | ! so far. |
---|
| 384 | ! NQU IWORK(14) Method order last used (successfully). |
---|
| 385 | ! LENRW IWORK(17) Length of RWORK actually required. This |
---|
| 386 | ! is defined on normal returns and on an |
---|
| 387 | ! illegal input return for insufficient |
---|
| 388 | ! storage. |
---|
| 389 | ! LENIW IWORK(18) Length of IWORK actually required. This |
---|
| 390 | ! is defined on normal returns and on an |
---|
| 391 | ! illegal input return for insufficient |
---|
| 392 | ! storage. |
---|
| 393 | ! |
---|
| 394 | ! LIW :IN Declared length of IWORK (in user's DIMENSION |
---|
| 395 | ! statement). |
---|
| 396 | ! |
---|
| 397 | ! JAC :EXT Name of subroutine for Jacobian matrix (MF = |
---|
| 398 | ! 21 or 24). If used, this name must be declared |
---|
| 399 | ! EXTERNAL in calling program. If not used, pass a |
---|
| 400 | ! dummy name. The form of JAC must be: |
---|
| 401 | ! |
---|
| 402 | ! SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD) |
---|
| 403 | ! INTEGER NEQ, ML, MU, NROWPD |
---|
| 404 | ! KPP_REAL T, Y(*), PD(NROWPD,*) |
---|
| 405 | ! |
---|
| 406 | ! See item c, under "Description" below for more |
---|
| 407 | ! information about JAC. |
---|
| 408 | ! |
---|
| 409 | ! MF :IN Method flag. Standard values are: |
---|
| 410 | ! 10 Nonstiff (Adams) method, no Jacobian used. |
---|
| 411 | ! 21 Stiff (BDF) method, user-supplied full Jacobian. |
---|
| 412 | ! 22 Stiff method, internally generated full |
---|
| 413 | ! Jacobian. |
---|
| 414 | ! 24 Stiff method, user-supplied banded Jacobian. |
---|
| 415 | ! 25 Stiff method, internally generated banded |
---|
| 416 | ! Jacobian. |
---|
| 417 | ! |
---|
| 418 | ! *Description: |
---|
| 419 | ! DLSODE solves the initial value problem for stiff or nonstiff |
---|
| 420 | ! systems of first-order ODE's, |
---|
| 421 | ! |
---|
| 422 | ! dy/dt = f(t,y) , |
---|
| 423 | ! |
---|
| 424 | ! or, in component form, |
---|
| 425 | ! |
---|
| 426 | ! dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) |
---|
| 427 | ! (i = 1, ..., NEQ) . |
---|
| 428 | ! |
---|
| 429 | ! DLSODE is a package based on the GEAR and GEARB packages, and on |
---|
| 430 | ! the October 23, 1978, version of the tentative ODEPACK user |
---|
| 431 | ! interface standard, with minor modifications. |
---|
| 432 | ! |
---|
| 433 | ! The steps in solving such a problem are as follows. |
---|
| 434 | ! |
---|
| 435 | ! a. First write a subroutine of the form |
---|
| 436 | ! |
---|
| 437 | ! SUBROUTINE F (NEQ, T, Y, YDOT) |
---|
| 438 | ! INTEGER NEQ |
---|
| 439 | ! KPP_REAL T, Y(*), YDOT(*) |
---|
| 440 | ! |
---|
| 441 | ! which supplies the vector function f by loading YDOT(i) with |
---|
| 442 | ! f(i). |
---|
| 443 | ! |
---|
| 444 | ! b. Next determine (or guess) whether or not the problem is stiff. |
---|
| 445 | ! Stiffness occurs when the Jacobian matrix df/dy has an |
---|
| 446 | ! eigenvalue whose real part is negative and large in magnitude |
---|
| 447 | ! compared to the reciprocal of the t span of interest. If the |
---|
| 448 | ! problem is nonstiff, use method flag MF = 10. If it is stiff, |
---|
| 449 | ! there are four standard choices for MF, and DLSODE requires the |
---|
| 450 | ! Jacobian matrix in some form. This matrix is regarded either |
---|
| 451 | ! as full (MF = 21 or 22), or banded (MF = 24 or 25). In the |
---|
| 452 | ! banded case, DLSODE requires two half-bandwidth parameters ML |
---|
| 453 | ! and MU. These are, respectively, the widths of the lower and |
---|
| 454 | ! upper parts of the band, excluding the main diagonal. Thus the |
---|
| 455 | ! band consists of the locations (i,j) with |
---|
| 456 | ! |
---|
| 457 | ! i - ML <= j <= i + MU , |
---|
| 458 | ! |
---|
| 459 | ! and the full bandwidth is ML + MU + 1 . |
---|
| 460 | ! |
---|
| 461 | ! c. If the problem is stiff, you are encouraged to supply the |
---|
| 462 | ! Jacobian directly (MF = 21 or 24), but if this is not feasible, |
---|
| 463 | ! DLSODE will compute it internally by difference quotients (MF = |
---|
| 464 | ! 22 or 25). If you are supplying the Jacobian, write a |
---|
| 465 | ! subroutine of the form |
---|
| 466 | ! |
---|
| 467 | ! SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD) |
---|
| 468 | ! INTEGER NEQ, ML, MU, NRWOPD |
---|
| 469 | ! KPP_REAL T, Y(*), PD(NROWPD,*) |
---|
| 470 | ! |
---|
| 471 | ! which provides df/dy by loading PD as follows: |
---|
| 472 | ! - For a full Jacobian (MF = 21), load PD(i,j) with df(i)/dy(j), |
---|
| 473 | ! the partial derivative of f(i) with respect to y(j). (Ignore |
---|
| 474 | ! the ML and MU arguments in this case.) |
---|
| 475 | ! - For a banded Jacobian (MF = 24), load PD(i-j+MU+1,j) with |
---|
| 476 | ! df(i)/dy(j); i.e., load the diagonal lines of df/dy into the |
---|
| 477 | ! rows of PD from the top down. |
---|
| 478 | ! - In either case, only nonzero elements need be loaded. |
---|
| 479 | ! |
---|
| 480 | ! d. Write a main program that calls subroutine DLSODE once for each |
---|
| 481 | ! point at which answers are desired. This should also provide |
---|
| 482 | ! for possible use of logical unit 6 for output of error messages |
---|
| 483 | ! by DLSODE. |
---|
| 484 | ! |
---|
| 485 | ! Before the first call to DLSODE, set ISTATE = 1, set Y and T to |
---|
| 486 | ! the initial values, and set TOUT to the first output point. To |
---|
| 487 | ! continue the integration after a successful return, simply |
---|
| 488 | ! reset TOUT and call DLSODE again. No other parameters need be |
---|
| 489 | ! reset. |
---|
| 490 | ! |
---|
| 491 | ! *Examples: |
---|
| 492 | ! The following is a simple example problem, with the coding needed |
---|
| 493 | ! for its solution by DLSODE. The problem is from chemical kinetics, |
---|
| 494 | ! and consists of the following three rate equations: |
---|
| 495 | ! |
---|
| 496 | ! dy1/dt = -.04*y1 + 1.E4*y2*y3 |
---|
| 497 | ! dy2/dt = .04*y1 - 1.E4*y2*y3 - 3.E7*y2**2 |
---|
| 498 | ! dy3/dt = 3.E7*y2**2 |
---|
| 499 | ! |
---|
| 500 | ! on the interval from t = 0.0 to t = 4.E10, with initial conditions |
---|
| 501 | ! y1 = 1.0, y2 = y3 = 0. The problem is stiff. |
---|
| 502 | ! |
---|
| 503 | ! The following coding solves this problem with DLSODE, using |
---|
| 504 | ! MF = 21 and printing results at t = .4, 4., ..., 4.E10. It uses |
---|
| 505 | ! ITOL = 2 and AbsTol much smaller for y2 than for y1 or y3 because y2 |
---|
| 506 | ! has much smaller values. At the end of the run, statistical |
---|
| 507 | ! quantities of interest are printed. |
---|
| 508 | ! |
---|
| 509 | ! EXTERNAL FEX, JEX |
---|
| 510 | ! INTEGER IOPT, IOUT, ISTATE, ITASK, ITOL, IWORK(23), LIW, LRW, |
---|
| 511 | ! * MF, NEQ |
---|
| 512 | ! KPP_REAL AbsTol(3), RelTol, RWORK(58), T, TOUT, Y(3) |
---|
| 513 | ! NEQ = 3 |
---|
| 514 | ! Y(1) = 1.D0 |
---|
| 515 | ! Y(2) = 0.D0 |
---|
| 516 | ! Y(3) = 0.D0 |
---|
| 517 | ! T = 0.D0 |
---|
| 518 | ! TOUT = .4D0 |
---|
| 519 | ! ITOL = 2 |
---|
| 520 | ! RelTol = 1.D-4 |
---|
| 521 | ! AbsTol(1) = 1.D-6 |
---|
| 522 | ! AbsTol(2) = 1.D-10 |
---|
| 523 | ! AbsTol(3) = 1.D-6 |
---|
| 524 | ! ITASK = 1 |
---|
| 525 | ! ISTATE = 1 |
---|
| 526 | ! IOPT = 0 |
---|
| 527 | ! LRW = 58 |
---|
| 528 | ! LIW = 23 |
---|
| 529 | ! MF = 21 |
---|
| 530 | ! DO 40 IOUT = 1,12 |
---|
| 531 | ! CALL DLSODE (FEX, NEQ, Y, T, TOUT, ITOL, RelTol, AbsTol, ITASK, |
---|
| 532 | ! * ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEX, MF) |
---|
| 533 | ! WRITE(6,20) T, Y(1), Y(2), Y(3) |
---|
| 534 | ! 20 FORMAT(' At t =',D12.4,' y =',3D14.6) |
---|
| 535 | ! IF (ISTATE .LT. 0) GO TO 80 |
---|
| 536 | ! 40 TOUT = TOUT*10.D0 |
---|
| 537 | ! WRITE(6,60) IWORK(11), IWORK(12), IWORK(13) |
---|
| 538 | ! 60 FORMAT(/' No. steps =',i4,', No. f-s =',i4,', No. J-s =',i4) |
---|
| 539 | ! STOP |
---|
| 540 | ! 80 WRITE(6,90) ISTATE |
---|
| 541 | ! 90 FORMAT(///' Error halt.. ISTATE =',I3) |
---|
| 542 | ! STOP |
---|
| 543 | ! END |
---|
| 544 | ! |
---|
| 545 | ! SUBROUTINE FEX (NEQ, T, Y, YDOT) |
---|
| 546 | ! INTEGER NEQ |
---|
| 547 | ! KPP_REAL T, Y(3), YDOT(3) |
---|
| 548 | ! YDOT(1) = -.04D0*Y(1) + 1.D4*Y(2)*Y(3) |
---|
| 549 | ! YDOT(3) = 3.D7*Y(2)*Y(2) |
---|
| 550 | ! YDOT(2) = -YDOT(1) - YDOT(3) |
---|
| 551 | ! RETURN |
---|
| 552 | ! END |
---|
| 553 | ! |
---|
| 554 | ! SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD) |
---|
| 555 | ! INTEGER NEQ, ML, MU, NRPD |
---|
| 556 | ! KPP_REAL T, Y(3), PD(NRPD,3) |
---|
| 557 | ! PD(1,1) = -.04D0 |
---|
| 558 | ! PD(1,2) = 1.D4*Y(3) |
---|
| 559 | ! PD(1,3) = 1.D4*Y(2) |
---|
| 560 | ! PD(2,1) = .04D0 |
---|
| 561 | ! PD(2,3) = -PD(1,3) |
---|
| 562 | ! PD(3,2) = 6.D7*Y(2) |
---|
| 563 | ! PD(2,2) = -PD(1,2) - PD(3,2) |
---|
| 564 | ! RETURN |
---|
| 565 | ! END |
---|
| 566 | ! |
---|
| 567 | ! The output from this program (on a Cray-1 in single precision) |
---|
| 568 | ! is as follows. |
---|
| 569 | ! |
---|
| 570 | ! At t = 4.0000e-01 y = 9.851726e-01 3.386406e-05 1.479357e-02 |
---|
| 571 | ! At t = 4.0000e+00 y = 9.055142e-01 2.240418e-05 9.446344e-02 |
---|
| 572 | ! At t = 4.0000e+01 y = 7.158050e-01 9.184616e-06 2.841858e-01 |
---|
| 573 | ! At t = 4.0000e+02 y = 4.504846e-01 3.222434e-06 5.495122e-01 |
---|
| 574 | ! At t = 4.0000e+03 y = 1.831701e-01 8.940379e-07 8.168290e-01 |
---|
| 575 | ! At t = 4.0000e+04 y = 3.897016e-02 1.621193e-07 9.610297e-01 |
---|
| 576 | ! At t = 4.0000e+05 y = 4.935213e-03 1.983756e-08 9.950648e-01 |
---|
| 577 | ! At t = 4.0000e+06 y = 5.159269e-04 2.064759e-09 9.994841e-01 |
---|
| 578 | ! At t = 4.0000e+07 y = 5.306413e-05 2.122677e-10 9.999469e-01 |
---|
| 579 | ! At t = 4.0000e+08 y = 5.494530e-06 2.197825e-11 9.999945e-01 |
---|
| 580 | ! At t = 4.0000e+09 y = 5.129458e-07 2.051784e-12 9.999995e-01 |
---|
| 581 | ! At t = 4.0000e+10 y = -7.170603e-08 -2.868241e-13 1.000000e+00 |
---|
| 582 | ! |
---|
| 583 | ! No. steps = 330, No. f-s = 405, No. J-s = 69 |
---|
| 584 | ! |
---|
| 585 | ! *Accuracy: |
---|
| 586 | ! The accuracy of the solution depends on the choice of tolerances |
---|
| 587 | ! RelTol and AbsTol. Actual (global) errors may exceed these local |
---|
| 588 | ! tolerances, so choose them conservatively. |
---|
| 589 | ! |
---|
| 590 | ! *Cautions: |
---|
| 591 | ! The work arrays should not be altered between calls to DLSODE for |
---|
| 592 | ! the same problem, except possibly for the conditional and optional |
---|
| 593 | ! inputs. |
---|
| 594 | ! |
---|
| 595 | ! *Portability: |
---|
| 596 | ! Since NEQ is dimensioned inside DLSODE, some compilers may object |
---|
| 597 | ! to a call to DLSODE with NEQ a scalar variable. In this event, |
---|
| 598 | ! use DIMENSION NEQ. Similar remarks apply to RelTol and AbsTol. |
---|
| 599 | ! |
---|
| 600 | ! Note to Cray users: |
---|
| 601 | ! For maximum efficiency, use the CFT77 compiler. Appropriate |
---|
| 602 | ! compiler optimization directives have been inserted for CFT77. |
---|
| 603 | ! |
---|
| 604 | ! *Reference: |
---|
| 605 | ! Alan C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE |
---|
| 606 | ! Solvers," in Scientific Computing, R. S. Stepleman, et al., Eds. |
---|
| 607 | ! (North-Holland, Amsterdam, 1983), pp. 55-64. |
---|
| 608 | ! |
---|
| 609 | ! *Long Description: |
---|
| 610 | ! The following complete description of the user interface to |
---|
| 611 | ! DLSODE consists of four parts: |
---|
| 612 | ! |
---|
| 613 | ! 1. The call sequence to subroutine DLSODE, which is a driver |
---|
| 614 | ! routine for the solver. This includes descriptions of both |
---|
| 615 | ! the call sequence arguments and user-supplied routines. |
---|
| 616 | ! Following these descriptions is a description of optional |
---|
| 617 | ! inputs available through the call sequence, and then a |
---|
| 618 | ! description of optional outputs in the work arrays. |
---|
| 619 | ! |
---|
| 620 | ! 2. Descriptions of other routines in the DLSODE package that may |
---|
| 621 | ! be (optionally) called by the user. These provide the ability |
---|
| 622 | ! to alter error message handling, save and restore the internal |
---|
| 623 | ! COMMON, and obtain specified derivatives of the solution y(t). |
---|
| 624 | ! |
---|
| 625 | ! 3. Descriptions of COMMON block to be declared in overlay or |
---|
| 626 | ! similar environments, or to be saved when doing an interrupt |
---|
| 627 | ! of the problem and continued solution later. |
---|
| 628 | ! |
---|
| 629 | ! 4. Description of two routines in the DLSODE package, either of |
---|
| 630 | ! which the user may replace with his own version, if desired. |
---|
| 631 | ! These relate to the measurement of errors. |
---|
| 632 | ! |
---|
| 633 | ! |
---|
| 634 | ! Part 1. Call Sequence |
---|
| 635 | ! ---------------------- |
---|
| 636 | ! |
---|
| 637 | ! Arguments |
---|
| 638 | ! --------- |
---|
| 639 | ! The call sequence parameters used for input only are |
---|
| 640 | ! |
---|
| 641 | ! F, NEQ, TOUT, ITOL, RelTol, AbsTol, ITASK, IOPT, LRW, LIW, JAC, MF, |
---|
| 642 | ! |
---|
| 643 | ! and those used for both input and output are |
---|
| 644 | ! |
---|
| 645 | ! Y, T, ISTATE. |
---|
| 646 | ! |
---|
| 647 | ! The work arrays RWORK and IWORK are also used for conditional and |
---|
| 648 | ! optional inputs and optional outputs. (The term output here |
---|
| 649 | ! refers to the return from subroutine DLSODE to the user's calling |
---|
| 650 | ! program.) |
---|
| 651 | ! |
---|
| 652 | ! The legality of input parameters will be thoroughly checked on the |
---|
| 653 | ! initial call for the problem, but not checked thereafter unless a |
---|
| 654 | ! change in input parameters is flagged by ISTATE = 3 on input. |
---|
| 655 | ! |
---|
| 656 | ! The descriptions of the call arguments are as follows. |
---|
| 657 | ! |
---|
| 658 | ! F The name of the user-supplied subroutine defining the ODE |
---|
| 659 | ! system. The system must be put in the first-order form |
---|
| 660 | ! dy/dt = f(t,y), where f is a vector-valued function of |
---|
| 661 | ! the scalar t and the vector y. Subroutine F is to compute |
---|
| 662 | ! the function f. It is to have the form |
---|
| 663 | ! |
---|
| 664 | ! SUBROUTINE F (NEQ, T, Y, YDOT) |
---|
| 665 | ! KPP_REAL T, Y(*), YDOT(*) |
---|
| 666 | ! |
---|
| 667 | ! where NEQ, T, and Y are input, and the array YDOT = |
---|
| 668 | ! f(T,Y) is output. Y and YDOT are arrays of length NEQ. |
---|
| 669 | ! Subroutine F should not alter Y(1),...,Y(NEQ). F must be |
---|
| 670 | ! declared EXTERNAL in the calling program. |
---|
| 671 | ! |
---|
| 672 | ! Subroutine F may access user-defined quantities in |
---|
| 673 | ! NEQ(2),... and/or in Y(NEQ+1),..., if NEQ is an array |
---|
| 674 | ! (dimensioned in F) and/or Y has length exceeding NEQ. |
---|
| 675 | ! See the descriptions of NEQ and Y below. |
---|
| 676 | ! |
---|
| 677 | ! If quantities computed in the F routine are needed |
---|
| 678 | ! externally to DLSODE, an extra call to F should be made |
---|
| 679 | ! for this purpose, for consistent and accurate results. |
---|
| 680 | ! If only the derivative dy/dt is needed, use DINTDY |
---|
| 681 | ! instead. |
---|
| 682 | ! |
---|
| 683 | ! NEQ The size of the ODE system (number of first-order |
---|
| 684 | ! ordinary differential equations). Used only for input. |
---|
| 685 | ! NEQ may be decreased, but not increased, during the |
---|
| 686 | ! problem. If NEQ is decreased (with ISTATE = 3 on input), |
---|
| 687 | ! the remaining components of Y should be left undisturbed, |
---|
| 688 | ! if these are to be accessed in F and/or JAC. |
---|
| 689 | ! |
---|
| 690 | ! Normally, NEQ is a scalar, and it is generally referred |
---|
| 691 | ! to as a scalar in this user interface description. |
---|
| 692 | ! However, NEQ may be an array, with NEQ set to the |
---|
| 693 | ! system size. (The DLSODE package accesses only NEQ.) |
---|
| 694 | ! In either case, this parameter is passed as the NEQ |
---|
| 695 | ! argument in all calls to F and JAC. Hence, if it is an |
---|
| 696 | ! array, locations NEQ(2),... may be used to store other |
---|
| 697 | ! integer data and pass it to F and/or JAC. Subroutines |
---|
| 698 | ! F and/or JAC must include NEQ in a DIMENSION statement |
---|
| 699 | ! in that case. |
---|
| 700 | ! |
---|
| 701 | ! Y A real array for the vector of dependent variables, of |
---|
| 702 | ! length NEQ or more. Used for both input and output on |
---|
| 703 | ! the first call (ISTATE = 1), and only for output on |
---|
| 704 | ! other calls. On the first call, Y must contain the |
---|
| 705 | ! vector of initial values. On output, Y contains the |
---|
| 706 | ! computed solution vector, evaluated at T. If desired, |
---|
| 707 | ! the Y array may be used for other purposes between |
---|
| 708 | ! calls to the solver. |
---|
| 709 | ! |
---|
| 710 | ! This array is passed as the Y argument in all calls to F |
---|
| 711 | ! and JAC. Hence its length may exceed NEQ, and locations |
---|
| 712 | ! Y(NEQ+1),... may be used to store other real data and |
---|
| 713 | ! pass it to F and/or JAC. (The DLSODE package accesses |
---|
| 714 | ! only Y(1),...,Y(NEQ).) |
---|
| 715 | ! |
---|
| 716 | ! T The independent variable. On input, T is used only on |
---|
| 717 | ! the first call, as the initial point of the integration. |
---|
| 718 | ! On output, after each call, T is the value at which a |
---|
| 719 | ! computed solution Y is evaluated (usually the same as |
---|
| 720 | ! TOUT). On an error return, T is the farthest point |
---|
| 721 | ! reached. |
---|
| 722 | ! |
---|
| 723 | ! TOUT The next value of T at which a computed solution is |
---|
| 724 | ! desired. Used only for input. |
---|
| 725 | ! |
---|
| 726 | ! When starting the problem (ISTATE = 1), TOUT may be equal |
---|
| 727 | ! to T for one call, then should not equal T for the next |
---|
| 728 | ! call. For the initial T, an input value of TOUT .NE. T |
---|
| 729 | ! is used in order to determine the direction of the |
---|
| 730 | ! integration (i.e., the algebraic sign of the step sizes) |
---|
| 731 | ! and the rough scale of the problem. Integration in |
---|
| 732 | ! either direction (forward or backward in T) is permitted. |
---|
| 733 | ! |
---|
| 734 | ! If ITASK = 2 or 5 (one-step modes), TOUT is ignored |
---|
| 735 | ! after the first call (i.e., the first call with |
---|
| 736 | ! TOUT .NE. T). Otherwise, TOUT is required on every call. |
---|
| 737 | ! |
---|
| 738 | ! If ITASK = 1, 3, or 4, the values of TOUT need not be |
---|
| 739 | ! monotone, but a value of TOUT which backs up is limited |
---|
| 740 | ! to the current internal T interval, whose endpoints are |
---|
| 741 | ! TCUR - HU and TCUR. (See "Optional Outputs" below for |
---|
| 742 | ! TCUR and HU.) |
---|
| 743 | ! |
---|
| 744 | ! |
---|
| 745 | ! ITOL An indicator for the type of error control. See |
---|
| 746 | ! description below under AbsTol. Used only for input. |
---|
| 747 | ! |
---|
| 748 | ! RelTol A relative error tolerance parameter, either a scalar or |
---|
| 749 | ! an array of length NEQ. See description below under |
---|
| 750 | ! AbsTol. Input only. |
---|
| 751 | ! |
---|
| 752 | ! AbsTol An absolute error tolerance parameter, either a scalar or |
---|
| 753 | ! an array of length NEQ. Input only. |
---|
| 754 | ! |
---|
| 755 | ! The input parameters ITOL, RelTol, and AbsTol determine the |
---|
| 756 | ! error control performed by the solver. The solver will |
---|
| 757 | ! control the vector e = (e(i)) of estimated local errors |
---|
| 758 | ! in Y, according to an inequality of the form |
---|
| 759 | ! |
---|
| 760 | ! rms-norm of ( e(i)/EWT(i) ) <= 1, |
---|
| 761 | ! |
---|
| 762 | ! where |
---|
| 763 | ! |
---|
| 764 | ! EWT(i) = RelTol(i)*ABS(Y(i)) + AbsTol(i), |
---|
| 765 | ! |
---|
| 766 | ! and the rms-norm (root-mean-square norm) here is |
---|
| 767 | ! |
---|
| 768 | ! rms-norm(v) = SQRT(sum v(i)**2 / NEQ). |
---|
| 769 | ! |
---|
| 770 | ! Here EWT = (EWT(i)) is a vector of weights which must |
---|
| 771 | ! always be positive, and the values of RelTol and AbsTol |
---|
| 772 | ! should all be nonnegative. The following table gives the |
---|
| 773 | ! types (scalar/array) of RelTol and AbsTol, and the |
---|
| 774 | ! corresponding form of EWT(i). |
---|
| 775 | ! |
---|
| 776 | ! ITOL RelTol AbsTol EWT(i) |
---|
| 777 | ! ---- ------ ------ ----------------------------- |
---|
| 778 | ! 1 scalar scalar RelTol*ABS(Y(i)) + AbsTol |
---|
| 779 | ! 2 scalar array RelTol*ABS(Y(i)) + AbsTol(i) |
---|
| 780 | ! 3 array scalar RelTol(i)*ABS(Y(i)) + AbsTol |
---|
| 781 | ! 4 array array RelTol(i)*ABS(Y(i)) + AbsTol(i) |
---|
| 782 | ! |
---|
| 783 | ! When either of these parameters is a scalar, it need not |
---|
| 784 | ! be dimensioned in the user's calling program. |
---|
| 785 | ! |
---|
| 786 | ! If none of the above choices (with ITOL, RelTol, and AbsTol |
---|
| 787 | ! fixed throughout the problem) is suitable, more general |
---|
| 788 | ! error controls can be obtained by substituting |
---|
| 789 | ! user-supplied routines for the setting of EWT and/or for |
---|
| 790 | ! the norm calculation. See Part 4 below. |
---|
| 791 | ! |
---|
| 792 | ! If global errors are to be estimated by making a repeated |
---|
| 793 | ! run on the same problem with smaller tolerances, then all |
---|
| 794 | ! components of RelTol and AbsTol (i.e., of EWT) should be |
---|
| 795 | ! scaled down uniformly. |
---|
| 796 | ! |
---|
| 797 | ! ITASK An index specifying the task to be performed. Input |
---|
| 798 | ! only. ITASK has the following values and meanings: |
---|
| 799 | ! 1 Normal computation of output values of y(t) at |
---|
| 800 | ! t = TOUT (by overshooting and interpolating). |
---|
| 801 | ! 2 Take one step only and return. |
---|
| 802 | ! 3 Stop at the first internal mesh point at or beyond |
---|
| 803 | ! t = TOUT and return. |
---|
| 804 | ! 4 Normal computation of output values of y(t) at |
---|
| 805 | ! t = TOUT but without overshooting t = TCRIT. TCRIT |
---|
| 806 | ! must be input as RWORK(1). TCRIT may be equal to or |
---|
| 807 | ! beyond TOUT, but not behind it in the direction of |
---|
| 808 | ! integration. This option is useful if the problem |
---|
| 809 | ! has a singularity at or beyond t = TCRIT. |
---|
| 810 | ! 5 Take one step, without passing TCRIT, and return. |
---|
| 811 | ! TCRIT must be input as RWORK(1). |
---|
| 812 | ! |
---|
| 813 | ! Note: If ITASK = 4 or 5 and the solver reaches TCRIT |
---|
| 814 | ! (within roundoff), it will return T = TCRIT (exactly) to |
---|
| 815 | ! indicate this (unless ITASK = 4 and TOUT comes before |
---|
| 816 | ! TCRIT, in which case answers at T = TOUT are returned |
---|
| 817 | ! first). |
---|
| 818 | ! |
---|
| 819 | ! ISTATE An index used for input and output to specify the state |
---|
| 820 | ! of the calculation. |
---|
| 821 | ! |
---|
| 822 | ! On input, the values of ISTATE are as follows: |
---|
| 823 | ! 1 This is the first call for the problem |
---|
| 824 | ! (initializations will be done). See "Note" below. |
---|
| 825 | ! 2 This is not the first call, and the calculation is to |
---|
| 826 | ! continue normally, with no change in any input |
---|
| 827 | ! parameters except possibly TOUT and ITASK. (If ITOL, |
---|
| 828 | ! RelTol, and/or AbsTol are changed between calls with |
---|
| 829 | ! ISTATE = 2, the new values will be used but not |
---|
| 830 | ! tested for legality.) |
---|
| 831 | ! 3 This is not the first call, and the calculation is to |
---|
| 832 | ! continue normally, but with a change in input |
---|
| 833 | ! parameters other than TOUT and ITASK. Changes are |
---|
| 834 | ! allowed in NEQ, ITOL, RelTol, AbsTol, IOPT, LRW, LIW, MF, |
---|
| 835 | ! ML, MU, and any of the optional inputs except H0. |
---|
| 836 | ! (See IWORK description for ML and MU.) |
---|
| 837 | ! |
---|
| 838 | ! Note: A preliminary call with TOUT = T is not counted as |
---|
| 839 | ! a first call here, as no initialization or checking of |
---|
| 840 | ! input is done. (Such a call is sometimes useful for the |
---|
| 841 | ! purpose of outputting the initial conditions.) Thus the |
---|
| 842 | ! first call for which TOUT .NE. T requires ISTATE = 1 on |
---|
| 843 | ! input. |
---|
| 844 | ! |
---|
| 845 | ! On output, ISTATE has the following values and meanings: |
---|
| 846 | ! 1 Nothing was done, as TOUT was equal to T with |
---|
| 847 | ! ISTATE = 1 on input. |
---|
| 848 | ! 2 The integration was performed successfully. |
---|
| 849 | ! -1 An excessive amount of work (more than MXSTEP steps) |
---|
| 850 | ! was done on this call, before completing the |
---|
| 851 | ! requested task, but the integration was otherwise |
---|
| 852 | ! successful as far as T. (MXSTEP is an optional input |
---|
| 853 | ! and is normally 500.) To continue, the user may |
---|
| 854 | ! simply reset ISTATE to a value >1 and call again (the |
---|
| 855 | ! excess work step counter will be reset to 0). In |
---|
| 856 | ! addition, the user may increase MXSTEP to avoid this |
---|
| 857 | ! error return; see "Optional Inputs" below. |
---|
| 858 | ! -2 Too much accuracy was requested for the precision of |
---|
| 859 | ! the machine being used. This was detected before |
---|
| 860 | ! completing the requested task, but the integration |
---|
| 861 | ! was successful as far as T. To continue, the |
---|
| 862 | ! tolerance parameters must be reset, and ISTATE must |
---|
| 863 | ! be set to 3. The optional output TOLSF may be used |
---|
| 864 | ! for this purpose. (Note: If this condition is |
---|
| 865 | ! detected before taking any steps, then an illegal |
---|
| 866 | ! input return (ISTATE = -3) occurs instead.) |
---|
| 867 | ! -3 Illegal input was detected, before taking any |
---|
| 868 | ! integration steps. See written message for details. |
---|
| 869 | ! (Note: If the solver detects an infinite loop of |
---|
| 870 | ! calls to the solver with illegal input, it will cause |
---|
| 871 | ! the run to stop.) |
---|
| 872 | ! -4 There were repeated error-test failures on one |
---|
| 873 | ! attempted step, before completing the requested task, |
---|
| 874 | ! but the integration was successful as far as T. The |
---|
| 875 | ! problem may have a singularity, or the input may be |
---|
| 876 | ! inappropriate. |
---|
| 877 | ! -5 There were repeated convergence-test failures on one |
---|
| 878 | ! attempted step, before completing the requested task, |
---|
| 879 | ! but the integration was successful as far as T. This |
---|
| 880 | ! may be caused by an inaccurate Jacobian matrix, if |
---|
| 881 | ! one is being used. |
---|
| 882 | ! -6 EWT(i) became zero for some i during the integration. |
---|
| 883 | ! Pure relative error control (AbsTol(i)=0.0) was |
---|
| 884 | ! requested on a variable which has now vanished. The |
---|
| 885 | ! integration was successful as far as T. |
---|
| 886 | ! |
---|
| 887 | ! Note: Since the normal output value of ISTATE is 2, it |
---|
| 888 | ! does not need to be reset for normal continuation. Also, |
---|
| 889 | ! since a negative input value of ISTATE will be regarded |
---|
| 890 | ! as illegal, a negative output value requires the user to |
---|
| 891 | ! change it, and possibly other inputs, before calling the |
---|
| 892 | ! solver again. |
---|
| 893 | ! |
---|
| 894 | ! IOPT An integer flag to specify whether any optional inputs |
---|
| 895 | ! are being used on this call. Input only. The optional |
---|
| 896 | ! inputs are listed under a separate heading below. |
---|
| 897 | ! 0 No optional inputs are being used. Default values |
---|
| 898 | ! will be used in all cases. |
---|
| 899 | ! 1 One or more optional inputs are being used. |
---|
| 900 | ! |
---|
| 901 | ! RWORK A real working array (double precision). The length of |
---|
| 902 | ! RWORK must be at least |
---|
| 903 | ! |
---|
| 904 | ! 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM |
---|
| 905 | ! |
---|
| 906 | ! where |
---|
| 907 | ! NYH = the initial value of NEQ, |
---|
| 908 | ! MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a |
---|
| 909 | ! smaller value is given as an optional input), |
---|
| 910 | ! LWM = 0 if MITER = 0, |
---|
| 911 | ! LWM = NEQ**2 + 2 if MITER = 1 or 2, |
---|
| 912 | ! LWM = NEQ + 2 if MITER = 3, and |
---|
| 913 | ! LWM = (2*ML + MU + 1)*NEQ + 2 |
---|
| 914 | ! if MITER = 4 or 5. |
---|
| 915 | ! (See the MF description below for METH and MITER.) |
---|
| 916 | ! |
---|
| 917 | ! Thus if MAXORD has its default value and NEQ is constant, |
---|
| 918 | ! this length is: |
---|
| 919 | ! 20 + 16*NEQ for MF = 10, |
---|
| 920 | ! 22 + 16*NEQ + NEQ**2 for MF = 11 or 12, |
---|
| 921 | ! 22 + 17*NEQ for MF = 13, |
---|
| 922 | ! 22 + 17*NEQ + (2*ML + MU)*NEQ for MF = 14 or 15, |
---|
| 923 | ! 20 + 9*NEQ for MF = 20, |
---|
| 924 | ! 22 + 9*NEQ + NEQ**2 for MF = 21 or 22, |
---|
| 925 | ! 22 + 10*NEQ for MF = 23, |
---|
| 926 | ! 22 + 10*NEQ + (2*ML + MU)*NEQ for MF = 24 or 25. |
---|
| 927 | ! |
---|
| 928 | ! The first 20 words of RWORK are reserved for conditional |
---|
| 929 | ! and optional inputs and optional outputs. |
---|
| 930 | ! |
---|
| 931 | ! The following word in RWORK is a conditional input: |
---|
| 932 | ! RWORK(1) = TCRIT, the critical value of t which the |
---|
| 933 | ! solver is not to overshoot. Required if ITASK |
---|
| 934 | ! is 4 or 5, and ignored otherwise. See ITASK. |
---|
| 935 | ! |
---|
| 936 | ! LRW The length of the array RWORK, as declared by the user. |
---|
| 937 | ! (This will be checked by the solver.) |
---|
| 938 | ! |
---|
| 939 | ! IWORK An integer work array. Its length must be at least |
---|
| 940 | ! 20 if MITER = 0 or 3 (MF = 10, 13, 20, 23), or |
---|
| 941 | ! 20 + NEQ otherwise (MF = 11, 12, 14, 15, 21, 22, 24, 25). |
---|
| 942 | ! (See the MF description below for MITER.) The first few |
---|
| 943 | ! words of IWORK are used for conditional and optional |
---|
| 944 | ! inputs and optional outputs. |
---|
| 945 | ! |
---|
| 946 | ! The following two words in IWORK are conditional inputs: |
---|
| 947 | ! IWORK(1) = ML These are the lower and upper half- |
---|
| 948 | ! IWORK(2) = MU bandwidths, respectively, of the banded |
---|
| 949 | ! Jacobian, excluding the main diagonal. |
---|
| 950 | ! The band is defined by the matrix locations |
---|
| 951 | ! (i,j) with i - ML <= j <= i + MU. ML and MU |
---|
| 952 | ! must satisfy 0 <= ML,MU <= NEQ - 1. These are |
---|
| 953 | ! required if MITER is 4 or 5, and ignored |
---|
| 954 | ! otherwise. ML and MU may in fact be the band |
---|
| 955 | ! parameters for a matrix to which df/dy is only |
---|
| 956 | ! approximately equal. |
---|
| 957 | ! |
---|
| 958 | ! LIW The length of the array IWORK, as declared by the user. |
---|
| 959 | ! (This will be checked by the solver.) |
---|
| 960 | ! |
---|
| 961 | ! Note: The work arrays must not be altered between calls to DLSODE |
---|
| 962 | ! for the same problem, except possibly for the conditional and |
---|
| 963 | ! optional inputs, and except for the last 3*NEQ words of RWORK. |
---|
| 964 | ! The latter space is used for internal scratch space, and so is |
---|
| 965 | ! available for use by the user outside DLSODE between calls, if |
---|
| 966 | ! desired (but not for use by F or JAC). |
---|
| 967 | ! |
---|
| 968 | ! JAC The name of the user-supplied routine (MITER = 1 or 4) to |
---|
| 969 | ! compute the Jacobian matrix, df/dy, as a function of the |
---|
| 970 | ! scalar t and the vector y. (See the MF description below |
---|
| 971 | ! for MITER.) It is to have the form |
---|
| 972 | ! |
---|
| 973 | ! SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD) |
---|
| 974 | ! KPP_REAL T, Y(*), PD(NROWPD,*) |
---|
| 975 | ! |
---|
| 976 | ! where NEQ, T, Y, ML, MU, and NROWPD are input and the |
---|
| 977 | ! array PD is to be loaded with partial derivatives |
---|
| 978 | ! (elements of the Jacobian matrix) on output. PD must be |
---|
| 979 | ! given a first dimension of NROWPD. T and Y have the same |
---|
| 980 | ! meaning as in subroutine F. |
---|
| 981 | ! |
---|
| 982 | ! In the full matrix case (MITER = 1), ML and MU are |
---|
| 983 | ! ignored, and the Jacobian is to be loaded into PD in |
---|
| 984 | ! columnwise manner, with df(i)/dy(j) loaded into PD(i,j). |
---|
| 985 | ! |
---|
| 986 | ! In the band matrix case (MITER = 4), the elements within |
---|
| 987 | ! the band are to be loaded into PD in columnwise manner, |
---|
| 988 | ! with diagonal lines of df/dy loaded into the rows of PD. |
---|
| 989 | ! Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j). ML |
---|
| 990 | ! and MU are the half-bandwidth parameters (see IWORK). |
---|
| 991 | ! The locations in PD in the two triangular areas which |
---|
| 992 | ! correspond to nonexistent matrix elements can be ignored |
---|
| 993 | ! or loaded arbitrarily, as they are overwritten by DLSODE. |
---|
| 994 | ! |
---|
| 995 | ! JAC need not provide df/dy exactly. A crude approximation |
---|
| 996 | ! (possibly with a smaller bandwidth) will do. |
---|
| 997 | ! |
---|
| 998 | ! In either case, PD is preset to zero by the solver, so |
---|
| 999 | ! that only the nonzero elements need be loaded by JAC. |
---|
| 1000 | ! Each call to JAC is preceded by a call to F with the same |
---|
| 1001 | ! arguments NEQ, T, and Y. Thus to gain some efficiency, |
---|
| 1002 | ! intermediate quantities shared by both calculations may |
---|
| 1003 | ! be saved in a user COMMON block by F and not recomputed |
---|
| 1004 | ! by JAC, if desired. Also, JAC may alter the Y array, if |
---|
| 1005 | ! desired. JAC must be declared EXTERNAL in the calling |
---|
| 1006 | ! program. |
---|
| 1007 | ! |
---|
| 1008 | ! Subroutine JAC may access user-defined quantities in |
---|
| 1009 | ! NEQ(2),... and/or in Y(NEQ+1),... if NEQ is an array |
---|
| 1010 | ! (dimensioned in JAC) and/or Y has length exceeding |
---|
| 1011 | ! NEQ. See the descriptions of NEQ and Y above. |
---|
| 1012 | ! |
---|
| 1013 | ! MF The method flag. Used only for input. The legal values |
---|
| 1014 | ! of MF are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, |
---|
| 1015 | ! and 25. MF has decimal digits METH and MITER: |
---|
| 1016 | ! MF = 10*METH + MITER . |
---|
| 1017 | ! |
---|
| 1018 | ! METH indicates the basic linear multistep method: |
---|
| 1019 | ! 1 Implicit Adams method. |
---|
| 1020 | ! 2 Method based on backward differentiation formulas |
---|
| 1021 | ! (BDF's). |
---|
| 1022 | ! |
---|
| 1023 | ! MITER indicates the corrector iteration method: |
---|
| 1024 | ! 0 Functional iteration (no Jacobian matrix is |
---|
| 1025 | ! involved). |
---|
| 1026 | ! 1 Chord iteration with a user-supplied full (NEQ by |
---|
| 1027 | ! NEQ) Jacobian. |
---|
| 1028 | ! 2 Chord iteration with an internally generated |
---|
| 1029 | ! (difference quotient) full Jacobian (using NEQ |
---|
| 1030 | ! extra calls to F per df/dy value). |
---|
| 1031 | ! 3 Chord iteration with an internally generated |
---|
| 1032 | ! diagonal Jacobian approximation (using one extra call |
---|
| 1033 | ! to F per df/dy evaluation). |
---|
| 1034 | ! 4 Chord iteration with a user-supplied banded Jacobian. |
---|
| 1035 | ! 5 Chord iteration with an internally generated banded |
---|
| 1036 | ! Jacobian (using ML + MU + 1 extra calls to F per |
---|
| 1037 | ! df/dy evaluation). |
---|
| 1038 | ! |
---|
| 1039 | ! If MITER = 1 or 4, the user must supply a subroutine JAC |
---|
| 1040 | ! (the name is arbitrary) as described above under JAC. |
---|
| 1041 | ! For other values of MITER, a dummy argument can be used. |
---|
| 1042 | ! |
---|
| 1043 | ! Optional Inputs |
---|
| 1044 | ! --------------- |
---|
| 1045 | ! The following is a list of the optional inputs provided for in the |
---|
| 1046 | ! call sequence. (See also Part 2.) For each such input variable, |
---|
| 1047 | ! this table lists its name as used in this documentation, its |
---|
| 1048 | ! location in the call sequence, its meaning, and the default value. |
---|
| 1049 | ! The use of any of these inputs requires IOPT = 1, and in that case |
---|
| 1050 | ! all of these inputs are examined. A value of zero for any of |
---|
| 1051 | ! these optional inputs will cause the default value to be used. |
---|
| 1052 | ! Thus to use a subset of the optional inputs, simply preload |
---|
| 1053 | ! locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively, |
---|
| 1054 | ! and then set those of interest to nonzero values. |
---|
| 1055 | ! |
---|
| 1056 | ! Name Location Meaning and default value |
---|
| 1057 | ! ------ --------- ----------------------------------------------- |
---|
| 1058 | ! H0 RWORK(5) Step size to be attempted on the first step. |
---|
| 1059 | ! The default value is determined by the solver. |
---|
| 1060 | ! HMAX RWORK(6) Maximum absolute step size allowed. The |
---|
| 1061 | ! default value is infinite. |
---|
| 1062 | ! HMIN RWORK(7) Minimum absolute step size allowed. The |
---|
| 1063 | ! default value is 0. (This lower bound is not |
---|
| 1064 | ! enforced on the final step before reaching |
---|
| 1065 | ! TCRIT when ITASK = 4 or 5.) |
---|
| 1066 | ! MAXORD IWORK(5) Maximum order to be allowed. The default value |
---|
| 1067 | ! is 12 if METH = 1, and 5 if METH = 2. (See the |
---|
| 1068 | ! MF description above for METH.) If MAXORD |
---|
| 1069 | ! exceeds the default value, it will be reduced |
---|
| 1070 | ! to the default value. If MAXORD is changed |
---|
| 1071 | ! during the problem, it may cause the current |
---|
| 1072 | ! order to be reduced. |
---|
| 1073 | ! MXSTEP IWORK(6) Maximum number of (internally defined) steps |
---|
| 1074 | ! allowed during one call to the solver. The |
---|
| 1075 | ! default value is 500. |
---|
| 1076 | ! MXHNIL IWORK(7) Maximum number of messages printed (per |
---|
| 1077 | ! problem) warning that T + H = T on a step |
---|
| 1078 | ! (H = step size). This must be positive to |
---|
| 1079 | ! result in a nondefault value. The default |
---|
| 1080 | ! value is 10. |
---|
| 1081 | ! |
---|
| 1082 | ! Optional Outputs |
---|
| 1083 | ! ---------------- |
---|
| 1084 | ! As optional additional output from DLSODE, the variables listed |
---|
| 1085 | ! below are quantities related to the performance of DLSODE which |
---|
| 1086 | ! are available to the user. These are communicated by way of the |
---|
| 1087 | ! work arrays, but also have internal mnemonic names as shown. |
---|
| 1088 | ! Except where stated otherwise, all of these outputs are defined on |
---|
| 1089 | ! any successful return from DLSODE, and on any return with ISTATE = |
---|
| 1090 | ! -1, -2, -4, -5, or -6. On an illegal input return (ISTATE = -3), |
---|
| 1091 | ! they will be unchanged from their existing values (if any), except |
---|
| 1092 | ! possibly for TOLSF, LENRW, and LENIW. On any error return, |
---|
| 1093 | ! outputs relevant to the error will be defined, as noted below. |
---|
| 1094 | ! |
---|
| 1095 | ! Name Location Meaning |
---|
| 1096 | ! ----- --------- ------------------------------------------------ |
---|
| 1097 | ! HU RWORK(11) Step size in t last used (successfully). |
---|
| 1098 | ! HCUR RWORK(12) Step size to be attempted on the next step. |
---|
| 1099 | ! TCUR RWORK(13) Current value of the independent variable which |
---|
| 1100 | ! the solver has actually reached, i.e., the |
---|
| 1101 | ! current internal mesh point in t. On output, |
---|
| 1102 | ! TCUR will always be at least as far as the |
---|
| 1103 | ! argument T, but may be farther (if interpolation |
---|
| 1104 | ! was done). |
---|
| 1105 | ! TOLSF RWORK(14) Tolerance scale factor, greater than 1.0, |
---|
| 1106 | ! computed when a request for too much accuracy |
---|
| 1107 | ! was detected (ISTATE = -3 if detected at the |
---|
| 1108 | ! start of the problem, ISTATE = -2 otherwise). |
---|
| 1109 | ! If ITOL is left unaltered but RelTol and AbsTol are |
---|
| 1110 | ! uniformly scaled up by a factor of TOLSF for the |
---|
| 1111 | ! next call, then the solver is deemed likely to |
---|
| 1112 | ! succeed. (The user may also ignore TOLSF and |
---|
| 1113 | ! alter the tolerance parameters in any other way |
---|
| 1114 | ! appropriate.) |
---|
| 1115 | ! NST IWORK(11) Number of steps taken for the problem so far. |
---|
| 1116 | ! NFE IWORK(12) Number of F evaluations for the problem so far. |
---|
| 1117 | ! NJE IWORK(13) Number of Jacobian evaluations (and of matrix LU |
---|
| 1118 | ! decompositions) for the problem so far. |
---|
| 1119 | ! NQU IWORK(14) Method order last used (successfully). |
---|
| 1120 | ! NQCUR IWORK(15) Order to be attempted on the next step. |
---|
| 1121 | ! IMXER IWORK(16) Index of the component of largest magnitude in |
---|
| 1122 | ! the weighted local error vector ( e(i)/EWT(i) ), |
---|
| 1123 | ! on an error return with ISTATE = -4 or -5. |
---|
| 1124 | ! LENRW IWORK(17) Length of RWORK actually required. This is |
---|
| 1125 | ! defined on normal returns and on an illegal |
---|
| 1126 | ! input return for insufficient storage. |
---|
| 1127 | ! LENIW IWORK(18) Length of IWORK actually required. This is |
---|
| 1128 | ! defined on normal returns and on an illegal |
---|
| 1129 | ! input return for insufficient storage. |
---|
| 1130 | ! |
---|
| 1131 | ! The following two arrays are segments of the RWORK array which may |
---|
| 1132 | ! also be of interest to the user as optional outputs. For each |
---|
| 1133 | ! array, the table below gives its internal name, its base address |
---|
| 1134 | ! in RWORK, and its description. |
---|
| 1135 | ! |
---|
| 1136 | ! Name Base address Description |
---|
| 1137 | ! ---- ------------ ---------------------------------------------- |
---|
| 1138 | ! YH 21 The Nordsieck history array, of size NYH by |
---|
| 1139 | ! (NQCUR + 1), where NYH is the initial value of |
---|
| 1140 | ! NEQ. For j = 0,1,...,NQCUR, column j + 1 of |
---|
| 1141 | ! YH contains HCUR**j/factorial(j) times the jth |
---|
| 1142 | ! derivative of the interpolating polynomial |
---|
| 1143 | ! currently representing the solution, evaluated |
---|
| 1144 | ! at t = TCUR. |
---|
| 1145 | ! ACOR LENRW-NEQ+1 Array of size NEQ used for the accumulated |
---|
| 1146 | ! corrections on each step, scaled on output to |
---|
| 1147 | ! represent the estimated local error in Y on |
---|
| 1148 | ! the last step. This is the vector e in the |
---|
| 1149 | ! description of the error control. It is |
---|
| 1150 | ! defined only on successful return from DLSODE. |
---|
| 1151 | ! |
---|
| 1152 | ! |
---|
| 1153 | ! Part 2. Other Callable Routines |
---|
| 1154 | ! -------------------------------- |
---|
| 1155 | ! |
---|
| 1156 | ! The following are optional calls which the user may make to gain |
---|
| 1157 | ! additional capabilities in conjunction with DLSODE. |
---|
| 1158 | ! |
---|
| 1159 | ! Form of call Function |
---|
| 1160 | ! ------------------------ ---------------------------------------- |
---|
| 1161 | ! CALL XSETUN(LUN) Set the logical unit number, LUN, for |
---|
| 1162 | ! output of messages from DLSODE, if the |
---|
| 1163 | ! default is not desired. The default |
---|
| 1164 | ! value of LUN is 6. This call may be made |
---|
| 1165 | ! at any time and will take effect |
---|
| 1166 | ! immediately. |
---|
| 1167 | ! CALL XSETF(MFLAG) Set a flag to control the printing of |
---|
| 1168 | ! messages by DLSODE. MFLAG = 0 means do |
---|
| 1169 | ! not print. (Danger: this risks losing |
---|
| 1170 | ! valuable information.) MFLAG = 1 means |
---|
| 1171 | ! print (the default). This call may be |
---|
| 1172 | ! made at any time and will take effect |
---|
| 1173 | ! immediately. |
---|
| 1174 | ! CALL DSRCOM(RSAV,ISAV,JOB) Saves and restores the contents of the |
---|
| 1175 | ! internal COMMON blocks used by DLSODE |
---|
| 1176 | ! (see Part 3 below). RSAV must be a |
---|
| 1177 | ! real array of length 218 or more, and |
---|
| 1178 | ! ISAV must be an integer array of length |
---|
| 1179 | ! 37 or more. JOB = 1 means save COMMON |
---|
| 1180 | ! into RSAV/ISAV. JOB = 2 means restore |
---|
| 1181 | ! COMMON from same. DSRCOM is useful if |
---|
| 1182 | ! one is interrupting a run and restarting |
---|
| 1183 | ! later, or alternating between two or |
---|
| 1184 | ! more problems solved with DLSODE. |
---|
| 1185 | ! CALL DINTDY(,,,,,) Provide derivatives of y, of various |
---|
| 1186 | ! (see below) orders, at a specified point t, if |
---|
| 1187 | ! desired. It may be called only after a |
---|
| 1188 | ! successful return from DLSODE. Detailed |
---|
| 1189 | ! instructions follow. |
---|
| 1190 | ! |
---|
| 1191 | ! Detailed instructions for using DINTDY |
---|
| 1192 | ! -------------------------------------- |
---|
| 1193 | ! The form of the CALL is: |
---|
| 1194 | ! |
---|
| 1195 | ! CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG) |
---|
| 1196 | ! |
---|
| 1197 | ! The input parameters are: |
---|
| 1198 | ! |
---|
| 1199 | ! T Value of independent variable where answers are |
---|
| 1200 | ! desired (normally the same as the T last returned by |
---|
| 1201 | ! DLSODE). For valid results, T must lie between |
---|
| 1202 | ! TCUR - HU and TCUR. (See "Optional Outputs" above |
---|
| 1203 | ! for TCUR and HU.) |
---|
| 1204 | ! K Integer order of the derivative desired. K must |
---|
| 1205 | ! satisfy 0 <= K <= NQCUR, where NQCUR is the current |
---|
| 1206 | ! order (see "Optional Outputs"). The capability |
---|
| 1207 | ! corresponding to K = 0, i.e., computing y(t), is |
---|
| 1208 | ! already provided by DLSODE directly. Since |
---|
| 1209 | ! NQCUR >= 1, the first derivative dy/dt is always |
---|
| 1210 | ! available with DINTDY. |
---|
| 1211 | ! RWORK(21) The base address of the history array YH. |
---|
| 1212 | ! NYH Column length of YH, equal to the initial value of NEQ. |
---|
| 1213 | ! |
---|
| 1214 | ! The output parameters are: |
---|
| 1215 | ! |
---|
| 1216 | ! DKY Real array of length NEQ containing the computed value |
---|
| 1217 | ! of the Kth derivative of y(t). |
---|
| 1218 | ! IFLAG Integer flag, returned as 0 if K and T were legal, |
---|
| 1219 | ! -1 if K was illegal, and -2 if T was illegal. |
---|
| 1220 | ! On an error return, a message is also written. |
---|
| 1221 | ! |
---|
| 1222 | ! |
---|
| 1223 | ! Part 3. Common Blocks |
---|
| 1224 | ! ---------------------- |
---|
| 1225 | ! |
---|
| 1226 | ! If DLSODE is to be used in an overlay situation, the user must |
---|
| 1227 | ! declare, in the primary overlay, the variables in: |
---|
| 1228 | ! (1) the call sequence to DLSODE, |
---|
| 1229 | ! (2) the internal COMMON block /DLS001/, of length 255 |
---|
| 1230 | ! (218 double precision words followed by 37 integer words). |
---|
| 1231 | ! |
---|
| 1232 | ! If DLSODE is used on a system in which the contents of internal |
---|
| 1233 | ! COMMON blocks are not preserved between calls, the user should |
---|
| 1234 | ! declare the above COMMON block in his main program to insure that |
---|
| 1235 | ! its contents are preserved. |
---|
| 1236 | ! |
---|
| 1237 | ! If the solution of a given problem by DLSODE is to be interrupted |
---|
| 1238 | ! and then later continued, as when restarting an interrupted run or |
---|
| 1239 | ! alternating between two or more problems, the user should save, |
---|
| 1240 | ! following the return from the last DLSODE call prior to the |
---|
| 1241 | ! interruption, the contents of the call sequence variables and the |
---|
| 1242 | ! internal COMMON block, and later restore these values before the |
---|
| 1243 | ! next DLSODE call for that problem. In addition, if XSETUN and/or |
---|
| 1244 | ! XSETF was called for non-default handling of error messages, then |
---|
| 1245 | ! these calls must be repeated. To save and restore the COMMON |
---|
| 1246 | ! block, use subroutine DSRCOM (see Part 2 above). |
---|
| 1247 | ! |
---|
| 1248 | ! |
---|
| 1249 | ! Part 4. Optionally Replaceable Solver Routines |
---|
| 1250 | ! ----------------------------------------------- |
---|
| 1251 | ! |
---|
| 1252 | ! Below are descriptions of two routines in the DLSODE package which |
---|
| 1253 | ! relate to the measurement of errors. Either routine can be |
---|
| 1254 | ! replaced by a user-supplied version, if desired. However, since |
---|
| 1255 | ! such a replacement may have a major impact on performance, it |
---|
| 1256 | ! should be done only when absolutely necessary, and only with great |
---|
| 1257 | ! caution. (Note: The means by which the package version of a |
---|
| 1258 | ! routine is superseded by the user's version may be system- |
---|
| 1259 | ! dependent.) |
---|
| 1260 | ! |
---|
| 1261 | ! DEWSET |
---|
| 1262 | ! ------ |
---|
| 1263 | ! The following subroutine is called just before each internal |
---|
| 1264 | ! integration step, and sets the array of error weights, EWT, as |
---|
| 1265 | ! described under ITOL/RelTol/AbsTol above: |
---|
| 1266 | ! |
---|
| 1267 | ! SUBROUTINE DEWSET (NEQ, ITOL, RelTol, AbsTol, YCUR, EWT) |
---|
| 1268 | ! |
---|
| 1269 | ! where NEQ, ITOL, RelTol, and AbsTol are as in the DLSODE call |
---|
| 1270 | ! sequence, YCUR contains the current dependent variable vector, |
---|
| 1271 | ! and EWT is the array of weights set by DEWSET. |
---|
| 1272 | ! |
---|
| 1273 | ! If the user supplies this subroutine, it must return in EWT(i) |
---|
| 1274 | ! (i = 1,...,NEQ) a positive quantity suitable for comparing errors |
---|
| 1275 | ! in Y(i) to. The EWT array returned by DEWSET is passed to the |
---|
| 1276 | ! DVNORM routine (see below), and also used by DLSODE in the |
---|
| 1277 | ! computation of the optional output IMXER, the diagonal Jacobian |
---|
| 1278 | ! approximation, and the increments for difference quotient |
---|
| 1279 | ! Jacobians. |
---|
| 1280 | ! |
---|
| 1281 | ! In the user-supplied version of DEWSET, it may be desirable to use |
---|
| 1282 | ! the current values of derivatives of y. Derivatives up to order NQ |
---|
| 1283 | ! are available from the history array YH, described above under |
---|
| 1284 | ! optional outputs. In DEWSET, YH is identical to the YCUR array, |
---|
| 1285 | ! extended to NQ + 1 columns with a column length of NYH and scale |
---|
| 1286 | ! factors of H**j/factorial(j). On the first call for the problem, |
---|
| 1287 | ! given by NST = 0, NQ is 1 and H is temporarily set to 1.0. |
---|
| 1288 | ! NYH is the initial value of NEQ. The quantities NQ, H, and NST |
---|
| 1289 | ! can be obtained by including in SEWSET the statements: |
---|
| 1290 | ! KPP_REAL RLS |
---|
| 1291 | ! COMMON /DLS001/ RLS(218),ILS(37) |
---|
| 1292 | ! NQ = ILS(33) |
---|
| 1293 | ! NST = ILS(34) |
---|
| 1294 | ! H = RLS(212) |
---|
| 1295 | ! Thus, for example, the current value of dy/dt can be obtained as |
---|
| 1296 | ! YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is unnecessary |
---|
| 1297 | ! when NST = 0). |
---|
| 1298 | ! |
---|
| 1299 | ! DVNORM |
---|
| 1300 | ! ------ |
---|
| 1301 | ! DVNORM is a real function routine which computes the weighted |
---|
| 1302 | ! root-mean-square norm of a vector v: |
---|
| 1303 | ! |
---|
| 1304 | ! d = DVNORM (n, v, w) |
---|
| 1305 | ! |
---|
| 1306 | ! where: |
---|
| 1307 | ! n = the length of the vector, |
---|
| 1308 | ! v = real array of length n containing the vector, |
---|
| 1309 | ! w = real array of length n containing weights, |
---|
| 1310 | ! d = SQRT( (1/n) * sum(v(i)*w(i))**2 ). |
---|
| 1311 | ! |
---|
| 1312 | ! DVNORM is called with n = NEQ and with w(i) = 1.0/EWT(i), where |
---|
| 1313 | ! EWT is as set by subroutine DEWSET. |
---|
| 1314 | ! |
---|
| 1315 | ! If the user supplies this function, it should return a nonnegative |
---|
| 1316 | ! value of DVNORM suitable for use in the error control in DLSODE. |
---|
| 1317 | ! None of the arguments should be altered by DVNORM. For example, a |
---|
| 1318 | ! user-supplied DVNORM routine might: |
---|
| 1319 | ! - Substitute a max-norm of (v(i)*w(i)) for the rms-norm, or |
---|
| 1320 | ! - Ignore some components of v in the norm, with the effect of |
---|
| 1321 | ! suppressing the error control on those components of Y. |
---|
| 1322 | ! --------------------------------------------------------------------- |
---|
| 1323 | !***ROUTINES CALLED DEWSET, DINTDY, DUMACH, DSTODE, DVNORM, XERRWD |
---|
| 1324 | !***COMMON BLOCKS DLS001 |
---|
| 1325 | !***REVISION HISTORY (YYYYMMDD) |
---|
| 1326 | ! 19791129 DATE WRITTEN |
---|
| 1327 | ! 19791213 Minor changes to declarations; DELP init. in STODE. |
---|
| 1328 | ! 19800118 Treat NEQ as array; integer declarations added throughout; |
---|
| 1329 | ! minor changes to prologue. |
---|
| 1330 | ! 19800306 Corrected TESCO(1,NQP1) setting in CFODE. |
---|
| 1331 | ! 19800519 Corrected access of YH on forced order reduction; |
---|
| 1332 | ! numerous corrections to prologues and other comments. |
---|
| 1333 | ! 19800617 In main driver, added loading of SQRT(UROUND) in RWORK; |
---|
| 1334 | ! minor corrections to main prologue. |
---|
| 1335 | ! 19800923 Added zero initialization of HU and NQU. |
---|
| 1336 | ! 19801218 Revised XERRWD routine; minor corrections to main prologue. |
---|
| 1337 | ! 19810401 Minor changes to comments and an error message. |
---|
| 1338 | ! 19810814 Numerous revisions: replaced EWT by 1/EWT; used flags |
---|
| 1339 | ! JCUR, ICF, IERPJ, IERSL between STODE and subordinates; |
---|
| 1340 | ! added tuning parameters CCMAX, MAXCOR, MSBP, MXNCF; |
---|
| 1341 | ! reorganized returns from STODE; reorganized type decls.; |
---|
| 1342 | ! fixed message length in XERRWD; changed default LUNIT to 6; |
---|
| 1343 | ! changed Common lengths; changed comments throughout. |
---|
| 1344 | ! 19870330 Major update by ACH: corrected comments throughout; |
---|
| 1345 | ! removed TRET from Common; rewrote EWSET with 4 loops; |
---|
| 1346 | ! fixed t test in INTDY; added Cray directives in STODE; |
---|
| 1347 | ! in STODE, fixed DELP init. and logic around PJAC call; |
---|
| 1348 | ! combined routines to save/restore Common; |
---|
| 1349 | ! passed LEVEL = 0 in error message calls (except run abort). |
---|
| 1350 | ! 19890426 Modified prologue to SLATEC/LDOC format. (FNF) |
---|
| 1351 | ! 19890501 Many improvements to prologue. (FNF) |
---|
| 1352 | ! 19890503 A few final corrections to prologue. (FNF) |
---|
| 1353 | ! 19890504 Minor cosmetic changes. (FNF) |
---|
| 1354 | ! 19890510 Corrected description of Y in Arguments section. (FNF) |
---|
| 1355 | ! 19890517 Minor corrections to prologue. (FNF) |
---|
| 1356 | ! 19920514 Updated with prologue edited 891025 by G. Shaw for manual. |
---|
| 1357 | ! 19920515 Converted source lines to upper case. (FNF) |
---|
| 1358 | ! 19920603 Revised XERRWD calls using mixed upper-lower case. (ACH) |
---|
| 1359 | ! 19920616 Revised prologue comment regarding CFT. (ACH) |
---|
| 1360 | ! 19921116 Revised prologue comments regarding Common. (ACH). |
---|
| 1361 | ! 19930326 Added comment about non-reentrancy. (FNF) |
---|
| 1362 | ! 19930723 Changed D1MACH to DUMACH. (FNF) |
---|
| 1363 | ! 19930801 Removed ILLIN and NTREP from Common (affects driver logic); |
---|
| 1364 | ! minor changes to prologue and internal comments; |
---|
| 1365 | ! changed Hollerith strings to quoted strings; |
---|
| 1366 | ! changed internal comments to mixed case; |
---|
| 1367 | ! replaced XERRWD with new version using character type; |
---|
| 1368 | ! changed dummy dimensions from 1 to *. (ACH) |
---|
| 1369 | ! 19930809 Changed to generic intrinsic names; changed names of |
---|
| 1370 | ! subprograms and Common blocks to DLSODE etc. (ACH) |
---|
| 1371 | ! 19930929 Eliminated use of REAL intrinsic; other minor changes. (ACH) |
---|
| 1372 | ! 20010412 Removed all 'own' variables from Common block /DLS001/ |
---|
| 1373 | ! (affects declarations in 6 routines). (ACH) |
---|
| 1374 | ! 20010509 Minor corrections to prologue. (ACH) |
---|
| 1375 | ! 20031105 Restored 'own' variables to Common block /DLS001/, to |
---|
| 1376 | ! enable interrupt/restart feature. (ACH) |
---|
| 1377 | ! 20031112 Added SAVE statements for data-loaded constants. |
---|
| 1378 | ! |
---|
| 1379 | !***END PROLOGUE DLSODE |
---|
| 1380 | ! |
---|
| 1381 | !*Internal Notes: |
---|
| 1382 | ! |
---|
| 1383 | ! Other Routines in the DLSODE Package. |
---|
| 1384 | ! |
---|
| 1385 | ! In addition to Subroutine DLSODE, the DLSODE package includes the |
---|
| 1386 | ! following subroutines and function routines: |
---|
| 1387 | ! DINTDY computes an interpolated value of the y vector at t = TOUT. |
---|
| 1388 | ! DSTODE is the core integrator, which does one step of the |
---|
| 1389 | ! integration and the associated error control. |
---|
| 1390 | ! DCFODE sets all method coefficients and test constants. |
---|
| 1391 | ! DPREPJ computes and preprocesses the Jacobian matrix J = df/dy |
---|
| 1392 | ! and the Newton iteration matrix P = I - h*l0*J. |
---|
| 1393 | ! DSOLSY manages solution of linear system in chord iteration. |
---|
| 1394 | ! DEWSET sets the error weight vector EWT before each step. |
---|
| 1395 | ! DVNORM computes the weighted R.M.S. norm of a vector. |
---|
| 1396 | ! DSRCOM is a user-callable routine to save and restore |
---|
| 1397 | ! the contents of the internal Common block. |
---|
| 1398 | ! DGEFA and DGESL are routines from LINPACK for solving full |
---|
| 1399 | ! systems of linear algebraic equations. |
---|
| 1400 | ! DGBFA and DGBSL are routines from LINPACK for solving banded |
---|
| 1401 | ! linear systems. |
---|
| 1402 | ! DUMACH computes the unit roundoff in a machine-independent manner. |
---|
| 1403 | ! XERRWD, XSETUN, XSETF, IXSAV, IUMACH handle the printing of all |
---|
| 1404 | ! error messages and warnings. XERRWD is machine-dependent. |
---|
| 1405 | ! Note: DVNORM, DUMACH, IXSAV, and IUMACH are function routines. |
---|
| 1406 | ! All the others are subroutines. |
---|
| 1407 | ! |
---|
| 1408 | !**End |
---|
| 1409 | ! |
---|
| 1410 | ! Declare externals. |
---|
| 1411 | ! Note: they are now internal |
---|
| 1412 | !EXTERNAL DPREPJ, DSOLSY |
---|
| 1413 | !KPP_REAL DUMACH, DVNORM |
---|
| 1414 | ! |
---|
| 1415 | ! Declare all other variables. |
---|
| 1416 | INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH, IOWNS, & |
---|
| 1417 | ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, & |
---|
| 1418 | LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, & |
---|
| 1419 | MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
---|
| 1420 | INTEGER I, I1, I2, IFLAG, IMXER, KGO, LF0, & |
---|
| 1421 | LENIW, LENRW, LENWM, ML, MORD(2), MU, MXHNL0, MXSTP0 |
---|
| 1422 | KPP_REAL ROWNS, & |
---|
| 1423 | CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND |
---|
| 1424 | KPP_REAL AbsTolI, AYI, BIG, EWTI, H0, HMAX, HMX, RH, RelTolI, & |
---|
| 1425 | TCRIT, TDIST, TNEXT, TOL, TOLSF, TP, SIZE, SUM, W0 |
---|
| 1426 | |
---|
| 1427 | LOGICAL IHIT |
---|
| 1428 | CHARACTER*80 MSG |
---|
| 1429 | SAVE MORD, MXSTP0, MXHNL0 |
---|
| 1430 | !----------------------------------------------------------------------- |
---|
| 1431 | ! The following internal Common block contains |
---|
| 1432 | ! (a) variables which are local to any subroutine but whose values must |
---|
| 1433 | ! be preserved between calls to the routine ("own" variables), and |
---|
| 1434 | ! (b) variables which are communicated between subroutines. |
---|
| 1435 | ! The block DLS001 is declared in subroutines DLSODE, DINTDY, DSTODE, |
---|
| 1436 | ! DPREPJ, and DSOLSY. |
---|
| 1437 | ! Groups of variables are replaced by dummy arrays in the Common |
---|
| 1438 | ! declarations in routines where those variables are not used. |
---|
| 1439 | !----------------------------------------------------------------------- |
---|
| 1440 | COMMON /DLS001/ ROWNS(209), & |
---|
| 1441 | CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, & |
---|
| 1442 | INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH, IOWNS(6), & |
---|
| 1443 | ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, & |
---|
| 1444 | LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, & |
---|
| 1445 | MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
---|
| 1446 | ! |
---|
| 1447 | DATA MORD(1),MORD(2)/12,5/, MXSTP0/500/, MXHNL0/10/ |
---|
| 1448 | !----------------------------------------------------------------------- |
---|
| 1449 | ! Block A. |
---|
| 1450 | ! This code block is executed on every call. |
---|
| 1451 | ! It tests ISTATE and ITASK for legality and branches appropriately. |
---|
| 1452 | ! If ISTATE .GT. 1 but the flag INIT shows that initialization has |
---|
| 1453 | ! not yet been done, an error return occurs. |
---|
| 1454 | ! If ISTATE = 1 and TOUT = T, return immediately. |
---|
| 1455 | !----------------------------------------------------------------------- |
---|
| 1456 | ! |
---|
| 1457 | !***FIRST EXECUTABLE STATEMENT DLSODE |
---|
| 1458 | IF (ISTATE .LT. 1 .OR. ISTATE .GT. 3) GO TO 601 |
---|
| 1459 | IF (ITASK .LT. 1 .OR. ITASK .GT. 5) GO TO 602 |
---|
| 1460 | IF (ISTATE .EQ. 1) GO TO 10 |
---|
| 1461 | IF (INIT .EQ. 0) GO TO 603 |
---|
| 1462 | IF (ISTATE .EQ. 2) GO TO 200 |
---|
| 1463 | GO TO 20 |
---|
| 1464 | 10 INIT = 0 |
---|
| 1465 | IF (TOUT .EQ. T) RETURN |
---|
| 1466 | !----------------------------------------------------------------------- |
---|
| 1467 | ! Block B. |
---|
| 1468 | ! The next code block is executed for the initial call (ISTATE = 1), |
---|
| 1469 | ! or for a continuation call with parameter changes (ISTATE = 3). |
---|
| 1470 | ! It contains checking of all inputs and various initializations. |
---|
| 1471 | ! |
---|
| 1472 | ! First check legality of the non-optional inputs NEQ, ITOL, IOPT, |
---|
| 1473 | ! MF, ML, and MU. |
---|
| 1474 | !----------------------------------------------------------------------- |
---|
| 1475 | 20 IF (NEQ .LE. 0) GO TO 604 |
---|
| 1476 | IF (ISTATE .EQ. 1) GO TO 25 |
---|
| 1477 | IF (NEQ .GT. N) GO TO 605 |
---|
| 1478 | 25 N = NEQ |
---|
| 1479 | IF (ITOL .LT. 1 .OR. ITOL .GT. 4) GO TO 606 |
---|
| 1480 | IF (IOPT .LT. 0 .OR. IOPT .GT. 1) GO TO 607 |
---|
| 1481 | METH = MF/10 |
---|
| 1482 | MITER = MF - 10*METH |
---|
| 1483 | IF (METH .LT. 1 .OR. METH .GT. 2) GO TO 608 |
---|
| 1484 | IF (MITER .LT. 0 .OR. MITER .GT. 5) GO TO 608 |
---|
| 1485 | IF (MITER .LE. 3) GO TO 30 |
---|
| 1486 | ML = IWORK(1) |
---|
| 1487 | MU = IWORK(2) |
---|
| 1488 | IF (ML .LT. 0 .OR. ML .GE. N) GO TO 609 |
---|
| 1489 | IF (MU .LT. 0 .OR. MU .GE. N) GO TO 610 |
---|
| 1490 | 30 CONTINUE |
---|
| 1491 | ! Next process and check the optional inputs. -------------------------- |
---|
| 1492 | IF (IOPT .EQ. 1) GO TO 40 |
---|
| 1493 | MAXORD = MORD(METH) |
---|
| 1494 | MXSTEP = MXSTP0 |
---|
| 1495 | MXHNIL = MXHNL0 |
---|
| 1496 | IF (ISTATE .EQ. 1) H0 = 0.0D0 |
---|
| 1497 | HMXI = 0.0D0 |
---|
| 1498 | HMIN = 0.0D0 |
---|
| 1499 | GO TO 60 |
---|
| 1500 | 40 MAXORD = IWORK(5) |
---|
| 1501 | IF (MAXORD .LT. 0) GO TO 611 |
---|
| 1502 | IF (MAXORD .EQ. 0) MAXORD = 100 |
---|
| 1503 | MAXORD = MIN(MAXORD,MORD(METH)) |
---|
| 1504 | MXSTEP = IWORK(6) |
---|
| 1505 | IF (MXSTEP .LT. 0) GO TO 612 |
---|
| 1506 | IF (MXSTEP .EQ. 0) MXSTEP = MXSTP0 |
---|
| 1507 | MXHNIL = IWORK(7) |
---|
| 1508 | IF (MXHNIL .LT. 0) GO TO 613 |
---|
| 1509 | IF (MXHNIL .EQ. 0) MXHNIL = MXHNL0 |
---|
| 1510 | IF (ISTATE .NE. 1) GO TO 50 |
---|
| 1511 | H0 = RWORK(5) |
---|
| 1512 | IF ((TOUT - T)*H0 .LT. 0.0D0) GO TO 614 |
---|
| 1513 | 50 HMAX = RWORK(6) |
---|
| 1514 | IF (HMAX .LT. 0.0D0) GO TO 615 |
---|
| 1515 | HMXI = 0.0D0 |
---|
| 1516 | IF (HMAX .GT. 0.0D0) HMXI = 1.0D0/HMAX |
---|
| 1517 | HMIN = RWORK(7) |
---|
| 1518 | IF (HMIN .LT. 0.0D0) GO TO 616 |
---|
| 1519 | !----------------------------------------------------------------------- |
---|
| 1520 | ! Set work array pointers and check lengths LRW and LIW. |
---|
| 1521 | ! Pointers to segments of RWORK and IWORK are named by prefixing L to |
---|
| 1522 | ! the name of the segment. E.g., the segment YH starts at RWORK(LYH). |
---|
| 1523 | ! Segments of RWORK (in order) are denoted YH, WM, EWT, SAVF, ACOR. |
---|
| 1524 | !----------------------------------------------------------------------- |
---|
| 1525 | 60 LYH = 21 |
---|
| 1526 | IF (ISTATE .EQ. 1) NYH = N |
---|
| 1527 | LWM = LYH + (MAXORD + 1)*NYH |
---|
| 1528 | IF (MITER .EQ. 0) LENWM = 0 |
---|
| 1529 | IF (MITER .EQ. 1 .OR. MITER .EQ. 2) LENWM = N*N + 2 |
---|
| 1530 | IF (MITER .EQ. 3) LENWM = N + 2 |
---|
| 1531 | IF (MITER .GE. 4) LENWM = (2*ML + MU + 1)*N + 2 |
---|
| 1532 | LEWT = LWM + LENWM |
---|
| 1533 | LSAVF = LEWT + N |
---|
| 1534 | LACOR = LSAVF + N |
---|
| 1535 | LENRW = LACOR + N - 1 |
---|
| 1536 | IWORK(17) = LENRW |
---|
| 1537 | LIWM = 1 |
---|
| 1538 | LENIW = 20 + N |
---|
| 1539 | IF (MITER .EQ. 0 .OR. MITER .EQ. 3) LENIW = 20 |
---|
| 1540 | IWORK(18) = LENIW |
---|
| 1541 | IF (LENRW .GT. LRW) GO TO 617 |
---|
| 1542 | IF (LENIW .GT. LIW) GO TO 618 |
---|
| 1543 | ! Check RelTol and AbsTol for legality. ------------------------------------ |
---|
| 1544 | RelTolI = RelTol(1) |
---|
| 1545 | AbsTolI = AbsTol(1) |
---|
| 1546 | DO 70 I = 1,N |
---|
| 1547 | IF (ITOL .GE. 3) RelTolI = RelTol(I) |
---|
| 1548 | IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) AbsTolI = AbsTol(I) |
---|
| 1549 | IF (RelTolI .LT. 0.0D0) GO TO 619 |
---|
| 1550 | IF (AbsTolI .LT. 0.0D0) GO TO 620 |
---|
| 1551 | 70 CONTINUE |
---|
| 1552 | IF (ISTATE .EQ. 1) GO TO 100 |
---|
| 1553 | ! If ISTATE = 3, set flag to signal parameter changes to DSTODE. ------- |
---|
| 1554 | JSTART = -1 |
---|
| 1555 | IF (NQ .LE. MAXORD) GO TO 90 |
---|
| 1556 | ! MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into SAVF. --------- |
---|
| 1557 | DO 80 I = 1,N |
---|
| 1558 | 80 RWORK(I+LSAVF-1) = RWORK(I+LWM-1) |
---|
| 1559 | ! Reload WM(1) = RWORK(LWM), since LWM may have changed. --------------- |
---|
| 1560 | 90 IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND) |
---|
| 1561 | IF (N .EQ. NYH) GO TO 200 |
---|
| 1562 | ! NEQ was reduced. Zero part of YH to avoid undefined references. ----- |
---|
| 1563 | I1 = LYH + L*NYH |
---|
| 1564 | I2 = LYH + (MAXORD + 1)*NYH - 1 |
---|
| 1565 | IF (I1 .GT. I2) GO TO 200 |
---|
| 1566 | DO 95 I = I1,I2 |
---|
| 1567 | 95 RWORK(I) = 0.0D0 |
---|
| 1568 | GO TO 200 |
---|
| 1569 | !----------------------------------------------------------------------- |
---|
| 1570 | ! Block C. |
---|
| 1571 | ! The next block is for the initial call only (ISTATE = 1). |
---|
| 1572 | ! It contains all remaining initializations, the initial call to F, |
---|
| 1573 | ! and the calculation of the initial step size. |
---|
| 1574 | ! The error weights in EWT are inverted after being loaded. |
---|
| 1575 | !----------------------------------------------------------------------- |
---|
| 1576 | 100 UROUND = DUMACH() |
---|
| 1577 | TN = T |
---|
| 1578 | IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 110 |
---|
| 1579 | TCRIT = RWORK(1) |
---|
| 1580 | IF ((TCRIT - TOUT)*(TOUT - T) .LT. 0.0D0) GO TO 625 |
---|
| 1581 | IF (H0 .NE. 0.0D0 .AND. (T + H0 - TCRIT)*H0 .GT. 0.0D0) & |
---|
| 1582 | H0 = TCRIT - T |
---|
| 1583 | 110 JSTART = 0 |
---|
| 1584 | IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND) |
---|
| 1585 | NHNIL = 0 |
---|
| 1586 | NST = 0 |
---|
| 1587 | NJE = 0 |
---|
| 1588 | NSLAST = 0 |
---|
| 1589 | HU = 0.0D0 |
---|
| 1590 | NQU = 0 |
---|
| 1591 | CCMAX = 0.3D0 |
---|
| 1592 | MAXCOR = 3 |
---|
| 1593 | MSBP = 20 |
---|
| 1594 | MXNCF = 10 |
---|
| 1595 | ! Initial call to F. (LF0 points to YH(*,2).) ------------------------- |
---|
| 1596 | LF0 = LYH + NYH |
---|
| 1597 | CALL F (NEQ, T, Y, RWORK(LF0)) |
---|
| 1598 | NFE = 1 |
---|
| 1599 | ! Load the initial value vector in YH. --------------------------------- |
---|
| 1600 | DO 115 I = 1,N |
---|
| 1601 | 115 RWORK(I+LYH-1) = Y(I) |
---|
| 1602 | ! Load and invert the EWT array. (H is temporarily set to 1.0.) ------- |
---|
| 1603 | NQ = 1 |
---|
| 1604 | H = 1.0D0 |
---|
| 1605 | CALL DEWSET (N, ITOL, RelTol, AbsTol, RWORK(LYH), RWORK(LEWT)) |
---|
| 1606 | DO 120 I = 1,N |
---|
| 1607 | IF (RWORK(I+LEWT-1) .LE. 0.0D0) GO TO 621 |
---|
| 1608 | 120 RWORK(I+LEWT-1) = 1.0D0/RWORK(I+LEWT-1) |
---|
| 1609 | !----------------------------------------------------------------------- |
---|
| 1610 | ! The coding below computes the step size, H0, to be attempted on the |
---|
| 1611 | ! first step, unless the user has supplied a value for this. |
---|
| 1612 | ! First check that TOUT - T differs significantly from zero. |
---|
| 1613 | ! A scalar tolerance quantity TOL is computed, as MAX(RelTol(I)) |
---|
| 1614 | ! if this is positive, or MAX(AbsTol(I)/ABS(Y(I))) otherwise, adjusted |
---|
| 1615 | ! so as to be between 100*UROUND and 1.0E-3. |
---|
| 1616 | ! Then the computed value H0 is given by.. |
---|
| 1617 | ! NEQ |
---|
| 1618 | ! H0**2 = TOL / ( w0**-2 + (1/NEQ) * SUM ( f(i)/ywt(i) )**2 ) |
---|
| 1619 | ! 1 |
---|
| 1620 | ! where w0 = MAX ( ABS(T), ABS(TOUT) ), |
---|
| 1621 | ! f(i) = i-th component of initial value of f, |
---|
| 1622 | ! ywt(i) = EWT(i)/TOL (a weight for y(i)). |
---|
| 1623 | ! The sign of H0 is inferred from the initial values of TOUT and T. |
---|
| 1624 | !----------------------------------------------------------------------- |
---|
| 1625 | IF (H0 .NE. 0.0D0) GO TO 180 |
---|
| 1626 | TDIST = ABS(TOUT - T) |
---|
| 1627 | W0 = MAX(ABS(T),ABS(TOUT)) |
---|
| 1628 | IF (TDIST .LT. 2.0D0*UROUND*W0) GO TO 622 |
---|
| 1629 | TOL = RelTol(1) |
---|
| 1630 | IF (ITOL .LE. 2) GO TO 140 |
---|
| 1631 | DO 130 I = 1,N |
---|
| 1632 | 130 TOL = MAX(TOL,RelTol(I)) |
---|
| 1633 | 140 IF (TOL .GT. 0.0D0) GO TO 160 |
---|
| 1634 | AbsTolI = AbsTol(1) |
---|
| 1635 | DO 150 I = 1,N |
---|
| 1636 | IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) AbsTolI = AbsTol(I) |
---|
| 1637 | AYI = ABS(Y(I)) |
---|
| 1638 | IF (AYI .NE. 0.0D0) TOL = MAX(TOL,AbsTolI/AYI) |
---|
| 1639 | 150 CONTINUE |
---|
| 1640 | 160 TOL = MAX(TOL,100.0D0*UROUND) |
---|
| 1641 | TOL = MIN(TOL,0.001D0) |
---|
| 1642 | SUM = DVNORM (N, RWORK(LF0), RWORK(LEWT)) |
---|
| 1643 | SUM = 1.0D0/(TOL*W0*W0) + TOL*SUM**2 |
---|
| 1644 | H0 = 1.0D0/SQRT(SUM) |
---|
| 1645 | H0 = MIN(H0,TDIST) |
---|
| 1646 | H0 = SIGN(H0,TOUT-T) |
---|
| 1647 | ! Adjust H0 if necessary to meet HMAX bound. --------------------------- |
---|
| 1648 | 180 RH = ABS(H0)*HMXI |
---|
| 1649 | IF (RH .GT. 1.0D0) H0 = H0/RH |
---|
| 1650 | ! Load H with H0 and scale YH(*,2) by H0. ------------------------------ |
---|
| 1651 | H = H0 |
---|
| 1652 | DO 190 I = 1,N |
---|
| 1653 | 190 RWORK(I+LF0-1) = H0*RWORK(I+LF0-1) |
---|
| 1654 | GO TO 270 |
---|
| 1655 | !----------------------------------------------------------------------- |
---|
| 1656 | ! Block D. |
---|
| 1657 | ! The next code block is for continuation calls only (ISTATE = 2 or 3) |
---|
| 1658 | ! and is to check stop conditions before taking a step. |
---|
| 1659 | !----------------------------------------------------------------------- |
---|
| 1660 | 200 NSLAST = NST |
---|
| 1661 | GO TO (210, 250, 220, 230, 240), ITASK |
---|
| 1662 | 210 IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250 |
---|
| 1663 | CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) |
---|
| 1664 | IF (IFLAG .NE. 0) GO TO 627 |
---|
| 1665 | T = TOUT |
---|
| 1666 | GO TO 420 |
---|
| 1667 | 220 TP = TN - HU*(1.0D0 + 100.0D0*UROUND) |
---|
| 1668 | IF ((TP - TOUT)*H .GT. 0.0D0) GO TO 623 |
---|
| 1669 | IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250 |
---|
| 1670 | GO TO 400 |
---|
| 1671 | 230 TCRIT = RWORK(1) |
---|
| 1672 | IF ((TN - TCRIT)*H .GT. 0.0D0) GO TO 624 |
---|
| 1673 | IF ((TCRIT - TOUT)*H .LT. 0.0D0) GO TO 625 |
---|
| 1674 | IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 245 |
---|
| 1675 | CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) |
---|
| 1676 | IF (IFLAG .NE. 0) GO TO 627 |
---|
| 1677 | T = TOUT |
---|
| 1678 | GO TO 420 |
---|
| 1679 | 240 TCRIT = RWORK(1) |
---|
| 1680 | IF ((TN - TCRIT)*H .GT. 0.0D0) GO TO 624 |
---|
| 1681 | 245 HMX = ABS(TN) + ABS(H) |
---|
| 1682 | IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX |
---|
| 1683 | IF (IHIT) GO TO 400 |
---|
| 1684 | TNEXT = TN + H*(1.0D0 + 4.0D0*UROUND) |
---|
| 1685 | IF ((TNEXT - TCRIT)*H .LE. 0.0D0) GO TO 250 |
---|
| 1686 | H = (TCRIT - TN)*(1.0D0 - 4.0D0*UROUND) |
---|
| 1687 | IF (ISTATE .EQ. 2) JSTART = -2 |
---|
| 1688 | !----------------------------------------------------------------------- |
---|
| 1689 | ! Block E. |
---|
| 1690 | ! The next block is normally executed for all calls and contains |
---|
| 1691 | ! the call to the one-step core integrator DSTODE. |
---|
| 1692 | ! |
---|
| 1693 | ! This is a looping point for the integration steps. |
---|
| 1694 | ! |
---|
| 1695 | ! First check for too many steps being taken, update EWT (if not at |
---|
| 1696 | ! start of problem), check for too much accuracy being requested, and |
---|
| 1697 | ! check for H below the roundoff level in T. |
---|
| 1698 | !----------------------------------------------------------------------- |
---|
| 1699 | 250 CONTINUE |
---|
| 1700 | IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500 |
---|
| 1701 | CALL DEWSET (N, ITOL, RelTol, AbsTol, RWORK(LYH), RWORK(LEWT)) |
---|
| 1702 | DO 260 I = 1,N |
---|
| 1703 | IF (RWORK(I+LEWT-1) .LE. 0.0D0) GO TO 510 |
---|
| 1704 | 260 RWORK(I+LEWT-1) = 1.0D0/RWORK(I+LEWT-1) |
---|
| 1705 | 270 TOLSF = UROUND*DVNORM (N, RWORK(LYH), RWORK(LEWT)) |
---|
| 1706 | IF (TOLSF .LE. 1.0D0) GO TO 280 |
---|
| 1707 | TOLSF = TOLSF*2.0D0 |
---|
| 1708 | IF (NST .EQ. 0) GO TO 626 |
---|
| 1709 | GO TO 520 |
---|
| 1710 | 280 IF ((TN + H) .NE. TN) GO TO 290 |
---|
| 1711 | NHNIL = NHNIL + 1 |
---|
| 1712 | IF (NHNIL .GT. MXHNIL) GO TO 290 |
---|
| 1713 | MSG = 'DLSODE- Warning..internal T (=R1) and H (=R2) are' |
---|
| 1714 | CALL XERRWD (MSG, 50, 101, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1715 | MSG=' such that in the machine, T + H = T on the next step ' |
---|
| 1716 | CALL XERRWD (MSG, 60, 101, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1717 | MSG = ' (H = step size). Solver will continue anyway' |
---|
| 1718 | CALL XERRWD (MSG, 50, 101, 0, 0, 0, 0, 2, TN, H) |
---|
| 1719 | IF (NHNIL .LT. MXHNIL) GO TO 290 |
---|
| 1720 | MSG = 'DLSODE- Above warning has been issued I1 times. ' |
---|
| 1721 | CALL XERRWD (MSG, 50, 102, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1722 | MSG = ' It will not be issued again for this problem' |
---|
| 1723 | CALL XERRWD (MSG, 50, 102, 0, 1, MXHNIL, 0, 0, 0.0D0, 0.0D0) |
---|
| 1724 | 290 CONTINUE |
---|
| 1725 | !----------------------------------------------------------------------- |
---|
| 1726 | ! CALL DSTODE(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,DPREPJ,DSOLSY) |
---|
| 1727 | !----------------------------------------------------------------------- |
---|
| 1728 | CALL DSTODE (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), RWORK(LEWT), & |
---|
| 1729 | RWORK(LSAVF), RWORK(LACOR), RWORK(LWM), IWORK(LIWM), & |
---|
| 1730 | F, JAC) |
---|
| 1731 | !F, JAC, DPREPJ, DSOLSY) |
---|
| 1732 | KGO = 1 - KFLAG |
---|
| 1733 | GO TO (300, 530, 540), KGO |
---|
| 1734 | !----------------------------------------------------------------------- |
---|
| 1735 | ! Block F. |
---|
| 1736 | ! The following block handles the case of a successful return from the |
---|
| 1737 | ! core integrator (KFLAG = 0). Test for stop conditions. |
---|
| 1738 | !----------------------------------------------------------------------- |
---|
| 1739 | 300 INIT = 1 |
---|
| 1740 | GO TO (310, 400, 330, 340, 350), ITASK |
---|
| 1741 | ! ITASK = 1. If TOUT has been reached, interpolate. ------------------- |
---|
| 1742 | 310 IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250 |
---|
| 1743 | CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) |
---|
| 1744 | T = TOUT |
---|
| 1745 | GO TO 420 |
---|
| 1746 | ! ITASK = 3. Jump to exit if TOUT was reached. ------------------------ |
---|
| 1747 | 330 IF ((TN - TOUT)*H .GE. 0.0D0) GO TO 400 |
---|
| 1748 | GO TO 250 |
---|
| 1749 | ! ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary. |
---|
| 1750 | 340 IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 345 |
---|
| 1751 | CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG) |
---|
| 1752 | T = TOUT |
---|
| 1753 | GO TO 420 |
---|
| 1754 | 345 HMX = ABS(TN) + ABS(H) |
---|
| 1755 | IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX |
---|
| 1756 | IF (IHIT) GO TO 400 |
---|
| 1757 | TNEXT = TN + H*(1.0D0 + 4.0D0*UROUND) |
---|
| 1758 | IF ((TNEXT - TCRIT)*H .LE. 0.0D0) GO TO 250 |
---|
| 1759 | H = (TCRIT - TN)*(1.0D0 - 4.0D0*UROUND) |
---|
| 1760 | JSTART = -2 |
---|
| 1761 | GO TO 250 |
---|
| 1762 | ! ITASK = 5. See if TCRIT was reached and jump to exit. --------------- |
---|
| 1763 | 350 HMX = ABS(TN) + ABS(H) |
---|
| 1764 | IHIT = ABS(TN - TCRIT) .LE. 100.0D0*UROUND*HMX |
---|
| 1765 | !----------------------------------------------------------------------- |
---|
| 1766 | ! Block G. |
---|
| 1767 | ! The following block handles all successful returns from DLSODE. |
---|
| 1768 | ! If ITASK .NE. 1, Y is loaded from YH and T is set accordingly. |
---|
| 1769 | ! ISTATE is set to 2, and the optional outputs are loaded into the |
---|
| 1770 | ! work arrays before returning. |
---|
| 1771 | !----------------------------------------------------------------------- |
---|
| 1772 | 400 DO 410 I = 1,N |
---|
| 1773 | 410 Y(I) = RWORK(I+LYH-1) |
---|
| 1774 | T = TN |
---|
| 1775 | IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 420 |
---|
| 1776 | IF (IHIT) T = TCRIT |
---|
| 1777 | 420 ISTATE = 2 |
---|
| 1778 | RWORK(11) = HU |
---|
| 1779 | RWORK(12) = H |
---|
| 1780 | RWORK(13) = TN |
---|
| 1781 | IWORK(11) = NST |
---|
| 1782 | IWORK(12) = NFE |
---|
| 1783 | IWORK(13) = NJE |
---|
| 1784 | IWORK(14) = NQU |
---|
| 1785 | IWORK(15) = NQ |
---|
| 1786 | RETURN |
---|
| 1787 | !----------------------------------------------------------------------- |
---|
| 1788 | ! Block H. |
---|
| 1789 | ! The following block handles all unsuccessful returns other than |
---|
| 1790 | ! those for illegal input. First the error message routine is called. |
---|
| 1791 | ! If there was an error test or convergence test failure, IMXER is set. |
---|
| 1792 | ! Then Y is loaded from YH and T is set to TN. The optional outputs |
---|
| 1793 | ! are loaded into the work arrays before returning. |
---|
| 1794 | !----------------------------------------------------------------------- |
---|
| 1795 | ! The maximum number of steps was taken before reaching TOUT. ---------- |
---|
| 1796 | 500 MSG = 'DLSODE- At current T (=R1), MXSTEP (=I1) steps ' |
---|
| 1797 | CALL XERRWD (MSG, 50, 201, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1798 | MSG = ' taken on this call before reaching TOUT ' |
---|
| 1799 | CALL XERRWD (MSG, 50, 201, 0, 1, MXSTEP, 0, 1, TN, 0.0D0) |
---|
| 1800 | ISTATE = -1 |
---|
| 1801 | GO TO 580 |
---|
| 1802 | ! EWT(I) .LE. 0.0 for some I (not at start of problem). ---------------- |
---|
| 1803 | 510 EWTI = RWORK(LEWT+I-1) |
---|
| 1804 | MSG = 'DLSODE- At T (=R1), EWT(I1) has become R2 .LE. 0.' |
---|
| 1805 | CALL XERRWD (MSG, 50, 202, 0, 1, I, 0, 2, TN, EWTI) |
---|
| 1806 | ISTATE = -6 |
---|
| 1807 | GO TO 580 |
---|
| 1808 | ! Too much accuracy requested for machine precision. ------------------- |
---|
| 1809 | 520 MSG = 'DLSODE- At T (=R1), too much accuracy requested ' |
---|
| 1810 | CALL XERRWD (MSG, 50, 203, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1811 | MSG = ' for precision of machine.. see TOLSF (=R2) ' |
---|
| 1812 | CALL XERRWD (MSG, 50, 203, 0, 0, 0, 0, 2, TN, TOLSF) |
---|
| 1813 | RWORK(14) = TOLSF |
---|
| 1814 | ISTATE = -2 |
---|
| 1815 | GO TO 580 |
---|
| 1816 | ! KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. ----- |
---|
| 1817 | 530 MSG = 'DLSODE- At T(=R1) and step size H(=R2), the error' |
---|
| 1818 | CALL XERRWD (MSG, 50, 204, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1819 | MSG = ' test failed repeatedly or with ABS(H) = HMIN' |
---|
| 1820 | CALL XERRWD (MSG, 50, 204, 0, 0, 0, 0, 2, TN, H) |
---|
| 1821 | ISTATE = -4 |
---|
| 1822 | GO TO 560 |
---|
| 1823 | ! KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ---- |
---|
| 1824 | 540 MSG = 'DLSODE- At T (=R1) and step size H (=R2), the ' |
---|
| 1825 | CALL XERRWD (MSG, 50, 205, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1826 | MSG = ' corrector convergence failed repeatedly ' |
---|
| 1827 | CALL XERRWD (MSG, 50, 205, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1828 | MSG = ' or with ABS(H) = HMIN ' |
---|
| 1829 | CALL XERRWD (MSG, 30, 205, 0, 0, 0, 0, 2, TN, H) |
---|
| 1830 | ISTATE = -5 |
---|
| 1831 | ! Compute IMXER if relevant. ------------------------------------------- |
---|
| 1832 | 560 BIG = 0.0D0 |
---|
| 1833 | IMXER = 1 |
---|
| 1834 | DO 570 I = 1,N |
---|
| 1835 | SIZE = ABS(RWORK(I+LACOR-1)*RWORK(I+LEWT-1)) |
---|
| 1836 | IF (BIG .GE. SIZE) GO TO 570 |
---|
| 1837 | BIG = SIZE |
---|
| 1838 | IMXER = I |
---|
| 1839 | 570 CONTINUE |
---|
| 1840 | IWORK(16) = IMXER |
---|
| 1841 | ! Set Y vector, T, and optional outputs. ------------------------------- |
---|
| 1842 | 580 DO 590 I = 1,N |
---|
| 1843 | 590 Y(I) = RWORK(I+LYH-1) |
---|
| 1844 | T = TN |
---|
| 1845 | RWORK(11) = HU |
---|
| 1846 | RWORK(12) = H |
---|
| 1847 | RWORK(13) = TN |
---|
| 1848 | IWORK(11) = NST |
---|
| 1849 | IWORK(12) = NFE |
---|
| 1850 | IWORK(13) = NJE |
---|
| 1851 | IWORK(14) = NQU |
---|
| 1852 | IWORK(15) = NQ |
---|
| 1853 | RETURN |
---|
| 1854 | !----------------------------------------------------------------------- |
---|
| 1855 | ! Block I. |
---|
| 1856 | ! The following block handles all error returns due to illegal input |
---|
| 1857 | ! (ISTATE = -3), as detected before calling the core integrator. |
---|
| 1858 | ! First the error message routine is called. If the illegal input |
---|
| 1859 | ! is a negative ISTATE, the run is aborted (apparent infinite loop). |
---|
| 1860 | !----------------------------------------------------------------------- |
---|
| 1861 | 601 MSG = 'DLSODE- ISTATE (=I1) illegal ' |
---|
| 1862 | CALL XERRWD (MSG, 30, 1, 0, 1, ISTATE, 0, 0, 0.0D0, 0.0D0) |
---|
| 1863 | IF (ISTATE .LT. 0) GO TO 800 |
---|
| 1864 | GO TO 700 |
---|
| 1865 | 602 MSG = 'DLSODE- ITASK (=I1) illegal ' |
---|
| 1866 | CALL XERRWD (MSG, 30, 2, 0, 1, ITASK, 0, 0, 0.0D0, 0.0D0) |
---|
| 1867 | GO TO 700 |
---|
| 1868 | 603 MSG = 'DLSODE- ISTATE .GT. 1 but DLSODE not initialized ' |
---|
| 1869 | CALL XERRWD (MSG, 50, 3, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1870 | GO TO 700 |
---|
| 1871 | 604 MSG = 'DLSODE- NEQ (=I1) .LT. 1 ' |
---|
| 1872 | CALL XERRWD (MSG, 30, 4, 0, 1, NEQ, 0, 0, 0.0D0, 0.0D0) |
---|
| 1873 | GO TO 700 |
---|
| 1874 | 605 MSG = 'DLSODE- ISTATE = 3 and NEQ increased (I1 to I2) ' |
---|
| 1875 | CALL XERRWD (MSG, 50, 5, 0, 2, N, NEQ, 0, 0.0D0, 0.0D0) |
---|
| 1876 | GO TO 700 |
---|
| 1877 | 606 MSG = 'DLSODE- ITOL (=I1) illegal ' |
---|
| 1878 | CALL XERRWD (MSG, 30, 6, 0, 1, ITOL, 0, 0, 0.0D0, 0.0D0) |
---|
| 1879 | GO TO 700 |
---|
| 1880 | 607 MSG = 'DLSODE- IOPT (=I1) illegal ' |
---|
| 1881 | CALL XERRWD (MSG, 30, 7, 0, 1, IOPT, 0, 0, 0.0D0, 0.0D0) |
---|
| 1882 | GO TO 700 |
---|
| 1883 | 608 MSG = 'DLSODE- MF (=I1) illegal ' |
---|
| 1884 | CALL XERRWD (MSG, 30, 8, 0, 1, MF, 0, 0, 0.0D0, 0.0D0) |
---|
| 1885 | GO TO 700 |
---|
| 1886 | 609 MSG = 'DLSODE- ML (=I1) illegal.. .LT.0 or .GE.NEQ (=I2)' |
---|
| 1887 | CALL XERRWD (MSG, 50, 9, 0, 2, ML, NEQ, 0, 0.0D0, 0.0D0) |
---|
| 1888 | GO TO 700 |
---|
| 1889 | 610 MSG = 'DLSODE- MU (=I1) illegal.. .LT.0 or .GE.NEQ (=I2)' |
---|
| 1890 | CALL XERRWD (MSG, 50, 10, 0, 2, MU, NEQ, 0, 0.0D0, 0.0D0) |
---|
| 1891 | GO TO 700 |
---|
| 1892 | 611 MSG = 'DLSODE- MAXORD (=I1) .LT. 0 ' |
---|
| 1893 | CALL XERRWD (MSG, 30, 11, 0, 1, MAXORD, 0, 0, 0.0D0, 0.0D0) |
---|
| 1894 | GO TO 700 |
---|
| 1895 | 612 MSG = 'DLSODE- MXSTEP (=I1) .LT. 0 ' |
---|
| 1896 | CALL XERRWD (MSG, 30, 12, 0, 1, MXSTEP, 0, 0, 0.0D0, 0.0D0) |
---|
| 1897 | GO TO 700 |
---|
| 1898 | 613 MSG = 'DLSODE- MXHNIL (=I1) .LT. 0 ' |
---|
| 1899 | CALL XERRWD (MSG, 30, 13, 0, 1, MXHNIL, 0, 0, 0.0D0, 0.0D0) |
---|
| 1900 | GO TO 700 |
---|
| 1901 | 614 MSG = 'DLSODE- TOUT (=R1) behind T (=R2) ' |
---|
| 1902 | CALL XERRWD (MSG, 40, 14, 0, 0, 0, 0, 2, TOUT, T) |
---|
| 1903 | MSG = ' Integration direction is given by H0 (=R1) ' |
---|
| 1904 | CALL XERRWD (MSG, 50, 14, 0, 0, 0, 0, 1, H0, 0.0D0) |
---|
| 1905 | GO TO 700 |
---|
| 1906 | 615 MSG = 'DLSODE- HMAX (=R1) .LT. 0.0 ' |
---|
| 1907 | CALL XERRWD (MSG, 30, 15, 0, 0, 0, 0, 1, HMAX, 0.0D0) |
---|
| 1908 | GO TO 700 |
---|
| 1909 | 616 MSG = 'DLSODE- HMIN (=R1) .LT. 0.0 ' |
---|
| 1910 | CALL XERRWD (MSG, 30, 16, 0, 0, 0, 0, 1, HMIN, 0.0D0) |
---|
| 1911 | GO TO 700 |
---|
| 1912 | 617 CONTINUE |
---|
| 1913 | MSG='DLSODE- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)' |
---|
| 1914 | CALL XERRWD (MSG, 60, 17, 0, 2, LENRW, LRW, 0, 0.0D0, 0.0D0) |
---|
| 1915 | GO TO 700 |
---|
| 1916 | 618 CONTINUE |
---|
| 1917 | MSG='DLSODE- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)' |
---|
| 1918 | CALL XERRWD (MSG, 60, 18, 0, 2, LENIW, LIW, 0, 0.0D0, 0.0D0) |
---|
| 1919 | GO TO 700 |
---|
| 1920 | 619 MSG = 'DLSODE- RelTol(I1) is R1 .LT. 0.0 ' |
---|
| 1921 | CALL XERRWD (MSG, 40, 19, 0, 1, I, 0, 1, RelTolI, 0.0D0) |
---|
| 1922 | GO TO 700 |
---|
| 1923 | 620 MSG = 'DLSODE- AbsTol(I1) is R1 .LT. 0.0 ' |
---|
| 1924 | CALL XERRWD (MSG, 40, 20, 0, 1, I, 0, 1, AbsTolI, 0.0D0) |
---|
| 1925 | GO TO 700 |
---|
| 1926 | 621 EWTI = RWORK(LEWT+I-1) |
---|
| 1927 | MSG = 'DLSODE- EWT(I1) is R1 .LE. 0.0 ' |
---|
| 1928 | CALL XERRWD (MSG, 40, 21, 0, 1, I, 0, 1, EWTI, 0.0D0) |
---|
| 1929 | GO TO 700 |
---|
| 1930 | 622 CONTINUE |
---|
| 1931 | MSG='DLSODE- TOUT (=R1) too close to T(=R2) to start integration' |
---|
| 1932 | CALL XERRWD (MSG, 60, 22, 0, 0, 0, 0, 2, TOUT, T) |
---|
| 1933 | GO TO 700 |
---|
| 1934 | 623 CONTINUE |
---|
| 1935 | MSG='DLSODE- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) ' |
---|
| 1936 | CALL XERRWD (MSG, 60, 23, 0, 1, ITASK, 0, 2, TOUT, TP) |
---|
| 1937 | GO TO 700 |
---|
| 1938 | 624 CONTINUE |
---|
| 1939 | MSG='DLSODE- ITASK = 4 OR 5 and TCRIT (=R1) behind TCUR (=R2) ' |
---|
| 1940 | CALL XERRWD (MSG, 60, 24, 0, 0, 0, 0, 2, TCRIT, TN) |
---|
| 1941 | GO TO 700 |
---|
| 1942 | 625 CONTINUE |
---|
| 1943 | MSG='DLSODE- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) ' |
---|
| 1944 | CALL XERRWD (MSG, 60, 25, 0, 0, 0, 0, 2, TCRIT, TOUT) |
---|
| 1945 | GO TO 700 |
---|
| 1946 | 626 MSG = 'DLSODE- At start of problem, too much accuracy ' |
---|
| 1947 | CALL XERRWD (MSG, 50, 26, 0, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1948 | MSG=' requested for precision of machine.. See TOLSF (=R1) ' |
---|
| 1949 | CALL XERRWD (MSG, 60, 26, 0, 0, 0, 0, 1, TOLSF, 0.0D0) |
---|
| 1950 | RWORK(14) = TOLSF |
---|
| 1951 | GO TO 700 |
---|
| 1952 | 627 MSG = 'DLSODE- Trouble in DINTDY. ITASK = I1, TOUT = R1' |
---|
| 1953 | CALL XERRWD (MSG, 50, 27, 0, 1, ITASK, 0, 1, TOUT, 0.0D0) |
---|
| 1954 | ! |
---|
| 1955 | 700 ISTATE = -3 |
---|
| 1956 | RETURN |
---|
| 1957 | ! |
---|
| 1958 | 800 MSG = 'DLSODE- Run aborted.. apparent infinite loop ' |
---|
| 1959 | CALL XERRWD (MSG, 50, 303, 2, 0, 0, 0, 0, 0.0D0, 0.0D0) |
---|
| 1960 | RETURN |
---|
| 1961 | !----------------------- END OF SUBROUTINE DLSODE ---------------------- |
---|
| 1962 | !END SUBROUTINE DLSODE |
---|
| 1963 | CONTAINS |
---|
| 1964 | |
---|
| 1965 | |
---|
| 1966 | !DECK DUMACH |
---|
| 1967 | KPP_REAL FUNCTION DUMACH () |
---|
| 1968 | !***BEGIN PROLOGUE DUMACH |
---|
| 1969 | !***PURPOSE Compute the unit roundoff of the machine. |
---|
| 1970 | !***CATEGORY R1 |
---|
| 1971 | !***TYPE KPP_REAL (RUMACH-S, DUMACH-D) |
---|
| 1972 | !***KEYWORDS MACHINE CONSTANTS |
---|
| 1973 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 1974 | !***DESCRIPTION |
---|
| 1975 | ! *Usage: |
---|
| 1976 | ! KPP_REAL A, DUMACH |
---|
| 1977 | ! A = DUMACH() |
---|
| 1978 | ! |
---|
| 1979 | ! *Function Return Values: |
---|
| 1980 | ! A : the unit roundoff of the machine. |
---|
| 1981 | ! |
---|
| 1982 | ! *Description: |
---|
| 1983 | ! The unit roundoff is defined as the smallest positive machine |
---|
| 1984 | ! number u such that 1.0 + u .ne. 1.0. This is computed by DUMACH |
---|
| 1985 | ! in a machine-independent manner. |
---|
| 1986 | ! |
---|
| 1987 | !***REFERENCES (NONE) |
---|
| 1988 | !***ROUTINES CALLED DUMSUM |
---|
| 1989 | !***REVISION HISTORY (YYYYMMDD) |
---|
| 1990 | ! 19930216 DATE WRITTEN |
---|
| 1991 | ! 19930818 Added SLATEC-format prologue. (FNF) |
---|
| 1992 | ! 20030707 Added DUMSUM to force normal storage of COMP. (ACH) |
---|
| 1993 | !***END PROLOGUE DUMACH |
---|
| 1994 | ! |
---|
| 1995 | KPP_REAL U, COMP |
---|
| 1996 | !***FIRST EXECUTABLE STATEMENT DUMACH |
---|
| 1997 | U = 1.0D0 |
---|
| 1998 | 10 U = U*0.5D0 |
---|
| 1999 | CALL DUMSUM(1.0D0, U, COMP) |
---|
| 2000 | IF (COMP .NE. 1.0D0) GO TO 10 |
---|
| 2001 | DUMACH = U*2.0D0 |
---|
| 2002 | RETURN |
---|
| 2003 | !----------------------- End of Function DUMACH ------------------------ |
---|
| 2004 | END FUNCTION DUMACH |
---|
| 2005 | |
---|
| 2006 | SUBROUTINE DUMSUM(A,B,C) |
---|
| 2007 | ! Routine to force normal storing of A + B, for DUMACH. |
---|
| 2008 | KPP_REAL A, B, C |
---|
| 2009 | C = A + B |
---|
| 2010 | RETURN |
---|
| 2011 | END SUBROUTINE DUMSUM |
---|
| 2012 | !DECK DCFODE |
---|
| 2013 | SUBROUTINE DCFODE (METH, ELCO, TESCO) |
---|
| 2014 | !***BEGIN PROLOGUE DCFODE |
---|
| 2015 | !***SUBSIDIARY |
---|
| 2016 | !***PURPOSE Set ODE integrator coefficients. |
---|
| 2017 | !***TYPE KPP_REAL (SCFODE-S, DCFODE-D) |
---|
| 2018 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 2019 | !***DESCRIPTION |
---|
| 2020 | ! |
---|
| 2021 | ! DCFODE is called by the integrator routine to set coefficients |
---|
| 2022 | ! needed there. The coefficients for the current method, as |
---|
| 2023 | ! given by the value of METH, are set for all orders and saved. |
---|
| 2024 | ! The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2. |
---|
| 2025 | ! (A smaller value of the maximum order is also allowed.) |
---|
| 2026 | ! DCFODE is called once at the beginning of the problem, |
---|
| 2027 | ! and is not called again unless and until METH is changed. |
---|
| 2028 | ! |
---|
| 2029 | ! The ELCO array contains the basic method coefficients. |
---|
| 2030 | ! The coefficients el(i), 1 .le. i .le. nq+1, for the method of |
---|
| 2031 | ! order nq are stored in ELCO(i,nq). They are given by a genetrating |
---|
| 2032 | ! polynomial, i.e., |
---|
| 2033 | ! l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq. |
---|
| 2034 | ! For the implicit Adams methods, l(x) is given by |
---|
| 2035 | ! dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0. |
---|
| 2036 | ! For the BDF methods, l(x) is given by |
---|
| 2037 | ! l(x) = (x+1)*(x+2)* ... *(x+nq)/K, |
---|
| 2038 | ! where K = factorial(nq)*(1 + 1/2 + ... + 1/nq). |
---|
| 2039 | ! |
---|
| 2040 | ! The TESCO array contains test constants used for the |
---|
| 2041 | ! local error test and the selection of step size and/or order. |
---|
| 2042 | ! At order nq, TESCO(k,nq) is used for the selection of step |
---|
| 2043 | ! size at order nq - 1 if k = 1, at order nq if k = 2, and at order |
---|
| 2044 | ! nq + 1 if k = 3. |
---|
| 2045 | ! |
---|
| 2046 | !***SEE ALSO DLSODE |
---|
| 2047 | !***ROUTINES CALLED (NONE) |
---|
| 2048 | !***REVISION HISTORY (YYMMDD) |
---|
| 2049 | ! 791129 DATE WRITTEN |
---|
| 2050 | ! 890501 Modified prologue to SLATEC/LDOC format. (FNF) |
---|
| 2051 | ! 890503 Minor cosmetic changes. (FNF) |
---|
| 2052 | ! 930809 Renamed to allow single/double precision versions. (ACH) |
---|
| 2053 | !***END PROLOGUE DCFODE |
---|
| 2054 | !**End |
---|
| 2055 | INTEGER METH |
---|
| 2056 | INTEGER I, IB, NQ, NQM1, NQP1 |
---|
| 2057 | KPP_REAL ELCO(13,12), TESCO(3,12), PC(12) |
---|
| 2058 | KPP_REAL AGAMQ, FNQ, FNQM1, PINT, RAGQ, RQFAC, RQ1FAC, TSIGN, XPIN |
---|
| 2059 | ! |
---|
| 2060 | !***FIRST EXECUTABLE STATEMENT DCFODE |
---|
| 2061 | GO TO (100, 200), METH |
---|
| 2062 | ! |
---|
| 2063 | 100 ELCO(1,1) = 1.0D0 |
---|
| 2064 | ELCO(2,1) = 1.0D0 |
---|
| 2065 | TESCO(1,1) = 0.0D0 |
---|
| 2066 | TESCO(2,1) = 2.0D0 |
---|
| 2067 | TESCO(1,2) = 1.0D0 |
---|
| 2068 | TESCO(3,12) = 0.0D0 |
---|
| 2069 | PC(1) = 1.0D0 |
---|
| 2070 | RQFAC = 1.0D0 |
---|
| 2071 | DO 140 NQ = 2,12 |
---|
| 2072 | !----------------------------------------------------------------------- |
---|
| 2073 | ! The PC array will contain the coefficients of the polynomial |
---|
| 2074 | ! p(x) = (x+1)*(x+2)*...*(x+nq-1). |
---|
| 2075 | ! Initially, p(x) = 1. |
---|
| 2076 | !----------------------------------------------------------------------- |
---|
| 2077 | RQ1FAC = RQFAC |
---|
| 2078 | RQFAC = RQFAC/NQ |
---|
| 2079 | NQM1 = NQ - 1 |
---|
| 2080 | FNQM1 = NQM1 |
---|
| 2081 | NQP1 = NQ + 1 |
---|
| 2082 | ! Form coefficients of p(x)*(x+nq-1). ---------------------------------- |
---|
| 2083 | PC(NQ) = 0.0D0 |
---|
| 2084 | DO 110 IB = 1,NQM1 |
---|
| 2085 | I = NQP1 - IB |
---|
| 2086 | 110 PC(I) = PC(I-1) + FNQM1*PC(I) |
---|
| 2087 | PC(1) = FNQM1*PC(1) |
---|
| 2088 | ! Compute integral, -1 to 0, of p(x) and x*p(x). ----------------------- |
---|
| 2089 | PINT = PC(1) |
---|
| 2090 | XPIN = PC(1)/2.0D0 |
---|
| 2091 | TSIGN = 1.0D0 |
---|
| 2092 | DO 120 I = 2,NQ |
---|
| 2093 | TSIGN = -TSIGN |
---|
| 2094 | PINT = PINT + TSIGN*PC(I)/I |
---|
| 2095 | 120 XPIN = XPIN + TSIGN*PC(I)/(I+1) |
---|
| 2096 | ! Store coefficients in ELCO and TESCO. -------------------------------- |
---|
| 2097 | ELCO(1,NQ) = PINT*RQ1FAC |
---|
| 2098 | ELCO(2,NQ) = 1.0D0 |
---|
| 2099 | DO 130 I = 2,NQ |
---|
| 2100 | 130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/I |
---|
| 2101 | AGAMQ = RQFAC*XPIN |
---|
| 2102 | RAGQ = 1.0D0/AGAMQ |
---|
| 2103 | TESCO(2,NQ) = RAGQ |
---|
| 2104 | IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/NQP1 |
---|
| 2105 | TESCO(3,NQM1) = RAGQ |
---|
| 2106 | 140 CONTINUE |
---|
| 2107 | RETURN |
---|
| 2108 | ! |
---|
| 2109 | 200 PC(1) = 1.0D0 |
---|
| 2110 | RQ1FAC = 1.0D0 |
---|
| 2111 | DO 230 NQ = 1,5 |
---|
| 2112 | !----------------------------------------------------------------------- |
---|
| 2113 | ! The PC array will contain the coefficients of the polynomial |
---|
| 2114 | ! p(x) = (x+1)*(x+2)*...*(x+nq). |
---|
| 2115 | ! Initially, p(x) = 1. |
---|
| 2116 | !----------------------------------------------------------------------- |
---|
| 2117 | FNQ = NQ |
---|
| 2118 | NQP1 = NQ + 1 |
---|
| 2119 | ! Form coefficients of p(x)*(x+nq). ------------------------------------ |
---|
| 2120 | PC(NQP1) = 0.0D0 |
---|
| 2121 | DO 210 IB = 1,NQ |
---|
| 2122 | I = NQ + 2 - IB |
---|
| 2123 | 210 PC(I) = PC(I-1) + FNQ*PC(I) |
---|
| 2124 | PC(1) = FNQ*PC(1) |
---|
| 2125 | ! Store coefficients in ELCO and TESCO. -------------------------------- |
---|
| 2126 | DO 220 I = 1,NQP1 |
---|
| 2127 | 220 ELCO(I,NQ) = PC(I)/PC(2) |
---|
| 2128 | ELCO(2,NQ) = 1.0D0 |
---|
| 2129 | TESCO(1,NQ) = RQ1FAC |
---|
| 2130 | TESCO(2,NQ) = NQP1/ELCO(1,NQ) |
---|
| 2131 | TESCO(3,NQ) = (NQ+2)/ELCO(1,NQ) |
---|
| 2132 | RQ1FAC = RQ1FAC/FNQ |
---|
| 2133 | 230 CONTINUE |
---|
| 2134 | RETURN |
---|
| 2135 | !----------------------- END OF SUBROUTINE DCFODE ---------------------- |
---|
| 2136 | END SUBROUTINE DCFODE |
---|
| 2137 | !DECK DINTDY |
---|
| 2138 | SUBROUTINE DINTDY (T, K, YH, NYH, DKY, IFLAG) |
---|
| 2139 | !***BEGIN PROLOGUE DINTDY |
---|
| 2140 | !***SUBSIDIARY |
---|
| 2141 | !***PURPOSE Interpolate solution derivatives. |
---|
| 2142 | !***TYPE KPP_REAL (SINTDY-S, DINTDY-D) |
---|
| 2143 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 2144 | !***DESCRIPTION |
---|
| 2145 | ! |
---|
| 2146 | ! DINTDY computes interpolated values of the K-th derivative of the |
---|
| 2147 | ! dependent variable vector y, and stores it in DKY. This routine |
---|
| 2148 | ! is called within the package with K = 0 and T = TOUT, but may |
---|
| 2149 | ! also be called by the user for any K up to the current order. |
---|
| 2150 | ! (See detailed instructions in the usage documentation.) |
---|
| 2151 | ! |
---|
| 2152 | ! The computed values in DKY are gotten by interpolation using the |
---|
| 2153 | ! Nordsieck history array YH. This array corresponds uniquely to a |
---|
| 2154 | ! vector-valued polynomial of degree NQCUR or less, and DKY is set |
---|
| 2155 | ! to the K-th derivative of this polynomial at T. |
---|
| 2156 | ! The formula for DKY is: |
---|
| 2157 | ! q |
---|
| 2158 | ! DKY(i) = sum c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1) |
---|
| 2159 | ! j=K |
---|
| 2160 | ! where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR. |
---|
| 2161 | ! The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are |
---|
| 2162 | ! communicated by COMMON. The above sum is done in reverse order. |
---|
| 2163 | ! IFLAG is returned negative if either K or T is out of bounds. |
---|
| 2164 | ! |
---|
| 2165 | !***SEE ALSO DLSODE |
---|
| 2166 | !***ROUTINES CALLED XERRWD |
---|
| 2167 | !***COMMON BLOCKS DLS001 |
---|
| 2168 | !***REVISION HISTORY (YYMMDD) |
---|
| 2169 | ! 791129 DATE WRITTEN |
---|
| 2170 | ! 890501 Modified prologue to SLATEC/LDOC format. (FNF) |
---|
| 2171 | ! 890503 Minor cosmetic changes. (FNF) |
---|
| 2172 | ! 930809 Renamed to allow single/double precision versions. (ACH) |
---|
| 2173 | ! 010418 Reduced size of Common block /DLS001/. (ACH) |
---|
| 2174 | ! 031105 Restored 'own' variables to Common block /DLS001/, to |
---|
| 2175 | ! enable interrupt/restart feature. (ACH) |
---|
| 2176 | ! 050427 Corrected roundoff decrement in TP. (ACH) |
---|
| 2177 | !***END PROLOGUE DINTDY |
---|
| 2178 | !**End |
---|
| 2179 | INTEGER K, NYH, IFLAG |
---|
| 2180 | KPP_REAL T, YH(NYH,*), DKY(*) |
---|
| 2181 | INTEGER IOWND, IOWNS, & |
---|
| 2182 | ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, & |
---|
| 2183 | LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, & |
---|
| 2184 | MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
---|
| 2185 | KPP_REAL ROWNS, & |
---|
| 2186 | CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND |
---|
| 2187 | COMMON /DLS001/ ROWNS(209), & |
---|
| 2188 | CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, & |
---|
| 2189 | IOWND(6), IOWNS(6), & |
---|
| 2190 | ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, & |
---|
| 2191 | LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, & |
---|
| 2192 | MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
---|
| 2193 | INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1 |
---|
| 2194 | KPP_REAL C, R, S, TP |
---|
| 2195 | CHARACTER*80 MSG |
---|
| 2196 | ! |
---|
| 2197 | !***FIRST EXECUTABLE STATEMENT DINTDY |
---|
| 2198 | IFLAG = 0 |
---|
| 2199 | IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80 |
---|
| 2200 | TP = TN - HU - 100.0D0*UROUND*SIGN(ABS(TN) + ABS(HU), HU) |
---|
| 2201 | IF ((T-TP)*(T-TN) .GT. 0.0D0) GO TO 90 |
---|
| 2202 | ! |
---|
| 2203 | S = (T - TN)/H |
---|
| 2204 | IC = 1 |
---|
| 2205 | IF (K .EQ. 0) GO TO 15 |
---|
| 2206 | JJ1 = L - K |
---|
| 2207 | DO 10 JJ = JJ1,NQ |
---|
| 2208 | 10 IC = IC*JJ |
---|
| 2209 | 15 C = IC |
---|
| 2210 | DO 20 I = 1,N |
---|
| 2211 | 20 DKY(I) = C*YH(I,L) |
---|
| 2212 | IF (K .EQ. NQ) GO TO 55 |
---|
| 2213 | JB2 = NQ - K |
---|
| 2214 | DO 50 JB = 1,JB2 |
---|
| 2215 | J = NQ - JB |
---|
| 2216 | JP1 = J + 1 |
---|
| 2217 | IC = 1 |
---|
| 2218 | IF (K .EQ. 0) GO TO 35 |
---|
| 2219 | JJ1 = JP1 - K |
---|
| 2220 | DO 30 JJ = JJ1,J |
---|
| 2221 | 30 IC = IC*JJ |
---|
| 2222 | 35 C = IC |
---|
| 2223 | DO 40 I = 1,N |
---|
| 2224 | 40 DKY(I) = C*YH(I,JP1) + S*DKY(I) |
---|
| 2225 | 50 CONTINUE |
---|
| 2226 | IF (K .EQ. 0) RETURN |
---|
| 2227 | 55 R = H**(-K) |
---|
| 2228 | DO 60 I = 1,N |
---|
| 2229 | 60 DKY(I) = R*DKY(I) |
---|
| 2230 | RETURN |
---|
| 2231 | ! |
---|
| 2232 | 80 MSG = 'DINTDY- K (=I1) illegal ' |
---|
| 2233 | CALL XERRWD (MSG, 30, 51, 0, 1, K, 0, 0, 0.0D0, 0.0D0) |
---|
| 2234 | IFLAG = -1 |
---|
| 2235 | RETURN |
---|
| 2236 | 90 MSG = 'DINTDY- T (=R1) illegal ' |
---|
| 2237 | CALL XERRWD (MSG, 30, 52, 0, 0, 0, 0, 1, T, 0.0D0) |
---|
| 2238 | MSG=' T not in interval TCUR - HU (= R1) to TCUR (=R2) ' |
---|
| 2239 | CALL XERRWD (MSG, 60, 52, 0, 0, 0, 0, 2, TP, TN) |
---|
| 2240 | IFLAG = -2 |
---|
| 2241 | RETURN |
---|
| 2242 | !----------------------- END OF SUBROUTINE DINTDY ---------------------- |
---|
| 2243 | END SUBROUTINE DINTDY |
---|
| 2244 | !DECK DPREPJ |
---|
| 2245 | SUBROUTINE DPREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC) |
---|
| 2246 | !***BEGIN PROLOGUE DPREPJ |
---|
| 2247 | !***SUBSIDIARY |
---|
| 2248 | !***PURPOSE Compute and process Newton iteration matrix. |
---|
| 2249 | !***TYPE KPP_REAL (SPREPJ-S, DPREPJ-D) |
---|
| 2250 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 2251 | !***DESCRIPTION |
---|
| 2252 | ! |
---|
| 2253 | ! DPREPJ is called by DSTODE to compute and process the matrix |
---|
| 2254 | ! P = I - h*el(1)*J , where J is an approximation to the Jacobian. |
---|
| 2255 | ! Here J is computed by the user-supplied routine JAC if |
---|
| 2256 | ! MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5. |
---|
| 2257 | ! If MITER = 3, a diagonal approximation to J is used. |
---|
| 2258 | ! J is stored in WM and replaced by P. If MITER .ne. 3, P is then |
---|
| 2259 | ! subjected to LU decomposition in preparation for later solution |
---|
| 2260 | ! of linear systems with P as coefficient matrix. This is done |
---|
| 2261 | ! by DGEFA if MITER = 1 or 2, and by DGBFA if MITER = 4 or 5. |
---|
| 2262 | ! |
---|
| 2263 | ! In addition to variables described in DSTODE and DLSODE prologues, |
---|
| 2264 | ! communication with DPREPJ uses the following: |
---|
| 2265 | ! Y = array containing predicted values on entry. |
---|
| 2266 | ! FTEM = work array of length N (ACOR in DSTODE). |
---|
| 2267 | ! SAVF = array containing f evaluated at predicted y. |
---|
| 2268 | ! WM = real work space for matrices. On output it contains the |
---|
| 2269 | ! inverse diagonal matrix if MITER = 3 and the LU decomposition |
---|
| 2270 | ! of P if MITER is 1, 2 , 4, or 5. |
---|
| 2271 | ! Storage of matrix elements starts at WM(3). |
---|
| 2272 | ! WM also contains the following matrix-related data: |
---|
| 2273 | ! WM(1) = SQRT(UROUND), used in numerical Jacobian increments. |
---|
| 2274 | ! WM(2) = H*EL0, saved for later use if MITER = 3. |
---|
| 2275 | ! IWM = integer work space containing pivot information, starting at |
---|
| 2276 | ! IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band |
---|
| 2277 | ! parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5. |
---|
| 2278 | ! EL0 = EL(1) (input). |
---|
| 2279 | ! IERPJ = output error flag, = 0 if no trouble, .gt. 0 if |
---|
| 2280 | ! P matrix found to be singular. |
---|
| 2281 | ! JCUR = output flag = 1 to indicate that the Jacobian matrix |
---|
| 2282 | ! (or approximation) is now current. |
---|
| 2283 | ! This routine also uses the COMMON variables EL0, H, TN, UROUND, |
---|
| 2284 | ! MITER, N, NFE, and NJE. |
---|
| 2285 | ! |
---|
| 2286 | !***SEE ALSO DLSODE |
---|
| 2287 | !***ROUTINES CALLED DGBFA, DGEFA, DVNORM |
---|
| 2288 | !***COMMON BLOCKS DLS001 |
---|
| 2289 | !***REVISION HISTORY (YYMMDD) |
---|
| 2290 | ! 791129 DATE WRITTEN |
---|
| 2291 | ! 890501 Modified prologue to SLATEC/LDOC format. (FNF) |
---|
| 2292 | ! 890504 Minor cosmetic changes. (FNF) |
---|
| 2293 | ! 930809 Renamed to allow single/double precision versions. (ACH) |
---|
| 2294 | ! 010418 Reduced size of Common block /DLS001/. (ACH) |
---|
| 2295 | ! 031105 Restored 'own' variables to Common block /DLS001/, to |
---|
| 2296 | ! enable interrupt/restart feature. (ACH) |
---|
| 2297 | !***END PROLOGUE DPREPJ |
---|
| 2298 | !**End |
---|
| 2299 | EXTERNAL F, JAC |
---|
| 2300 | INTEGER NEQ, NYH, IWM(*) |
---|
| 2301 | KPP_REAL Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*), WM(*) |
---|
| 2302 | INTEGER IOWND, IOWNS, & |
---|
| 2303 | ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, & |
---|
| 2304 | LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, & |
---|
| 2305 | MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
---|
| 2306 | KPP_REAL ROWNS, & |
---|
| 2307 | CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND |
---|
| 2308 | COMMON /DLS001/ ROWNS(209), & |
---|
| 2309 | CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, & |
---|
| 2310 | IOWND(6), IOWNS(6), & |
---|
| 2311 | ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, & |
---|
| 2312 | LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, & |
---|
| 2313 | MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
---|
| 2314 | INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP, & |
---|
| 2315 | MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1 |
---|
| 2316 | KPP_REAL CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ |
---|
| 2317 | !KPP_REAL DVNORM |
---|
| 2318 | ! |
---|
| 2319 | !***FIRST EXECUTABLE STATEMENT DPREPJ |
---|
| 2320 | NJE = NJE + 1 |
---|
| 2321 | IERPJ = 0 |
---|
| 2322 | JCUR = 1 |
---|
| 2323 | HL0 = H*EL0 |
---|
| 2324 | CON = -HL0 |
---|
| 2325 | |
---|
| 2326 | #ifdef FULL_ALGEBRA |
---|
| 2327 | LENP = N*N |
---|
| 2328 | DO i = 1,LENP |
---|
| 2329 | WM(i+2) = 0.0D0 |
---|
| 2330 | END DO |
---|
| 2331 | CALL JAC_CHEM (NEQ, TN, Y, WM(3)) |
---|
| 2332 | DO I = 1,LENP |
---|
| 2333 | WM(I+2) = WM(I+2)*CON |
---|
| 2334 | END DO |
---|
| 2335 | ! Add identity matrix |
---|
| 2336 | J = 3 |
---|
| 2337 | NP1 = N + 1 |
---|
| 2338 | DO I = 1,N |
---|
| 2339 | WM(J) = WM(J) + 1.0D0 |
---|
| 2340 | J = J + NP1 |
---|
| 2341 | END DO |
---|
| 2342 | ! Do LU decomposition on P |
---|
| 2343 | CALL DGETRF(N,N,WM(3),N,IWM(21),IER) |
---|
| 2344 | #else |
---|
| 2345 | CALL JAC_CHEM (NEQ, TN, Y, WM(3)) |
---|
| 2346 | DO i = 1,LU_NONZERO |
---|
| 2347 | WM(i+2) = WM(i+2)*CON |
---|
| 2348 | END DO |
---|
| 2349 | ! Add identity matrix |
---|
| 2350 | DO i = 1,N |
---|
| 2351 | j = 2+LU_DIAG(i) |
---|
| 2352 | WM(j) = WM(j) + 1.0D0 |
---|
| 2353 | END DO |
---|
| 2354 | ! Do LU decomposition on P |
---|
| 2355 | CALL KppDecomp(WM(3),IER) |
---|
| 2356 | #endif |
---|
| 2357 | IF (IER .NE. 0) IERPJ = 1 |
---|
| 2358 | RETURN |
---|
| 2359 | !----------------------- END OF SUBROUTINE DPREPJ ---------------------- |
---|
| 2360 | END SUBROUTINE DPREPJ |
---|
| 2361 | !DECK DSOLSY |
---|
| 2362 | SUBROUTINE DSOLSY (WM, IWM, X, TEM) |
---|
| 2363 | !***BEGIN PROLOGUE DSOLSY |
---|
| 2364 | !***SUBSIDIARY |
---|
| 2365 | !***PURPOSE ODEPACK linear system solver. |
---|
| 2366 | !***TYPE KPP_REAL (SSOLSY-S, DSOLSY-D) |
---|
| 2367 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 2368 | !***DESCRIPTION |
---|
| 2369 | ! |
---|
| 2370 | ! This routine manages the solution of the linear system arising from |
---|
| 2371 | ! a chord iteration. It is called if MITER .ne. 0. |
---|
| 2372 | ! If MITER is 1 or 2, it calls DGESL to accomplish this. |
---|
| 2373 | ! If MITER = 3 it updates the coefficient h*EL0 in the diagonal |
---|
| 2374 | ! matrix, and then computes the solution. |
---|
| 2375 | ! If MITER is 4 or 5, it calls DGBSL. |
---|
| 2376 | ! Communication with DSOLSY uses the following variables: |
---|
| 2377 | ! WM = real work space containing the inverse diagonal matrix if |
---|
| 2378 | ! MITER = 3 and the LU decomposition of the matrix otherwise. |
---|
| 2379 | ! Storage of matrix elements starts at WM(3). |
---|
| 2380 | ! WM also contains the following matrix-related data: |
---|
| 2381 | ! WM(1) = SQRT(UROUND) (not used here), |
---|
| 2382 | ! WM(2) = HL0, the previous value of h*EL0, used if MITER = 3. |
---|
| 2383 | ! IWM = integer work space containing pivot information, starting at |
---|
| 2384 | ! IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band |
---|
| 2385 | ! parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5. |
---|
| 2386 | ! X = the right-hand side vector on input, and the solution vector |
---|
| 2387 | ! on output, of length N. |
---|
| 2388 | ! TEM = vector of work space of length N, not used in this version. |
---|
| 2389 | ! IERSL = output flag (in COMMON). IERSL = 0 if no trouble occurred. |
---|
| 2390 | ! IERSL = 1 if a singular matrix arose with MITER = 3. |
---|
| 2391 | ! This routine also uses the COMMON variables EL0, H, MITER, and N. |
---|
| 2392 | ! |
---|
| 2393 | !***SEE ALSO DLSODE |
---|
| 2394 | !***ROUTINES CALLED DGBSL, DGESL |
---|
| 2395 | !***COMMON BLOCKS DLS001 |
---|
| 2396 | !***REVISION HISTORY (YYMMDD) |
---|
| 2397 | ! 791129 DATE WRITTEN |
---|
| 2398 | ! 890501 Modified prologue to SLATEC/LDOC format. (FNF) |
---|
| 2399 | ! 890503 Minor cosmetic changes. (FNF) |
---|
| 2400 | ! 930809 Renamed to allow single/double precision versions. (ACH) |
---|
| 2401 | ! 010418 Reduced size of Common block /DLS001/. (ACH) |
---|
| 2402 | ! 031105 Restored 'own' variables to Common block /DLS001/, to |
---|
| 2403 | ! enable interrupt/restart feature. (ACH) |
---|
| 2404 | !***END PROLOGUE DSOLSY |
---|
| 2405 | !**End |
---|
| 2406 | INTEGER IWM(*) |
---|
| 2407 | KPP_REAL WM(*), X(*), TEM(*) |
---|
| 2408 | INTEGER IOWND, IOWNS, & |
---|
| 2409 | ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, & |
---|
| 2410 | LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, & |
---|
| 2411 | MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
---|
| 2412 | KPP_REAL ROWNS, & |
---|
| 2413 | CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND |
---|
| 2414 | COMMON /DLS001/ ROWNS(209), & |
---|
| 2415 | CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, & |
---|
| 2416 | IOWND(6), IOWNS(6), & |
---|
| 2417 | ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, & |
---|
| 2418 | LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, & |
---|
| 2419 | MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
---|
| 2420 | INTEGER I, MEBAND, ML, MU |
---|
| 2421 | KPP_REAL DI, HL0, PHL0, R |
---|
| 2422 | ! |
---|
| 2423 | !***FIRST EXECUTABLE STATEMENT DSOLSY |
---|
| 2424 | IERSL = 0 |
---|
| 2425 | #ifdef FULL_ALGEBRA |
---|
| 2426 | CALL DGETRS ('N',N,1,WM(3),N,IWM(21),X,N,0) |
---|
| 2427 | #else |
---|
| 2428 | CALL KppSolve(WM(3),X) |
---|
| 2429 | #endif |
---|
| 2430 | RETURN |
---|
| 2431 | !----------------------- END OF SUBROUTINE DSOLSY ---------------------- |
---|
| 2432 | END SUBROUTINE DSOLSY |
---|
| 2433 | !DECK DSRCOM |
---|
| 2434 | SUBROUTINE DSRCOM (RSAV, ISAV, JOB) |
---|
| 2435 | !***BEGIN PROLOGUE DSRCOM |
---|
| 2436 | !***SUBSIDIARY |
---|
| 2437 | !***PURPOSE Save/restore ODEPACK COMMON blocks. |
---|
| 2438 | !***TYPE KPP_REAL (SSRCOM-S, DSRCOM-D) |
---|
| 2439 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 2440 | !***DESCRIPTION |
---|
| 2441 | ! |
---|
| 2442 | ! This routine saves or restores (depending on JOB) the contents of |
---|
| 2443 | ! the COMMON block DLS001, which is used internally |
---|
| 2444 | ! by one or more ODEPACK solvers. |
---|
| 2445 | ! |
---|
| 2446 | ! RSAV = real array of length 218 or more. |
---|
| 2447 | ! ISAV = integer array of length 37 or more. |
---|
| 2448 | ! JOB = flag indicating to save or restore the COMMON blocks: |
---|
| 2449 | ! JOB = 1 if COMMON is to be saved (written to RSAV/ISAV) |
---|
| 2450 | ! JOB = 2 if COMMON is to be restored (read from RSAV/ISAV) |
---|
| 2451 | ! A call with JOB = 2 presumes a prior call with JOB = 1. |
---|
| 2452 | ! |
---|
| 2453 | !***SEE ALSO DLSODE |
---|
| 2454 | !***ROUTINES CALLED (NONE) |
---|
| 2455 | !***COMMON BLOCKS DLS001 |
---|
| 2456 | !***REVISION HISTORY (YYMMDD) |
---|
| 2457 | ! 791129 DATE WRITTEN |
---|
| 2458 | ! 890501 Modified prologue to SLATEC/LDOC format. (FNF) |
---|
| 2459 | ! 890503 Minor cosmetic changes. (FNF) |
---|
| 2460 | ! 921116 Deleted treatment of block /EH0001/. (ACH) |
---|
| 2461 | ! 930801 Reduced Common block length by 2. (ACH) |
---|
| 2462 | ! 930809 Renamed to allow single/double precision versions. (ACH) |
---|
| 2463 | ! 010418 Reduced Common block length by 209+12. (ACH) |
---|
| 2464 | ! 031105 Restored 'own' variables to Common block /DLS001/, to |
---|
| 2465 | ! enable interrupt/restart feature. (ACH) |
---|
| 2466 | ! 031112 Added SAVE statement for data-loaded constants. |
---|
| 2467 | !***END PROLOGUE DSRCOM |
---|
| 2468 | !**End |
---|
| 2469 | INTEGER ISAV(*), JOB |
---|
| 2470 | INTEGER ILS |
---|
| 2471 | INTEGER I, LENILS, LENRLS |
---|
| 2472 | KPP_REAL RSAV(*), RLS |
---|
| 2473 | SAVE LENRLS, LENILS |
---|
| 2474 | COMMON /DLS001/ RLS(218), ILS(37) |
---|
| 2475 | DATA LENRLS/218/, LENILS/37/ |
---|
| 2476 | ! |
---|
| 2477 | !***FIRST EXECUTABLE STATEMENT DSRCOM |
---|
| 2478 | IF (JOB .EQ. 2) GO TO 100 |
---|
| 2479 | ! |
---|
| 2480 | DO 10 I = 1,LENRLS |
---|
| 2481 | 10 RSAV(I) = RLS(I) |
---|
| 2482 | DO 20 I = 1,LENILS |
---|
| 2483 | 20 ISAV(I) = ILS(I) |
---|
| 2484 | RETURN |
---|
| 2485 | ! |
---|
| 2486 | 100 CONTINUE |
---|
| 2487 | DO 110 I = 1,LENRLS |
---|
| 2488 | 110 RLS(I) = RSAV(I) |
---|
| 2489 | DO 120 I = 1,LENILS |
---|
| 2490 | 120 ILS(I) = ISAV(I) |
---|
| 2491 | RETURN |
---|
| 2492 | !----------------------- END OF SUBROUTINE DSRCOM ---------------------- |
---|
| 2493 | END SUBROUTINE DSRCOM |
---|
| 2494 | !DECK DSTODE |
---|
| 2495 | SUBROUTINE DSTODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, & |
---|
| 2496 | WM, IWM, F, JAC) |
---|
| 2497 | !WM, IWM, F, JAC, PJAC, SLVS) |
---|
| 2498 | !***BEGIN PROLOGUE DSTODE |
---|
| 2499 | !***SUBSIDIARY |
---|
| 2500 | !***PURPOSE Performs one step of an ODEPACK integration. |
---|
| 2501 | !***TYPE KPP_REAL (SSTODE-S, DSTODE-D) |
---|
| 2502 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 2503 | !***DESCRIPTION |
---|
| 2504 | ! |
---|
| 2505 | ! DSTODE performs one step of the integration of an initial value |
---|
| 2506 | ! problem for a system of ordinary differential equations. |
---|
| 2507 | ! Note: DSTODE is independent of the value of the iteration method |
---|
| 2508 | ! indicator MITER, when this is .ne. 0, and hence is independent |
---|
| 2509 | ! of the type of chord method used, or the Jacobian structure. |
---|
| 2510 | ! Communication with DSTODE is done with the following variables: |
---|
| 2511 | ! |
---|
| 2512 | ! NEQ = integer array containing problem size in NEQ, and |
---|
| 2513 | ! passed as the NEQ argument in all calls to F and JAC. |
---|
| 2514 | ! Y = an array of length .ge. N used as the Y argument in |
---|
| 2515 | ! all calls to F and JAC. |
---|
| 2516 | ! YH = an NYH by LMAX array containing the dependent variables |
---|
| 2517 | ! and their approximate scaled derivatives, where |
---|
| 2518 | ! LMAX = MAXORD + 1. YH(i,j+1) contains the approximate |
---|
| 2519 | ! j-th derivative of y(i), scaled by h**j/factorial(j) |
---|
| 2520 | ! (j = 0,1,...,NQ). on entry for the first step, the first |
---|
| 2521 | ! two columns of YH must be set from the initial values. |
---|
| 2522 | ! NYH = a constant integer .ge. N, the first dimension of YH. |
---|
| 2523 | ! YH1 = a one-dimensional array occupying the same space as YH. |
---|
| 2524 | ! EWT = an array of length N containing multiplicative weights |
---|
| 2525 | ! for local error measurements. Local errors in Y(i) are |
---|
| 2526 | ! compared to 1.0/EWT(i) in various error tests. |
---|
| 2527 | ! SAVF = an array of working storage, of length N. |
---|
| 2528 | ! Also used for input of YH(*,MAXORD+2) when JSTART = -1 |
---|
| 2529 | ! and MAXORD .lt. the current order NQ. |
---|
| 2530 | ! ACOR = a work array of length N, used for the accumulated |
---|
| 2531 | ! corrections. On a successful return, ACOR(i) contains |
---|
| 2532 | ! the estimated one-step local error in Y(i). |
---|
| 2533 | ! WM,IWM = real and integer work arrays associated with matrix |
---|
| 2534 | ! operations in chord iteration (MITER .ne. 0). |
---|
| 2535 | ! PJAC = name of routine to evaluate and preprocess Jacobian matrix |
---|
| 2536 | ! and P = I - h*el0*JAC, if a chord method is being used. |
---|
| 2537 | ! SLVS = name of routine to solve linear system in chord iteration. |
---|
| 2538 | ! CCMAX = maximum relative change in h*el0 before PJAC is called. |
---|
| 2539 | ! H = the step size to be attempted on the next step. |
---|
| 2540 | ! H is altered by the error control algorithm during the |
---|
| 2541 | ! problem. H can be either positive or negative, but its |
---|
| 2542 | ! sign must remain constant throughout the problem. |
---|
| 2543 | ! HMIN = the minimum absolute value of the step size h to be used. |
---|
| 2544 | ! HMXI = inverse of the maximum absolute value of h to be used. |
---|
| 2545 | ! HMXI = 0.0 is allowed and corresponds to an infinite hmax. |
---|
| 2546 | ! HMIN and HMXI may be changed at any time, but will not |
---|
| 2547 | ! take effect until the next change of h is considered. |
---|
| 2548 | ! TN = the independent variable. TN is updated on each step taken. |
---|
| 2549 | ! JSTART = an integer used for input only, with the following |
---|
| 2550 | ! values and meanings: |
---|
| 2551 | ! 0 perform the first step. |
---|
| 2552 | ! .gt.0 take a new step continuing from the last. |
---|
| 2553 | ! -1 take the next step with a new value of H, MAXORD, |
---|
| 2554 | ! N, METH, MITER, and/or matrix parameters. |
---|
| 2555 | ! -2 take the next step with a new value of H, |
---|
| 2556 | ! but with other inputs unchanged. |
---|
| 2557 | ! On return, JSTART is set to 1 to facilitate continuation. |
---|
| 2558 | ! KFLAG = a completion code with the following meanings: |
---|
| 2559 | ! 0 the step was succesful. |
---|
| 2560 | ! -1 the requested error could not be achieved. |
---|
| 2561 | ! -2 corrector convergence could not be achieved. |
---|
| 2562 | ! -3 fatal error in PJAC or SLVS. |
---|
| 2563 | ! A return with KFLAG = -1 or -2 means either |
---|
| 2564 | ! abs(H) = HMIN or 10 consecutive failures occurred. |
---|
| 2565 | ! On a return with KFLAG negative, the values of TN and |
---|
| 2566 | ! the YH array are as of the beginning of the last |
---|
| 2567 | ! step, and H is the last step size attempted. |
---|
| 2568 | ! MAXORD = the maximum order of integration method to be allowed. |
---|
| 2569 | ! MAXCOR = the maximum number of corrector iterations allowed. |
---|
| 2570 | ! MSBP = maximum number of steps between PJAC calls (MITER .gt. 0). |
---|
| 2571 | ! MXNCF = maximum number of convergence failures allowed. |
---|
| 2572 | ! METH/MITER = the method flags. See description in driver. |
---|
| 2573 | ! N = the number of first-order differential equations. |
---|
| 2574 | ! The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD, |
---|
| 2575 | ! MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON. |
---|
| 2576 | ! |
---|
| 2577 | !***SEE ALSO DLSODE |
---|
| 2578 | !***ROUTINES CALLED DCFODE, DVNORM |
---|
| 2579 | !***COMMON BLOCKS DLS001 |
---|
| 2580 | !***REVISION HISTORY (YYMMDD) |
---|
| 2581 | ! 791129 DATE WRITTEN |
---|
| 2582 | ! 890501 Modified prologue to SLATEC/LDOC format. (FNF) |
---|
| 2583 | ! 890503 Minor cosmetic changes. (FNF) |
---|
| 2584 | ! 930809 Renamed to allow single/double precision versions. (ACH) |
---|
| 2585 | ! 010418 Reduced size of Common block /DLS001/. (ACH) |
---|
| 2586 | ! 031105 Restored 'own' variables to Common block /DLS001/, to |
---|
| 2587 | ! enable interrupt/restart feature. (ACH) |
---|
| 2588 | !***END PROLOGUE DSTODE |
---|
| 2589 | !**End |
---|
| 2590 | EXTERNAL F, JAC !, PJAC, SLVS |
---|
| 2591 | INTEGER NEQ, NYH, IWM(*) |
---|
| 2592 | KPP_REAL Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*), & |
---|
| 2593 | ACOR(*), WM(*) |
---|
| 2594 | INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, & |
---|
| 2595 | ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, & |
---|
| 2596 | LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, & |
---|
| 2597 | MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
---|
| 2598 | INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ |
---|
| 2599 | KPP_REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, & |
---|
| 2600 | CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND |
---|
| 2601 | KPP_REAL DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, & |
---|
| 2602 | EXUP,R, RH, RHDN, RHSM, RHUP, TOLD |
---|
| 2603 | !KPP_REAL DVNORM |
---|
| 2604 | COMMON /DLS001/ CONIT, CRATE, EL(13), ELCO(13,12), & |
---|
| 2605 | HOLD, RMAX, TESCO(3,12), & |
---|
| 2606 | CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, & |
---|
| 2607 | IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, & |
---|
| 2608 | ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, & |
---|
| 2609 | LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, & |
---|
| 2610 | MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU |
---|
| 2611 | ! |
---|
| 2612 | !***FIRST EXECUTABLE STATEMENT DSTODE |
---|
| 2613 | KFLAG = 0 |
---|
| 2614 | TOLD = TN |
---|
| 2615 | NCF = 0 |
---|
| 2616 | IERPJ = 0 |
---|
| 2617 | IERSL = 0 |
---|
| 2618 | JCUR = 0 |
---|
| 2619 | ICF = 0 |
---|
| 2620 | DELP = 0.0D0 |
---|
| 2621 | IF (JSTART .GT. 0) GO TO 200 |
---|
| 2622 | IF (JSTART .EQ. -1) GO TO 100 |
---|
| 2623 | IF (JSTART .EQ. -2) GO TO 160 |
---|
| 2624 | !----------------------------------------------------------------------- |
---|
| 2625 | ! On the first call, the order is set to 1, and other variables are |
---|
| 2626 | ! initialized. RMAX is the maximum ratio by which H can be increased |
---|
| 2627 | ! in a single step. It is initially 1.E4 to compensate for the small |
---|
| 2628 | ! initial H, but then is normally equal to 10. If a failure |
---|
| 2629 | ! occurs (in corrector convergence or error test), RMAX is set to 2 |
---|
| 2630 | ! for the next increase. |
---|
| 2631 | !----------------------------------------------------------------------- |
---|
| 2632 | LMAX = MAXORD + 1 |
---|
| 2633 | NQ = 1 |
---|
| 2634 | L = 2 |
---|
| 2635 | IALTH = 2 |
---|
| 2636 | RMAX = 10000.0D0 |
---|
| 2637 | RC = 0.0D0 |
---|
| 2638 | EL0 = 1.0D0 |
---|
| 2639 | CRATE = 0.7D0 |
---|
| 2640 | HOLD = H |
---|
| 2641 | MEO = METH |
---|
| 2642 | NSLP = 0 |
---|
| 2643 | IPUP = MITER |
---|
| 2644 | IRET = 3 |
---|
| 2645 | GO TO 140 |
---|
| 2646 | !----------------------------------------------------------------------- |
---|
| 2647 | ! The following block handles preliminaries needed when JSTART = -1. |
---|
| 2648 | ! IPUP is set to MITER to force a matrix update. |
---|
| 2649 | ! If an order increase is about to be considered (IALTH = 1), |
---|
| 2650 | ! IALTH is reset to 2 to postpone consideration one more step. |
---|
| 2651 | ! If the caller has changed METH, DCFODE is called to reset |
---|
| 2652 | ! the coefficients of the method. |
---|
| 2653 | ! If the caller has changed MAXORD to a value less than the current |
---|
| 2654 | ! order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly. |
---|
| 2655 | ! If H is to be changed, YH must be rescaled. |
---|
| 2656 | ! If H or METH is being changed, IALTH is reset to L = NQ + 1 |
---|
| 2657 | ! to prevent further changes in H for that many steps. |
---|
| 2658 | !----------------------------------------------------------------------- |
---|
| 2659 | 100 IPUP = MITER |
---|
| 2660 | LMAX = MAXORD + 1 |
---|
| 2661 | IF (IALTH .EQ. 1) IALTH = 2 |
---|
| 2662 | IF (METH .EQ. MEO) GO TO 110 |
---|
| 2663 | CALL DCFODE (METH, ELCO, TESCO) |
---|
| 2664 | MEO = METH |
---|
| 2665 | IF (NQ .GT. MAXORD) GO TO 120 |
---|
| 2666 | IALTH = L |
---|
| 2667 | IRET = 1 |
---|
| 2668 | GO TO 150 |
---|
| 2669 | 110 IF (NQ .LE. MAXORD) GO TO 160 |
---|
| 2670 | 120 NQ = MAXORD |
---|
| 2671 | L = LMAX |
---|
| 2672 | DO 125 I = 1,L |
---|
| 2673 | 125 EL(I) = ELCO(I,NQ) |
---|
| 2674 | NQNYH = NQ*NYH |
---|
| 2675 | RC = RC*EL(1)/EL0 |
---|
| 2676 | EL0 = EL(1) |
---|
| 2677 | CONIT = 0.5D0/(NQ+2) |
---|
| 2678 | DDN = DVNORM (N, SAVF, EWT)/TESCO(1,L) |
---|
| 2679 | EXDN = 1.0D0/L |
---|
| 2680 | RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) |
---|
| 2681 | RH = MIN(RHDN,1.0D0) |
---|
| 2682 | IREDO = 3 |
---|
| 2683 | IF (H .EQ. HOLD) GO TO 170 |
---|
| 2684 | RH = MIN(RH,ABS(H/HOLD)) |
---|
| 2685 | H = HOLD |
---|
| 2686 | GO TO 175 |
---|
| 2687 | !----------------------------------------------------------------------- |
---|
| 2688 | ! DCFODE is called to get all the integration coefficients for the |
---|
| 2689 | ! current METH. Then the EL vector and related constants are reset |
---|
| 2690 | ! whenever the order NQ is changed, or at the start of the problem. |
---|
| 2691 | !----------------------------------------------------------------------- |
---|
| 2692 | 140 CALL DCFODE (METH, ELCO, TESCO) |
---|
| 2693 | 150 DO 155 I = 1,L |
---|
| 2694 | 155 EL(I) = ELCO(I,NQ) |
---|
| 2695 | NQNYH = NQ*NYH |
---|
| 2696 | RC = RC*EL(1)/EL0 |
---|
| 2697 | EL0 = EL(1) |
---|
| 2698 | CONIT = 0.5D0/(NQ+2) |
---|
| 2699 | GO TO (160, 170, 200), IRET |
---|
| 2700 | !----------------------------------------------------------------------- |
---|
| 2701 | ! If H is being changed, the H ratio RH is checked against |
---|
| 2702 | ! RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to |
---|
| 2703 | ! L = NQ + 1 to prevent a change of H for that many steps, unless |
---|
| 2704 | ! forced by a convergence or error test failure. |
---|
| 2705 | !----------------------------------------------------------------------- |
---|
| 2706 | 160 IF (H .EQ. HOLD) GO TO 200 |
---|
| 2707 | RH = H/HOLD |
---|
| 2708 | H = HOLD |
---|
| 2709 | IREDO = 3 |
---|
| 2710 | GO TO 175 |
---|
| 2711 | 170 RH = MAX(RH,HMIN/ABS(H)) |
---|
| 2712 | 175 RH = MIN(RH,RMAX) |
---|
| 2713 | RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH) |
---|
| 2714 | R = 1.0D0 |
---|
| 2715 | DO 180 J = 2,L |
---|
| 2716 | R = R*RH |
---|
| 2717 | DO 180 I = 1,N |
---|
| 2718 | 180 YH(I,J) = YH(I,J)*R |
---|
| 2719 | H = H*RH |
---|
| 2720 | RC = RC*RH |
---|
| 2721 | IALTH = L |
---|
| 2722 | IF (IREDO .EQ. 0) GO TO 690 |
---|
| 2723 | !----------------------------------------------------------------------- |
---|
| 2724 | ! This section computes the predicted values by effectively |
---|
| 2725 | ! multiplying the YH array by the Pascal Triangle matrix. |
---|
| 2726 | ! RC is the ratio of new to old values of the coefficient H*EL(1). |
---|
| 2727 | ! When RC differs from 1 by more than CCMAX, IPUP is set to MITER |
---|
| 2728 | ! to force PJAC to be called, if a Jacobian is involved. |
---|
| 2729 | ! In any case, PJAC is called at least every MSBP steps. |
---|
| 2730 | !----------------------------------------------------------------------- |
---|
| 2731 | 200 IF (ABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER |
---|
| 2732 | IF (NST .GE. NSLP+MSBP) IPUP = MITER |
---|
| 2733 | TN = TN + H |
---|
| 2734 | I1 = NQNYH + 1 |
---|
| 2735 | DO 215 JB = 1,NQ |
---|
| 2736 | I1 = I1 - NYH |
---|
| 2737 | !dir$ ivdep |
---|
| 2738 | DO 210 I = I1,NQNYH |
---|
| 2739 | 210 YH1(I) = YH1(I) + YH1(I+NYH) |
---|
| 2740 | 215 CONTINUE |
---|
| 2741 | !----------------------------------------------------------------------- |
---|
| 2742 | ! Up to MAXCOR corrector iterations are taken. A convergence test is |
---|
| 2743 | ! made on the R.M.S. norm of each correction, weighted by the error |
---|
| 2744 | ! weight vector EWT. The sum of the corrections is accumulated in the |
---|
| 2745 | ! vector ACOR(i). The YH array is not altered in the corrector loop. |
---|
| 2746 | !----------------------------------------------------------------------- |
---|
| 2747 | 220 M = 0 |
---|
| 2748 | DO 230 I = 1,N |
---|
| 2749 | 230 Y(I) = YH(I,1) |
---|
| 2750 | CALL F (NEQ, TN, Y, SAVF) |
---|
| 2751 | NFE = NFE + 1 |
---|
| 2752 | IF (IPUP .LE. 0) GO TO 250 |
---|
| 2753 | !----------------------------------------------------------------------- |
---|
| 2754 | ! If indicated, the matrix P = I - h*el(1)*J is reevaluated and |
---|
| 2755 | ! preprocessed before starting the corrector iteration. IPUP is set |
---|
| 2756 | ! to 0 as an indicator that this has been done. |
---|
| 2757 | !----------------------------------------------------------------------- |
---|
| 2758 | !CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC) |
---|
| 2759 | CALL DPREPJ(NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC) |
---|
| 2760 | IPUP = 0 |
---|
| 2761 | RC = 1.0D0 |
---|
| 2762 | NSLP = NST |
---|
| 2763 | CRATE = 0.7D0 |
---|
| 2764 | IF (IERPJ .NE. 0) GO TO 430 |
---|
| 2765 | 250 DO 260 I = 1,N |
---|
| 2766 | 260 ACOR(I) = 0.0D0 |
---|
| 2767 | 270 IF (MITER .NE. 0) GO TO 350 |
---|
| 2768 | !----------------------------------------------------------------------- |
---|
| 2769 | ! In the case of functional iteration, update Y directly from |
---|
| 2770 | ! the result of the last function evaluation. |
---|
| 2771 | !----------------------------------------------------------------------- |
---|
| 2772 | DO 290 I = 1,N |
---|
| 2773 | SAVF(I) = H*SAVF(I) - YH(I,2) |
---|
| 2774 | 290 Y(I) = SAVF(I) - ACOR(I) |
---|
| 2775 | DEL = DVNORM (N, Y, EWT) |
---|
| 2776 | DO 300 I = 1,N |
---|
| 2777 | Y(I) = YH(I,1) + EL(1)*SAVF(I) |
---|
| 2778 | 300 ACOR(I) = SAVF(I) |
---|
| 2779 | GO TO 400 |
---|
| 2780 | !----------------------------------------------------------------------- |
---|
| 2781 | ! In the case of the chord method, compute the corrector error, |
---|
| 2782 | ! and solve the linear system with that as right-hand side and |
---|
| 2783 | ! P as coefficient matrix. |
---|
| 2784 | !----------------------------------------------------------------------- |
---|
| 2785 | 350 DO 360 I = 1,N |
---|
| 2786 | 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I)) |
---|
| 2787 | !CALL SLVS (WM, IWM, Y, SAVF) |
---|
| 2788 | CALL DSOLSY(WM, IWM, Y, SAVF) |
---|
| 2789 | IF (IERSL .LT. 0) GO TO 430 |
---|
| 2790 | IF (IERSL .GT. 0) GO TO 410 |
---|
| 2791 | DEL = DVNORM (N, Y, EWT) |
---|
| 2792 | DO 380 I = 1,N |
---|
| 2793 | ACOR(I) = ACOR(I) + Y(I) |
---|
| 2794 | 380 Y(I) = YH(I,1) + EL(1)*ACOR(I) |
---|
| 2795 | !----------------------------------------------------------------------- |
---|
| 2796 | ! Test for convergence. If M.gt.0, an estimate of the convergence |
---|
| 2797 | ! rate constant is stored in CRATE, and this is used in the test. |
---|
| 2798 | !----------------------------------------------------------------------- |
---|
| 2799 | 400 IF (M .NE. 0) CRATE = MAX(0.2D0*CRATE,DEL/DELP) |
---|
| 2800 | DCON = DEL*MIN(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT) |
---|
| 2801 | IF (DCON .LE. 1.0D0) GO TO 450 |
---|
| 2802 | M = M + 1 |
---|
| 2803 | IF (M .EQ. MAXCOR) GO TO 410 |
---|
| 2804 | IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410 |
---|
| 2805 | DELP = DEL |
---|
| 2806 | CALL F (NEQ, TN, Y, SAVF) |
---|
| 2807 | NFE = NFE + 1 |
---|
| 2808 | GO TO 270 |
---|
| 2809 | !----------------------------------------------------------------------- |
---|
| 2810 | ! The corrector iteration failed to converge. |
---|
| 2811 | ! If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for |
---|
| 2812 | ! the next try. Otherwise the YH array is retracted to its values |
---|
| 2813 | ! before prediction, and H is reduced, if possible. If H cannot be |
---|
| 2814 | ! reduced or MXNCF failures have occurred, exit with KFLAG = -2. |
---|
| 2815 | !----------------------------------------------------------------------- |
---|
| 2816 | 410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430 |
---|
| 2817 | ICF = 1 |
---|
| 2818 | IPUP = MITER |
---|
| 2819 | GO TO 220 |
---|
| 2820 | 430 ICF = 2 |
---|
| 2821 | NCF = NCF + 1 |
---|
| 2822 | RMAX = 2.0D0 |
---|
| 2823 | TN = TOLD |
---|
| 2824 | I1 = NQNYH + 1 |
---|
| 2825 | DO 445 JB = 1,NQ |
---|
| 2826 | I1 = I1 - NYH |
---|
| 2827 | !dir$ ivdep |
---|
| 2828 | DO 440 I = I1,NQNYH |
---|
| 2829 | 440 YH1(I) = YH1(I) - YH1(I+NYH) |
---|
| 2830 | 445 CONTINUE |
---|
| 2831 | IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680 |
---|
| 2832 | IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 670 |
---|
| 2833 | IF (NCF .EQ. MXNCF) GO TO 670 |
---|
| 2834 | RH = 0.25D0 |
---|
| 2835 | IPUP = MITER |
---|
| 2836 | IREDO = 1 |
---|
| 2837 | GO TO 170 |
---|
| 2838 | !----------------------------------------------------------------------- |
---|
| 2839 | ! The corrector has converged. JCUR is set to 0 |
---|
| 2840 | ! to signal that the Jacobian involved may need updating later. |
---|
| 2841 | ! The local error test is made and control passes to statement 500 |
---|
| 2842 | ! if it fails. |
---|
| 2843 | !----------------------------------------------------------------------- |
---|
| 2844 | 450 JCUR = 0 |
---|
| 2845 | IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) |
---|
| 2846 | IF (M .GT. 0) DSM = DVNORM (N, ACOR, EWT)/TESCO(2,NQ) |
---|
| 2847 | IF (DSM .GT. 1.0D0) GO TO 500 |
---|
| 2848 | !----------------------------------------------------------------------- |
---|
| 2849 | ! After a successful step, update the YH array. |
---|
| 2850 | ! Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1. |
---|
| 2851 | ! If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for |
---|
| 2852 | ! use in a possible order increase on the next step. |
---|
| 2853 | ! If a change in H is considered, an increase or decrease in order |
---|
| 2854 | ! by one is considered also. A change in H is made only if it is by a |
---|
| 2855 | ! factor of at least 1.1. If not, IALTH is set to 3 to prevent |
---|
| 2856 | ! testing for that many steps. |
---|
| 2857 | !----------------------------------------------------------------------- |
---|
| 2858 | KFLAG = 0 |
---|
| 2859 | IREDO = 0 |
---|
| 2860 | NST = NST + 1 |
---|
| 2861 | HU = H |
---|
| 2862 | NQU = NQ |
---|
| 2863 | DO 470 J = 1,L |
---|
| 2864 | DO 470 I = 1,N |
---|
| 2865 | 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I) |
---|
| 2866 | IALTH = IALTH - 1 |
---|
| 2867 | IF (IALTH .EQ. 0) GO TO 520 |
---|
| 2868 | IF (IALTH .GT. 1) GO TO 700 |
---|
| 2869 | IF (L .EQ. LMAX) GO TO 700 |
---|
| 2870 | DO 490 I = 1,N |
---|
| 2871 | 490 YH(I,LMAX) = ACOR(I) |
---|
| 2872 | GO TO 700 |
---|
| 2873 | !----------------------------------------------------------------------- |
---|
| 2874 | ! The error test failed. KFLAG keeps track of multiple failures. |
---|
| 2875 | ! Restore TN and the YH array to their previous values, and prepare |
---|
| 2876 | ! to try the step again. Compute the optimum step size for this or |
---|
| 2877 | ! one lower order. After 2 or more failures, H is forced to decrease |
---|
| 2878 | ! by a factor of 0.2 or less. |
---|
| 2879 | !----------------------------------------------------------------------- |
---|
| 2880 | 500 KFLAG = KFLAG - 1 |
---|
| 2881 | TN = TOLD |
---|
| 2882 | I1 = NQNYH + 1 |
---|
| 2883 | DO 515 JB = 1,NQ |
---|
| 2884 | I1 = I1 - NYH |
---|
| 2885 | !dir$ ivdep |
---|
| 2886 | DO 510 I = I1,NQNYH |
---|
| 2887 | 510 YH1(I) = YH1(I) - YH1(I+NYH) |
---|
| 2888 | 515 CONTINUE |
---|
| 2889 | RMAX = 2.0D0 |
---|
| 2890 | IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 660 |
---|
| 2891 | IF (KFLAG .LE. -3) GO TO 640 |
---|
| 2892 | IREDO = 2 |
---|
| 2893 | RHUP = 0.0D0 |
---|
| 2894 | GO TO 540 |
---|
| 2895 | !----------------------------------------------------------------------- |
---|
| 2896 | ! Regardless of the success or failure of the step, factors |
---|
| 2897 | ! RHDN, RHSM, and RHUP are computed, by which H could be multiplied |
---|
| 2898 | ! at order NQ - 1, order NQ, or order NQ + 1, respectively. |
---|
| 2899 | ! In the case of failure, RHUP = 0.0 to avoid an order increase. |
---|
| 2900 | ! The largest of these is determined and the new order chosen |
---|
| 2901 | ! accordingly. If the order is to be increased, we compute one |
---|
| 2902 | ! additional scaled derivative. |
---|
| 2903 | !----------------------------------------------------------------------- |
---|
| 2904 | 520 RHUP = 0.0D0 |
---|
| 2905 | IF (L .EQ. LMAX) GO TO 540 |
---|
| 2906 | DO 530 I = 1,N |
---|
| 2907 | 530 SAVF(I) = ACOR(I) - YH(I,LMAX) |
---|
| 2908 | DUP = DVNORM (N, SAVF, EWT)/TESCO(3,NQ) |
---|
| 2909 | EXUP = 1.0D0/(L+1) |
---|
| 2910 | RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0) |
---|
| 2911 | 540 EXSM = 1.0D0/L |
---|
| 2912 | RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0) |
---|
| 2913 | RHDN = 0.0D0 |
---|
| 2914 | IF (NQ .EQ. 1) GO TO 560 |
---|
| 2915 | DDN = DVNORM (N, YH(1,L), EWT)/TESCO(1,NQ) |
---|
| 2916 | EXDN = 1.0D0/NQ |
---|
| 2917 | RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) |
---|
| 2918 | 560 IF (RHSM .GE. RHUP) GO TO 570 |
---|
| 2919 | IF (RHUP .GT. RHDN) GO TO 590 |
---|
| 2920 | GO TO 580 |
---|
| 2921 | 570 IF (RHSM .LT. RHDN) GO TO 580 |
---|
| 2922 | NEWQ = NQ |
---|
| 2923 | RH = RHSM |
---|
| 2924 | GO TO 620 |
---|
| 2925 | 580 NEWQ = NQ - 1 |
---|
| 2926 | RH = RHDN |
---|
| 2927 | IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0 |
---|
| 2928 | GO TO 620 |
---|
| 2929 | 590 NEWQ = L |
---|
| 2930 | RH = RHUP |
---|
| 2931 | IF (RH .LT. 1.1D0) GO TO 610 |
---|
| 2932 | R = EL(L)/L |
---|
| 2933 | DO 600 I = 1,N |
---|
| 2934 | 600 YH(I,NEWQ+1) = ACOR(I)*R |
---|
| 2935 | GO TO 630 |
---|
| 2936 | 610 IALTH = 3 |
---|
| 2937 | GO TO 700 |
---|
| 2938 | 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610 |
---|
| 2939 | IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0) |
---|
| 2940 | !----------------------------------------------------------------------- |
---|
| 2941 | ! If there is a change of order, reset NQ, l, and the coefficients. |
---|
| 2942 | ! In any case H is reset according to RH and the YH array is rescaled. |
---|
| 2943 | ! Then exit from 690 if the step was OK, or redo the step otherwise. |
---|
| 2944 | !----------------------------------------------------------------------- |
---|
| 2945 | IF (NEWQ .EQ. NQ) GO TO 170 |
---|
| 2946 | 630 NQ = NEWQ |
---|
| 2947 | L = NQ + 1 |
---|
| 2948 | IRET = 2 |
---|
| 2949 | GO TO 150 |
---|
| 2950 | !----------------------------------------------------------------------- |
---|
| 2951 | ! Control reaches this section if 3 or more failures have occured. |
---|
| 2952 | ! If 10 failures have occurred, exit with KFLAG = -1. |
---|
| 2953 | ! It is assumed that the derivatives that have accumulated in the |
---|
| 2954 | ! YH array have errors of the wrong order. Hence the first |
---|
| 2955 | ! derivative is recomputed, and the order is set to 1. Then |
---|
| 2956 | ! H is reduced by a factor of 10, and the step is retried, |
---|
| 2957 | ! until it succeeds or H reaches HMIN. |
---|
| 2958 | !----------------------------------------------------------------------- |
---|
| 2959 | 640 IF (KFLAG .EQ. -10) GO TO 660 |
---|
| 2960 | RH = 0.1D0 |
---|
| 2961 | RH = MAX(HMIN/ABS(H),RH) |
---|
| 2962 | H = H*RH |
---|
| 2963 | DO 645 I = 1,N |
---|
| 2964 | 645 Y(I) = YH(I,1) |
---|
| 2965 | CALL F (NEQ, TN, Y, SAVF) |
---|
| 2966 | NFE = NFE + 1 |
---|
| 2967 | DO 650 I = 1,N |
---|
| 2968 | 650 YH(I,2) = H*SAVF(I) |
---|
| 2969 | IPUP = MITER |
---|
| 2970 | IALTH = 5 |
---|
| 2971 | IF (NQ .EQ. 1) GO TO 200 |
---|
| 2972 | NQ = 1 |
---|
| 2973 | L = 2 |
---|
| 2974 | IRET = 3 |
---|
| 2975 | GO TO 150 |
---|
| 2976 | !----------------------------------------------------------------------- |
---|
| 2977 | ! All returns are made through this section. H is saved in HOLD |
---|
| 2978 | ! to allow the caller to change H on the next step. |
---|
| 2979 | !----------------------------------------------------------------------- |
---|
| 2980 | 660 KFLAG = -1 |
---|
| 2981 | GO TO 720 |
---|
| 2982 | 670 KFLAG = -2 |
---|
| 2983 | GO TO 720 |
---|
| 2984 | 680 KFLAG = -3 |
---|
| 2985 | GO TO 720 |
---|
| 2986 | 690 RMAX = 10.0D0 |
---|
| 2987 | 700 R = 1.0D0/TESCO(2,NQU) |
---|
| 2988 | DO 710 I = 1,N |
---|
| 2989 | 710 ACOR(I) = ACOR(I)*R |
---|
| 2990 | 720 HOLD = H |
---|
| 2991 | JSTART = 1 |
---|
| 2992 | RETURN |
---|
| 2993 | !----------------------- END OF SUBROUTINE DSTODE ---------------------- |
---|
| 2994 | END SUBROUTINE DSTODE |
---|
| 2995 | !DECK DEWSET |
---|
| 2996 | SUBROUTINE DEWSET (N, ITOL, RelTol, AbsTol, YCUR, EWT) |
---|
| 2997 | !***BEGIN PROLOGUE DEWSET |
---|
| 2998 | !***SUBSIDIARY |
---|
| 2999 | !***PURPOSE Set error weight vector. |
---|
| 3000 | !***TYPE KPP_REAL (SEWSET-S, DEWSET-D) |
---|
| 3001 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 3002 | !***DESCRIPTION |
---|
| 3003 | ! |
---|
| 3004 | ! This subroutine sets the error weight vector EWT according to |
---|
| 3005 | ! EWT(i) = RelTol(i)*ABS(YCUR(i)) + AbsTol(i), i = 1,...,N, |
---|
| 3006 | ! with the subscript on RelTol and/or AbsTol possibly replaced by 1 above, |
---|
| 3007 | ! depending on the value of ITOL. |
---|
| 3008 | ! |
---|
| 3009 | !***SEE ALSO DLSODE |
---|
| 3010 | !***ROUTINES CALLED (NONE) |
---|
| 3011 | !***REVISION HISTORY (YYMMDD) |
---|
| 3012 | ! 791129 DATE WRITTEN |
---|
| 3013 | ! 890501 Modified prologue to SLATEC/LDOC format. (FNF) |
---|
| 3014 | ! 890503 Minor cosmetic changes. (FNF) |
---|
| 3015 | ! 930809 Renamed to allow single/double precision versions. (ACH) |
---|
| 3016 | !***END PROLOGUE DEWSET |
---|
| 3017 | !**End |
---|
| 3018 | INTEGER N, ITOL |
---|
| 3019 | INTEGER I |
---|
| 3020 | KPP_REAL RelTol(*), AbsTol(*), YCUR(N), EWT(N) |
---|
| 3021 | ! |
---|
| 3022 | !***FIRST EXECUTABLE STATEMENT DEWSET |
---|
| 3023 | GO TO (10, 20, 30, 40), ITOL |
---|
| 3024 | 10 CONTINUE |
---|
| 3025 | DO 15 I = 1,N |
---|
| 3026 | 15 EWT(I) = RelTol(1)*ABS(YCUR(I)) + AbsTol(1) |
---|
| 3027 | RETURN |
---|
| 3028 | 20 CONTINUE |
---|
| 3029 | DO 25 I = 1,N |
---|
| 3030 | 25 EWT(I) = RelTol(1)*ABS(YCUR(I)) + AbsTol(I) |
---|
| 3031 | RETURN |
---|
| 3032 | 30 CONTINUE |
---|
| 3033 | DO 35 I = 1,N |
---|
| 3034 | 35 EWT(I) = RelTol(I)*ABS(YCUR(I)) + AbsTol(1) |
---|
| 3035 | RETURN |
---|
| 3036 | 40 CONTINUE |
---|
| 3037 | DO 45 I = 1,N |
---|
| 3038 | 45 EWT(I) = RelTol(I)*ABS(YCUR(I)) + AbsTol(I) |
---|
| 3039 | RETURN |
---|
| 3040 | !----------------------- END OF SUBROUTINE DEWSET ---------------------- |
---|
| 3041 | END SUBROUTINE DEWSET |
---|
| 3042 | !DECK DVNORM |
---|
| 3043 | KPP_REAL FUNCTION DVNORM (N, V, W) |
---|
| 3044 | !***BEGIN PROLOGUE DVNORM |
---|
| 3045 | !***SUBSIDIARY |
---|
| 3046 | !***PURPOSE Weighted root-mean-square vector norm. |
---|
| 3047 | !***TYPE KPP_REAL (SVNORM-S, DVNORM-D) |
---|
| 3048 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 3049 | !***DESCRIPTION |
---|
| 3050 | ! |
---|
| 3051 | ! This function routine computes the weighted root-mean-square norm |
---|
| 3052 | ! of the vector of length N contained in the array V, with weights |
---|
| 3053 | ! contained in the array W of length N: |
---|
| 3054 | ! DVNORM = SQRT( (1/N) * SUM( V(i)*W(i) )**2 ) |
---|
| 3055 | ! |
---|
| 3056 | !***SEE ALSO DLSODE |
---|
| 3057 | !***ROUTINES CALLED (NONE) |
---|
| 3058 | !***REVISION HISTORY (YYMMDD) |
---|
| 3059 | ! 791129 DATE WRITTEN |
---|
| 3060 | ! 890501 Modified prologue to SLATEC/LDOC format. (FNF) |
---|
| 3061 | ! 890503 Minor cosmetic changes. (FNF) |
---|
| 3062 | ! 930809 Renamed to allow single/double precision versions. (ACH) |
---|
| 3063 | !***END PROLOGUE DVNORM |
---|
| 3064 | !**End |
---|
| 3065 | INTEGER N, I |
---|
| 3066 | KPP_REAL V(N), W(N), SUM |
---|
| 3067 | ! |
---|
| 3068 | !***FIRST EXECUTABLE STATEMENT DVNORM |
---|
| 3069 | SUM = 0.0D0 |
---|
| 3070 | DO 10 I = 1,N |
---|
| 3071 | 10 SUM = SUM + (V(I)*W(I))**2 |
---|
| 3072 | DVNORM = SQRT(SUM/N) |
---|
| 3073 | RETURN |
---|
| 3074 | !----------------------- END OF FUNCTION DVNORM ------------------------ |
---|
| 3075 | END FUNCTION DVNORM |
---|
| 3076 | !DECK XERRWD |
---|
| 3077 | SUBROUTINE XERRWD (MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2) |
---|
| 3078 | !***BEGIN PROLOGUE XERRWD |
---|
| 3079 | !***SUBSIDIARY |
---|
| 3080 | !***PURPOSE Write error message with values. |
---|
| 3081 | !***CATEGORY R3C |
---|
| 3082 | !***TYPE KPP_REAL (XERRWV-S, XERRWD-D) |
---|
| 3083 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 3084 | !***DESCRIPTION |
---|
| 3085 | ! |
---|
| 3086 | ! Subroutines XERRWD, XSETF, XSETUN, and the function routine IXSAV, |
---|
| 3087 | ! as given here, constitute a simplified version of the SLATEC error |
---|
| 3088 | ! handling package. |
---|
| 3089 | ! |
---|
| 3090 | ! All arguments are input arguments. |
---|
| 3091 | ! |
---|
| 3092 | ! MSG = The message (character array). |
---|
| 3093 | ! NMES = The length of MSG (number of characters). |
---|
| 3094 | ! NERR = The error number (not used). |
---|
| 3095 | ! LEVEL = The error level.. |
---|
| 3096 | ! 0 or 1 means recoverable (control returns to caller). |
---|
| 3097 | ! 2 means fatal (run is aborted--see note below). |
---|
| 3098 | ! NI = Number of integers (0, 1, or 2) to be printed with message. |
---|
| 3099 | ! I1,I2 = Integers to be printed, depending on NI. |
---|
| 3100 | ! NR = Number of reals (0, 1, or 2) to be printed with message. |
---|
| 3101 | ! R1,R2 = Reals to be printed, depending on NR. |
---|
| 3102 | ! |
---|
| 3103 | ! Note.. this routine is machine-dependent and specialized for use |
---|
| 3104 | ! in limited context, in the following ways.. |
---|
| 3105 | ! 1. The argument MSG is assumed to be of type CHARACTER, and |
---|
| 3106 | ! the message is printed with a format of (1X,A). |
---|
| 3107 | ! 2. The message is assumed to take only one line. |
---|
| 3108 | ! Multi-line messages are generated by repeated calls. |
---|
| 3109 | ! 3. If LEVEL = 2, control passes to the statement STOP |
---|
| 3110 | ! to abort the run. This statement may be machine-dependent. |
---|
| 3111 | ! 4. R1 and R2 are assumed to be in double precision and are printed |
---|
| 3112 | ! in D21.13 format. |
---|
| 3113 | ! |
---|
| 3114 | !***ROUTINES CALLED IXSAV |
---|
| 3115 | !***REVISION HISTORY (YYMMDD) |
---|
| 3116 | ! 920831 DATE WRITTEN |
---|
| 3117 | ! 921118 Replaced MFLGSV/LUNSAV by IXSAV. (ACH) |
---|
| 3118 | ! 930329 Modified prologue to SLATEC format. (FNF) |
---|
| 3119 | ! 930407 Changed MSG from CHARACTER*1 array to variable. (FNF) |
---|
| 3120 | ! 930922 Minor cosmetic change. (FNF) |
---|
| 3121 | !***END PROLOGUE XERRWD |
---|
| 3122 | ! |
---|
| 3123 | !*Internal Notes: |
---|
| 3124 | ! |
---|
| 3125 | ! For a different default logical unit number, IXSAV (or a subsidiary |
---|
| 3126 | ! routine that it calls) will need to be modified. |
---|
| 3127 | ! For a different run-abort command, change the statement following |
---|
| 3128 | ! statement 100 at the end. |
---|
| 3129 | !----------------------------------------------------------------------- |
---|
| 3130 | ! Subroutines called by XERRWD.. None |
---|
| 3131 | ! Function routine called by XERRWD.. IXSAV |
---|
| 3132 | !----------------------------------------------------------------------- |
---|
| 3133 | !**End |
---|
| 3134 | ! |
---|
| 3135 | ! Declare arguments. |
---|
| 3136 | ! |
---|
| 3137 | KPP_REAL R1, R2 |
---|
| 3138 | INTEGER NMES, NERR, LEVEL, NI, I1, I2, NR |
---|
| 3139 | CHARACTER*(*) MSG |
---|
| 3140 | ! |
---|
| 3141 | ! Declare local variables. |
---|
| 3142 | ! |
---|
| 3143 | INTEGER LUNIT, MESFLG !, IXSAV |
---|
| 3144 | ! |
---|
| 3145 | ! Get logical unit number and message print flag. |
---|
| 3146 | ! |
---|
| 3147 | !***FIRST EXECUTABLE STATEMENT XERRWD |
---|
| 3148 | LUNIT = IXSAV (1, 0, .FALSE.) |
---|
| 3149 | MESFLG = IXSAV (2, 0, .FALSE.) |
---|
| 3150 | IF (MESFLG .EQ. 0) GO TO 100 |
---|
| 3151 | ! |
---|
| 3152 | ! Write the message. |
---|
| 3153 | ! |
---|
| 3154 | WRITE (LUNIT,10) MSG |
---|
| 3155 | 10 FORMAT(1X,A) |
---|
| 3156 | IF (NI .EQ. 1) WRITE (LUNIT, 20) I1 |
---|
| 3157 | 20 FORMAT(6X,'In above message, I1 =',I10) |
---|
| 3158 | IF (NI .EQ. 2) WRITE (LUNIT, 30) I1,I2 |
---|
| 3159 | 30 FORMAT(6X,'In above message, I1 =',I10,3X,'I2 =',I10) |
---|
| 3160 | IF (NR .EQ. 1) WRITE (LUNIT, 40) R1 |
---|
| 3161 | 40 FORMAT(6X,'In above message, R1 =',D21.13) |
---|
| 3162 | IF (NR .EQ. 2) WRITE (LUNIT, 50) R1,R2 |
---|
| 3163 | 50 FORMAT(6X,'In above, R1 =',D21.13,3X,'R2 =',D21.13) |
---|
| 3164 | ! |
---|
| 3165 | ! Abort the run if LEVEL = 2. |
---|
| 3166 | ! |
---|
| 3167 | 100 IF (LEVEL .NE. 2) RETURN |
---|
| 3168 | STOP |
---|
| 3169 | !----------------------- End of Subroutine XERRWD ---------------------- |
---|
| 3170 | END SUBROUTINE XERRWD |
---|
| 3171 | !DECK XSETF |
---|
| 3172 | SUBROUTINE XSETF (MFLAG) |
---|
| 3173 | !***BEGIN PROLOGUE XSETF |
---|
| 3174 | !***PURPOSE Reset the error print control flag. |
---|
| 3175 | !***CATEGORY R3A |
---|
| 3176 | !***TYPE ALL (XSETF-A) |
---|
| 3177 | !***KEYWORDS ERROR CONTROL |
---|
| 3178 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 3179 | !***DESCRIPTION |
---|
| 3180 | ! |
---|
| 3181 | ! XSETF sets the error print control flag to MFLAG: |
---|
| 3182 | ! MFLAG=1 means print all messages (the default). |
---|
| 3183 | ! MFLAG=0 means no printing. |
---|
| 3184 | ! |
---|
| 3185 | !***SEE ALSO XERRWD, XERRWV |
---|
| 3186 | !***REFERENCES (NONE) |
---|
| 3187 | !***ROUTINES CALLED IXSAV |
---|
| 3188 | !***REVISION HISTORY (YYMMDD) |
---|
| 3189 | ! 921118 DATE WRITTEN |
---|
| 3190 | ! 930329 Added SLATEC format prologue. (FNF) |
---|
| 3191 | ! 930407 Corrected SEE ALSO section. (FNF) |
---|
| 3192 | ! 930922 Made user-callable, and other cosmetic changes. (FNF) |
---|
| 3193 | !***END PROLOGUE XSETF |
---|
| 3194 | ! |
---|
| 3195 | ! Subroutines called by XSETF.. None |
---|
| 3196 | ! Function routine called by XSETF.. IXSAV |
---|
| 3197 | !----------------------------------------------------------------------- |
---|
| 3198 | !**End |
---|
| 3199 | INTEGER MFLAG, JUNK !, IXSAV |
---|
| 3200 | ! |
---|
| 3201 | !***FIRST EXECUTABLE STATEMENT XSETF |
---|
| 3202 | IF (MFLAG .EQ. 0 .OR. MFLAG .EQ. 1) JUNK = IXSAV (2,MFLAG,.TRUE.) |
---|
| 3203 | RETURN |
---|
| 3204 | !----------------------- End of Subroutine XSETF ----------------------- |
---|
| 3205 | END SUBROUTINE XSETF |
---|
| 3206 | !DECK XSETUN |
---|
| 3207 | SUBROUTINE XSETUN (LUN) |
---|
| 3208 | !***BEGIN PROLOGUE XSETUN |
---|
| 3209 | !***PURPOSE Reset the logical unit number for error messages. |
---|
| 3210 | !***CATEGORY R3B |
---|
| 3211 | !***TYPE ALL (XSETUN-A) |
---|
| 3212 | !***KEYWORDS ERROR CONTROL |
---|
| 3213 | !***DESCRIPTION |
---|
| 3214 | ! |
---|
| 3215 | ! XSETUN sets the logical unit number for error messages to LUN. |
---|
| 3216 | ! |
---|
| 3217 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 3218 | !***SEE ALSO XERRWD, XERRWV |
---|
| 3219 | !***REFERENCES (NONE) |
---|
| 3220 | !***ROUTINES CALLED IXSAV |
---|
| 3221 | !***REVISION HISTORY (YYMMDD) |
---|
| 3222 | ! 921118 DATE WRITTEN |
---|
| 3223 | ! 930329 Added SLATEC format prologue. (FNF) |
---|
| 3224 | ! 930407 Corrected SEE ALSO section. (FNF) |
---|
| 3225 | ! 930922 Made user-callable, and other cosmetic changes. (FNF) |
---|
| 3226 | !***END PROLOGUE XSETUN |
---|
| 3227 | ! |
---|
| 3228 | ! Subroutines called by XSETUN.. None |
---|
| 3229 | ! Function routine called by XSETUN.. IXSAV |
---|
| 3230 | !----------------------------------------------------------------------- |
---|
| 3231 | !**End |
---|
| 3232 | INTEGER LUN, JUNK !, IXSAV |
---|
| 3233 | ! |
---|
| 3234 | !***FIRST EXECUTABLE STATEMENT XSETUN |
---|
| 3235 | IF (LUN .GT. 0) JUNK = IXSAV (1,LUN,.TRUE.) |
---|
| 3236 | RETURN |
---|
| 3237 | !----------------------- End of Subroutine XSETUN ---------------------- |
---|
| 3238 | END SUBROUTINE XSETUN |
---|
| 3239 | !DECK IXSAV |
---|
| 3240 | INTEGER FUNCTION IXSAV (IPAR, IVALUE, ISET) |
---|
| 3241 | !***BEGIN PROLOGUE IXSAV |
---|
| 3242 | !***SUBSIDIARY |
---|
| 3243 | !***PURPOSE Save and recall error message control parameters. |
---|
| 3244 | !***CATEGORY R3C |
---|
| 3245 | !***TYPE ALL (IXSAV-A) |
---|
| 3246 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 3247 | !***DESCRIPTION |
---|
| 3248 | ! |
---|
| 3249 | ! IXSAV saves and recalls one of two error message parameters: |
---|
| 3250 | ! LUNIT, the logical unit number to which messages are printed, and |
---|
| 3251 | ! MESFLG, the message print flag. |
---|
| 3252 | ! This is a modification of the SLATEC library routine J4SAVE. |
---|
| 3253 | ! |
---|
| 3254 | ! Saved local variables.. |
---|
| 3255 | ! LUNIT = Logical unit number for messages. The default is obtained |
---|
| 3256 | ! by a call to IUMACH (may be machine-dependent). |
---|
| 3257 | ! MESFLG = Print control flag.. |
---|
| 3258 | ! 1 means print all messages (the default). |
---|
| 3259 | ! 0 means no printing. |
---|
| 3260 | ! |
---|
| 3261 | ! On input.. |
---|
| 3262 | ! IPAR = Parameter indicator (1 for LUNIT, 2 for MESFLG). |
---|
| 3263 | ! IVALUE = The value to be set for the parameter, if ISET = .TRUE. |
---|
| 3264 | ! ISET = Logical flag to indicate whether to read or write. |
---|
| 3265 | ! If ISET = .TRUE., the parameter will be given |
---|
| 3266 | ! the value IVALUE. If ISET = .FALSE., the parameter |
---|
| 3267 | ! will be unchanged, and IVALUE is a dummy argument. |
---|
| 3268 | ! |
---|
| 3269 | ! On return.. |
---|
| 3270 | ! IXSAV = The (old) value of the parameter. |
---|
| 3271 | ! |
---|
| 3272 | !***SEE ALSO XERRWD, XERRWV |
---|
| 3273 | !***ROUTINES CALLED IUMACH |
---|
| 3274 | !***REVISION HISTORY (YYMMDD) |
---|
| 3275 | ! 921118 DATE WRITTEN |
---|
| 3276 | ! 930329 Modified prologue to SLATEC format. (FNF) |
---|
| 3277 | ! 930915 Added IUMACH call to get default output unit. (ACH) |
---|
| 3278 | ! 930922 Minor cosmetic changes. (FNF) |
---|
| 3279 | ! 010425 Type declaration for IUMACH added. (ACH) |
---|
| 3280 | !***END PROLOGUE IXSAV |
---|
| 3281 | ! |
---|
| 3282 | ! Subroutines called by IXSAV.. None |
---|
| 3283 | ! Function routine called by IXSAV.. IUMACH |
---|
| 3284 | !----------------------------------------------------------------------- |
---|
| 3285 | !**End |
---|
| 3286 | LOGICAL ISET |
---|
| 3287 | INTEGER IPAR, IVALUE |
---|
| 3288 | !----------------------------------------------------------------------- |
---|
| 3289 | INTEGER LUNIT, MESFLG!, IUMACH |
---|
| 3290 | !----------------------------------------------------------------------- |
---|
| 3291 | ! The following Fortran-77 declaration is to cause the values of the |
---|
| 3292 | ! listed (local) variables to be saved between calls to this routine. |
---|
| 3293 | !----------------------------------------------------------------------- |
---|
| 3294 | SAVE LUNIT, MESFLG |
---|
| 3295 | DATA LUNIT/-1/, MESFLG/1/ |
---|
| 3296 | ! |
---|
| 3297 | !***FIRST EXECUTABLE STATEMENT IXSAV |
---|
| 3298 | IF (IPAR .EQ. 1) THEN |
---|
| 3299 | IF (LUNIT .EQ. -1) LUNIT = IUMACH() |
---|
| 3300 | IXSAV = LUNIT |
---|
| 3301 | IF (ISET) LUNIT = IVALUE |
---|
| 3302 | ENDIF |
---|
| 3303 | ! |
---|
| 3304 | IF (IPAR .EQ. 2) THEN |
---|
| 3305 | IXSAV = MESFLG |
---|
| 3306 | IF (ISET) MESFLG = IVALUE |
---|
| 3307 | ENDIF |
---|
| 3308 | ! |
---|
| 3309 | RETURN |
---|
| 3310 | !----------------------- End of Function IXSAV ------------------------- |
---|
| 3311 | END FUNCTION IXSAV |
---|
| 3312 | !DECK IUMACH |
---|
| 3313 | INTEGER FUNCTION IUMACH() |
---|
| 3314 | !***BEGIN PROLOGUE IUMACH |
---|
| 3315 | !***PURPOSE Provide standard output unit number. |
---|
| 3316 | !***CATEGORY R1 |
---|
| 3317 | !***TYPE INTEGER (IUMACH-I) |
---|
| 3318 | !***KEYWORDS MACHINE CONSTANTS |
---|
| 3319 | !***AUTHOR Hindmarsh, Alan C., (LLNL) |
---|
| 3320 | !***DESCRIPTION |
---|
| 3321 | ! *Usage: |
---|
| 3322 | ! INTEGER LOUT, IUMACH |
---|
| 3323 | ! LOUT = IUMACH() |
---|
| 3324 | ! |
---|
| 3325 | ! *Function Return Values: |
---|
| 3326 | ! LOUT : the standard logical unit for Fortran output. |
---|
| 3327 | ! |
---|
| 3328 | !***REFERENCES (NONE) |
---|
| 3329 | !***ROUTINES CALLED (NONE) |
---|
| 3330 | !***REVISION HISTORY (YYMMDD) |
---|
| 3331 | ! 930915 DATE WRITTEN |
---|
| 3332 | ! 930922 Made user-callable, and other cosmetic changes. (FNF) |
---|
| 3333 | !***END PROLOGUE IUMACH |
---|
| 3334 | ! |
---|
| 3335 | !*Internal Notes: |
---|
| 3336 | ! The built-in value of 6 is standard on a wide range of Fortran |
---|
| 3337 | ! systems. This may be machine-dependent. |
---|
| 3338 | !**End |
---|
| 3339 | !***FIRST EXECUTABLE STATEMENT IUMACH |
---|
| 3340 | IUMACH = 6 |
---|
| 3341 | ! |
---|
| 3342 | RETURN |
---|
| 3343 | !----------------------- End of Function IUMACH ------------------------ |
---|
| 3344 | END FUNCTION IUMACH |
---|
| 3345 | |
---|
| 3346 | !---- END OF SUBROUTINE DLSODE AND ITS INTERNAL PROCEDURES |
---|
| 3347 | END SUBROUTINE DLSODE |
---|
| 3348 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 3349 | SUBROUTINE FUN_CHEM(N, T, V, FCT) |
---|
| 3350 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 3351 | |
---|
| 3352 | USE KPP_ROOT_Parameters |
---|
| 3353 | USE KPP_ROOT_Global |
---|
| 3354 | USE KPP_ROOT_Function, ONLY: Fun |
---|
| 3355 | USE KPP_ROOT_Rates |
---|
| 3356 | |
---|
| 3357 | IMPLICIT NONE |
---|
| 3358 | |
---|
| 3359 | INTEGER :: N |
---|
| 3360 | KPP_REAL :: V(NVAR), FCT(NVAR), T, TOLD |
---|
| 3361 | |
---|
| 3362 | ! TOLD = TIME |
---|
| 3363 | ! TIME = T |
---|
| 3364 | ! CALL Update_SUN() |
---|
| 3365 | ! CALL Update_RCONST() |
---|
| 3366 | ! CALL Update_PHOTO() |
---|
| 3367 | ! TIME = TOLD |
---|
| 3368 | |
---|
| 3369 | CALL Fun(V, FIX, RCONST, FCT) |
---|
| 3370 | |
---|
| 3371 | !Nfun=Nfun+1 |
---|
| 3372 | |
---|
| 3373 | END SUBROUTINE FUN_CHEM |
---|
| 3374 | |
---|
| 3375 | |
---|
| 3376 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 3377 | SUBROUTINE JAC_CHEM (N, T, V, JF) |
---|
| 3378 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 3379 | |
---|
| 3380 | USE KPP_ROOT_Parameters |
---|
| 3381 | USE KPP_ROOT_Global |
---|
| 3382 | USE KPP_ROOT_JacobianSP |
---|
| 3383 | USE KPP_ROOT_Jacobian, ONLY: Jac_SP |
---|
| 3384 | USE KPP_ROOT_Rates |
---|
| 3385 | |
---|
| 3386 | IMPLICIT NONE |
---|
| 3387 | |
---|
| 3388 | KPP_REAL :: V(NVAR), T, TOLD |
---|
| 3389 | INTEGER :: I, J, N, ML, MU, NROWPD |
---|
| 3390 | #ifdef FULL_ALGEBRA |
---|
| 3391 | KPP_REAL :: JV(LU_NONZERO), JF(NVAR,NVAR) |
---|
| 3392 | #else |
---|
| 3393 | KPP_REAL :: JF(LU_NONZERO) |
---|
| 3394 | #endif |
---|
| 3395 | |
---|
| 3396 | ! TOLD = TIME |
---|
| 3397 | ! TIME = T |
---|
| 3398 | ! CALL Update_SUN() |
---|
| 3399 | ! CALL Update_RCONST() |
---|
| 3400 | ! CALL Update_PHOTO() |
---|
| 3401 | ! TIME = TOLD |
---|
| 3402 | |
---|
| 3403 | #ifdef FULL_ALGEBRA |
---|
| 3404 | CALL Jac_SP(V, FIX, RCONST, JV) |
---|
| 3405 | DO j=1,NVAR |
---|
| 3406 | DO i=1,NVAR |
---|
| 3407 | JF(i,j) = 0.0d0 |
---|
| 3408 | END DO |
---|
| 3409 | END DO |
---|
| 3410 | DO i=1,LU_NONZERO |
---|
| 3411 | JF(LU_IROW(i),LU_ICOL(i)) = JV(i) |
---|
| 3412 | END DO |
---|
| 3413 | #else |
---|
| 3414 | CALL Jac_SP(V, FIX, RCONST, JF) |
---|
| 3415 | #endif |
---|
| 3416 | !Njac=Njac+1 |
---|
| 3417 | |
---|
| 3418 | END SUBROUTINE JAC_CHEM |
---|
| 3419 | |
---|
| 3420 | |
---|
| 3421 | END MODULE KPP_ROOT_Integrator |
---|