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

Last change on this file was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk

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