source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/kpp_lsode.f90 @ 4794

Last change on this file since 4794 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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