I believe this was a bug, no idea how it was even working before

[WRF.git] / chem / module_cmu_dvode_solver.F

!**********************************************************************************  
! This computer software was prepared by Battelle Memorial Institute, hereinafter
! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of 
! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
!
! MOSAIC module: see module_mosaic_driver.F for references and terms of use
!**********************************************************************************  

!-----------------------------------------------------------------------
!   module_cmu_dvode_solver - from vode.f & vode_subs.f on 27-oct-2005
!           (vode.f - downloaded from www.netlib.org on 28-jul-2004)
!           (vode_subs.f - created on 28-jul-2004
!               by downloading following from www.netlib.org
!               1.  daxpy, dcopy, ddot, dnrm2, dscal, idamax from blas
!               2.  dgbfa, dbgsl, dgefa, dgesl from linpack)
!
!       first converted to lowercase
!       then converted to fortran-90
!       then converted to module
!           for this step, had to comment out declarations of any
!               functions that are part of the module
!       also changed common block names to reduce potential for conflicts
!           dvod01 --> dvod_cmn01
!           dvod02 --> dvod_cmn02
!-----------------------------------------------------------------------
! 27-oct-2005 - do not declare functions that are in the module
! 04-mar-2008 - eliminated common blocks;
!    halt if istate .ne. 1, because repeated 1-step execution
!       will not work now that common blocks are gone
!-----------------------------------------------------------------------

      module module_cmu_dvode_solver



      type dvode_cmn_vars
         integer ::   &
            icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
            l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
            locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
            n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
            nslp, nyh,   &
            ncfn, netf, nfe, nje, nlu, nni, nqu, nst
         double precision ::   &
            acnrm, ccmxj, conp, crate, drc, el(13),   &
            eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
            rc, rl1, tau(13), tq(5), tn, uround,   &
            hu,   &
            etaq, etaqm1
!        common /dvod_cmn01/ acnrm, ccmxj, conp, crate, drc, el(13),   &
!                     eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!                     rc, rl1, tau(13), tq(5), tn, uround,   &
!                     icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!                     l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!                     locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!                     n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!                     nslp, nyh
!        common /dvod_cmn02/ hu, ncfn, netf, nfe, nje, nlu, nni, nqu, nst
      end type dvode_cmn_vars



      contains



!-----------------------------------------------------------------------
!   vode.f - downloaded from www.netlib.org on 28-jul-2004
!-----------------------------------------------------------------------

!deck dvode
      subroutine dvode (f, neq, y, t, tout, itol, rtol, atol, itask,   &
                  istate, iopt, rwork, lrw, iwork, liw, jac, mf,   &
                  rpar, ipar)
      implicit none
      external f, jac
      double precision y, t, tout, rtol, atol, rwork, rpar
      integer neq, itol, itask, istate, iopt, lrw, iwork, liw,   &
              mf, ipar
      dimension y(*), rtol(*), atol(*), rwork(lrw), iwork(liw),   &
                rpar(*), ipar(*)
!-----------------------------------------------------------------------
! dvode: variable-coefficient ordinary differential equation solver,
! with fixed-leading-coefficient implementation.
! this version is in double precision.
!
! dvode solves the initial value problem for stiff or nonstiff
! systems of first order odes,
!     dy/dt = f(t,y) ,  or, in component form,
!     dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(neq)) (i = 1,...,neq).
! dvode is a package based on the episode and episodeb packages, and
! on the odepack user interface standard, with minor modifications.
!-----------------------------------------------------------------------
! authors:
!               peter n. brown and alan c. hindmarsh
!               center for applied scientific computing, l-561
!               lawrence livermore national laboratory
!               livermore, ca 94551
! and
!               george d. byrne
!               illinois institute of technology
!               chicago, il 60616
!-----------------------------------------------------------------------
! references:
!
! 1. p. n. brown, g. d. byrne, and a. c. hindmarsh, 'vode: a variable
!    coefficient ode solver,' siam j. sci. stat. comput., 10 (1989),
!    pp. 1038-1051.  also, llnl report ucrl-98412, june 1988.
! 2. g. d. byrne and a. c. hindmarsh, 'a polyalgorithm for the
!    numerical solution of ordinary differential equations,'
!    acm trans. math. software, 1 (1975), pp. 71-96.
! 3. a. c. hindmarsh and g. d. byrne, 'episode: an effective package
!    for the integration of systems of ordinary differential
!    equations,' llnl report ucid-30112, rev. 1, april 1977.
! 4. g. d. byrne and a. c. hindmarsh, 'episodeb: an experimental
!    package for the integration of systems of ordinary differential
!    equations with banded jacobians,' llnl report ucid-30132, april
!    1976.
! 5. a. c. hindmarsh, 'odepack, a systematized collection of ode
!    solvers,' in scientific computing, r. s. stepleman et al., eds.,
!    north-holland, amsterdam, 1983, pp. 55-64.
! 6. k. r. jackson and r. sacks-davis, 'an alternative implementation
!    of variable step-size multistep formulas for stiff odes,' acm
!    trans. math. software, 6 (1980), pp. 295-318.
!-----------------------------------------------------------------------
! summary of usage.
!
! communication between the user and the dvode package, for normal
! situations, is summarized here.  this summary describes only a subset
! of the full set of options available.  see the full description for
! details, including optional communication, nonstandard options,
! and instructions for special situations.  see also the example
! problem (with program and output) following this summary.
!
! a. first provide a subroutine of the form:
!           subroutine f (neq, t, y, ydot, rpar, ipar)
!           double precision t, y(neq), ydot(neq), rpar
! which supplies the vector function f by loading ydot(i) with f(i).
!
! b. next determine (or guess) whether or not the problem is stiff.
! stiffness occurs when the jacobian matrix df/dy has an eigenvalue
! whose real part is negative and large in magnitude, compared to the
! reciprocal of the t span of interest.  if the problem is nonstiff,
! use a method flag mf = 10.  if it is stiff, there are four standard
! choices for mf (21, 22, 24, 25), and dvode requires the jacobian
! matrix in some form.  in these cases (mf .gt. 0), dvode will use a
! saved copy of the jacobian matrix.  if this is undesirable because of
! storage limitations, set mf to the corresponding negative value
! (-21, -22, -24, -25).  (see full description of mf below.)
! the jacobian matrix is regarded either as full (mf = 21 or 22),
! or banded (mf = 24 or 25).  in the banded case, dvode requires two
! half-bandwidth parameters ml and mu.  these are, respectively, the
! widths of the lower and upper parts of the band, excluding the main
! diagonal.  thus the band consists of the locations (i,j) with
! i-ml .le. j .le. i+mu, and the full bandwidth is ml+mu+1.
!
! c. if the problem is stiff, you are encouraged to supply the jacobian
! directly (mf = 21 or 24), but if this is not feasible, dvode will
! compute it internally by difference quotients (mf = 22 or 25).
! if you are supplying the jacobian, provide a subroutine of the form:
!           subroutine jac (neq, t, y, ml, mu, pd, nrowpd, rpar, ipar)
!           double precision t, y(neq), pd(nrowpd,neq), rpar
! which supplies df/dy by loading pd as follows:
!     for a full jacobian (mf = 21), load pd(i,j) with df(i)/dy(j),
! the partial derivative of f(i) with respect to y(j).  (ignore the
! ml and mu arguments in this case.)
!     for a banded jacobian (mf = 24), load pd(i-j+mu+1,j) with
! df(i)/dy(j), i.e. load the diagonal lines of df/dy into the rows of
! pd from the top down.
!     in either case, only nonzero elements need be loaded.
!
! d. write a main program which calls subroutine dvode once for
! each point at which answers are desired.  this should also provide
! for possible use of logical unit 6 for output of error messages
! by dvode.  on the first call to dvode, supply arguments as follows:
! f      = name of subroutine for right-hand side vector f.
!          this name must be declared external in calling program.
! neq    = number of first order odes.
! y      = array of initial values, of length neq.
! t      = the initial value of the independent variable.
! tout   = first point where output is desired (.ne. t).
! itol   = 1 or 2 according as atol (below) is a scalar or array.
! rtol   = relative tolerance parameter (scalar).
! atol   = absolute tolerance parameter (scalar or array).
!          the estimated local error in y(i) will be controlled so as
!          to be roughly less (in magnitude) than
!             ewt(i) = rtol*abs(y(i)) + atol     if itol = 1, or
!             ewt(i) = rtol*abs(y(i)) + atol(i)  if itol = 2.
!          thus the local error test passes if, in each component,
!          either the absolute error is less than atol (or atol(i)),
!          or the relative error is less than rtol.
!          use rtol = 0.0 for pure absolute error control, and
!          use atol = 0.0 (or atol(i) = 0.0) for pure relative error
!          control.  caution: actual (global) errors may exceed these
!          local tolerances, so choose them conservatively.
! itask  = 1 for normal computation of output values of y at t = tout.
! istate = integer flag (input and output).  set istate = 1.
! iopt   = 0 to indicate no optional input used.
! rwork  = real work array of length at least:
!             20 + 16*neq                      for mf = 10,
!             22 +  9*neq + 2*neq**2           for mf = 21 or 22,
!             22 + 11*neq + (3*ml + 2*mu)*neq  for mf = 24 or 25.
! lrw    = declared length of rwork (in user's dimension statement).
! iwork  = integer work array of length at least:
!             30        for mf = 10,
!             30 + neq  for mf = 21, 22, 24, or 25.
!          if mf = 24 or 25, input in iwork(1),iwork(2) the lower
!          and upper half-bandwidths ml,mu.
! liw    = declared length of iwork (in user's dimension statement).
! jac    = name of subroutine for jacobian matrix (mf = 21 or 24).
!          if used, this name must be declared external in calling
!          program.  if not used, pass a dummy name.
! mf     = method flag.  standard values are:
!          10 for nonstiff (adams) method, no jacobian used.
!          21 for stiff (bdf) method, user-supplied full jacobian.
!          22 for stiff method, internally generated full jacobian.
!          24 for stiff method, user-supplied banded jacobian.
!          25 for stiff method, internally generated banded jacobian.
! rpar,ipar = user-defined real and integer arrays passed to f and jac.
! note that the main program must declare arrays y, rwork, iwork,
! and possibly atol, rpar, and ipar.
!
! e. the output from the first call (or any call) is:
!      y = array of computed values of y(t) vector.
!      t = corresponding value of independent variable (normally tout).
! istate = 2  if dvode was successful, negative otherwise.
!          -1 means excess work done on this call. (perhaps wrong mf.)
!          -2 means excess accuracy requested. (tolerances too small.)
!          -3 means illegal input detected. (see printed message.)
!          -4 means repeated error test failures. (check all input.)
!          -5 means repeated convergence failures. (perhaps bad
!             jacobian supplied or wrong choice of mf or tolerances.)
!          -6 means error weight became zero during problem. (solution
!             component i vanished, and atol or atol(i) = 0.)
!
! f. to continue the integration after a successful return, simply
! reset tout and call dvode again.  no other parameters need be reset.
!
!-----------------------------------------------------------------------
! example problem
!
! the following is a simple example problem, with the coding
! needed for its solution by dvode.  the problem is from chemical
! kinetics, and consists of the following three rate equations:
!     dy1/dt = -.04*y1 + 1.e4*y2*y3
!     dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
!     dy3/dt = 3.e7*y2**2
! on the interval from t = 0.0 to t = 4.e10, with initial conditions
! y1 = 1.0, y2 = y3 = 0.  the problem is stiff.
!
! the following coding solves this problem with dvode, using mf = 21
! and printing results at t = .4, 4., ..., 4.e10.  it uses
! itol = 2 and atol much smaller for y2 than y1 or y3 because
! y2 has much smaller values.
! at the end of the run, statistical quantities of interest are
! printed. (see optional output in the full description below.)
! to generate fortran source code, replace c in column 1 with a blank
! in the coding below.
!
!     external fex, jex
!     double precision atol, rpar, rtol, rwork, t, tout, y
!     dimension y(3), atol(3), rwork(67), iwork(33)
!     neq = 3
!     y(1) = 1.0d0
!     y(2) = 0.0d0
!     y(3) = 0.0d0
!     t = 0.0d0
!     tout = 0.4d0
!     itol = 2
!     rtol = 1.d-4
!     atol(1) = 1.d-8
!     atol(2) = 1.d-14
!     atol(3) = 1.d-6
!     itask = 1
!     istate = 1
!     iopt = 0
!     lrw = 67
!     liw = 33
!     mf = 21
!     do 40 iout = 1,12
!       call dvode(fex,neq,y,t,tout,itol,rtol,atol,itask,istate,
!    1            iopt,rwork,lrw,iwork,liw,jex,mf,rpar,ipar)
!       write(6,20)t,y(1),y(2),y(3)
! 20    format(' at t =',d12.4,'   y =',3d14.6)
!       if (istate .lt. 0) go to 80
! 40    tout = tout*10.
!     write(6,60) iwork(11),iwork(12),iwork(13),iwork(19),
!    1            iwork(20),iwork(21),iwork(22)
! 60  format(/' no. steps =',i4,'   no. f-s =',i4,
!    1       '   no. j-s =',i4,'   no. lu-s =',i4/
!    2       '  no. nonlinear iterations =',i4/
!    3       '  no. nonlinear convergence failures =',i4/
!    4       '  no. error test failures =',i4/)
!     stop
! 80  write(6,90)istate
! 90  format(///' error halt: istate =',i3)
!     stop
!     end
!
!     subroutine fex (neq, t, y, ydot, rpar, ipar)
!     double precision rpar, t, y, ydot
!     dimension y(neq), ydot(neq)
!     ydot(1) = -.04d0*y(1) + 1.d4*y(2)*y(3)
!     ydot(3) = 3.d7*y(2)*y(2)
!     ydot(2) = -ydot(1) - ydot(3)
!     return
!     end
!
!     subroutine jex (neq, t, y, ml, mu, pd, nrpd, rpar, ipar)
!     double precision pd, rpar, t, y
!     dimension y(neq), pd(nrpd,neq)
!     pd(1,1) = -.04d0
!     pd(1,2) = 1.d4*y(3)
!     pd(1,3) = 1.d4*y(2)
!     pd(2,1) = .04d0
!     pd(2,3) = -pd(1,3)
!     pd(3,2) = 6.d7*y(2)
!     pd(2,2) = -pd(1,2) - pd(3,2)
!     return
!     end
!
! the following output was obtained from the above program on a
! cray-1 computer with the cft compiler.
!
! at t =  4.0000e-01   y =  9.851680e-01  3.386314e-05  1.479817e-02
! at t =  4.0000e+00   y =  9.055255e-01  2.240539e-05  9.445214e-02
! at t =  4.0000e+01   y =  7.158108e-01  9.184883e-06  2.841800e-01
! at t =  4.0000e+02   y =  4.505032e-01  3.222940e-06  5.494936e-01
! at t =  4.0000e+03   y =  1.832053e-01  8.942690e-07  8.167938e-01
! at t =  4.0000e+04   y =  3.898560e-02  1.621875e-07  9.610142e-01
! at t =  4.0000e+05   y =  4.935882e-03  1.984013e-08  9.950641e-01
! at t =  4.0000e+06   y =  5.166183e-04  2.067528e-09  9.994834e-01
! at t =  4.0000e+07   y =  5.201214e-05  2.080593e-10  9.999480e-01
! at t =  4.0000e+08   y =  5.213149e-06  2.085271e-11  9.999948e-01
! at t =  4.0000e+09   y =  5.183495e-07  2.073399e-12  9.999995e-01
! at t =  4.0000e+10   y =  5.450996e-08  2.180399e-13  9.999999e-01
!
! no. steps = 595   no. f-s = 832   no. j-s =  13   no. lu-s = 112
!  no. nonlinear iterations = 831
!  no. nonlinear convergence failures =   0
!  no. error test failures =  22
!-----------------------------------------------------------------------
! full description of user interface to dvode.
!
! the user interface to dvode consists of the following parts.
!
! i.   the call sequence to subroutine dvode, which is a driver
!      routine for the solver.  this includes descriptions of both
!      the call sequence arguments and of user-supplied routines.
!      following these descriptions is
!        * a description of optional input available through the
!          call sequence,
!        * a description of optional output (in the work arrays), and
!        * instructions for interrupting and restarting a solution.
!
! ii.  descriptions of other routines in the dvode package that may be
!      (optionally) called by the user.  these provide the ability to
!      alter error message handling, save and restore the internal
!      common, and obtain specified derivatives of the solution y(t).
!
! iii. descriptions of common blocks to be declared in overlay
!      or similar environments.
!
! iv.  description of two routines in the dvode package, either of
!      which the user may replace with his own version, if desired.
!      these relate to the measurement of errors.
!
!-----------------------------------------------------------------------
! part i.  call sequence.
!
! the call sequence parameters used for input only are
!     f, neq, tout, itol, rtol, atol, itask, iopt, lrw, liw, jac, mf,
! and those used for both input and output are
!     y, t, istate.
! the work arrays rwork and iwork are also used for conditional and
! optional input and optional output.  (the term output here refers
! to the return from subroutine dvode to the user's calling program.)
!
! the legality of input parameters will be thoroughly checked on the
! initial call for the problem, but not checked thereafter unless a
! change in input parameters is flagged by istate = 3 in the input.
!
! the descriptions of the call arguments are as follows.
!
! f      = the name of the user-supplied subroutine defining the
!          ode system.  the system must be put in the first-order
!          form dy/dt = f(t,y), where f is a vector-valued function
!          of the scalar t and the vector y.  subroutine f is to
!          compute the function f.  it is to have the form
!               subroutine f (neq, t, y, ydot, rpar, ipar)
!               double precision t, y(neq), ydot(neq), rpar
!          where neq, t, and y are input, and the array ydot = f(t,y)
!          is output.  y and ydot are arrays of length neq.
!          subroutine f should not alter y(1),...,y(neq).
!          f must be declared external in the calling program.
!
!          subroutine f may access user-defined real and integer
!          work arrays rpar and ipar, which are to be dimensioned
!          in the main program.
!
!          if quantities computed in the f routine are needed
!          externally to dvode, an extra call to f should be made
!          for this purpose, for consistent and accurate results.
!          if only the derivative dy/dt is needed, use dvindy instead.
!
! neq    = the size of the ode system (number of first order
!          ordinary differential equations).  used only for input.
!          neq may not be increased during the problem, but
!          can be decreased (with istate = 3 in the input).
!
! y      = a real array for the vector of dependent variables, of
!          length neq or more.  used for both input and output on the
!          first call (istate = 1), and only for output on other calls.
!          on the first call, y must contain the vector of initial
!          values.  in the output, y contains the computed solution
!          evaluated at t.  if desired, the y array may be used
!          for other purposes between calls to the solver.
!
!          this array is passed as the y argument in all calls to
!          f and jac.
!
! t      = the independent variable.  in the input, t is used only on
!          the first call, as the initial point of the integration.
!          in the output, after each call, t is the value at which a
!          computed solution y is evaluated (usually the same as tout).
!          on an error return, t is the farthest point reached.
!
! tout   = the next value of t at which a computed solution is desired.
!          used only for input.
!
!          when starting the problem (istate = 1), tout may be equal
!          to t for one call, then should .ne. t for the next call.
!          for the initial t, an input value of tout .ne. t is used
!          in order to determine the direction of the integration
!          (i.e. the algebraic sign of the step sizes) and the rough
!          scale of the problem.  integration in either direction
!          (forward or backward in t) is permitted.
!
!          if itask = 2 or 5 (one-step modes), tout is ignored after
!          the first call (i.e. the first call with tout .ne. t).
!          otherwise, tout is required on every call.
!
!          if itask = 1, 3, or 4, the values of tout need not be
!          monotone, but a value of tout which backs up is limited
!          to the current internal t interval, whose endpoints are
!          tcur - hu and tcur.  (see optional output, below, for
!          tcur and hu.)
!
! itol   = an indicator for the type of error control.  see
!          description below under atol.  used only for input.
!
! rtol   = a relative error tolerance parameter, either a scalar or
!          an array of length neq.  see description below under atol.
!          input only.
!
! atol   = an absolute error tolerance parameter, either a scalar or
!          an array of length neq.  input only.
!
!          the input parameters itol, rtol, and atol determine
!          the error control performed by the solver.  the solver will
!          control the vector e = (e(i)) of estimated local errors
!          in y, according to an inequality of the form
!                      rms-norm of ( e(i)/ewt(i) )   .le.   1,
!          where       ewt(i) = rtol(i)*abs(y(i)) + atol(i),
!          and the rms-norm (root-mean-square norm) here is
!          rms-norm(v) = sqrt(sum v(i)**2 / neq).  here ewt = (ewt(i))
!          is a vector of weights which must always be positive, and
!          the values of rtol and atol should all be non-negative.
!          the following table gives the types (scalar/array) of
!          rtol and atol, and the corresponding form of ewt(i).
!
!             itol    rtol       atol          ewt(i)
!              1     scalar     scalar     rtol*abs(y(i)) + atol
!              2     scalar     array      rtol*abs(y(i)) + atol(i)
!              3     array      scalar     rtol(i)*abs(y(i)) + atol
!              4     array      array      rtol(i)*abs(y(i)) + atol(i)
!
!          when either of these parameters is a scalar, it need not
!          be dimensioned in the user's calling program.
!
!          if none of the above choices (with itol, rtol, and atol
!          fixed throughout the problem) is suitable, more general
!          error controls can be obtained by substituting
!          user-supplied routines for the setting of ewt and/or for
!          the norm calculation.  see part iv below.
!
!          if global errors are to be estimated by making a repeated
!          run on the same problem with smaller tolerances, then all
!          components of rtol and atol (i.e. of ewt) should be scaled
!          down uniformly.
!
! itask  = an index specifying the task to be performed.
!          input only.  itask has the following values and meanings.
!          1  means normal computation of output values of y(t) at
!             t = tout (by overshooting and interpolating).
!          2  means take one step only and return.
!          3  means stop at the first internal mesh point at or
!             beyond t = tout and return.
!          4  means normal computation of output values of y(t) at
!             t = tout but without overshooting t = tcrit.
!             tcrit must be input as rwork(1).  tcrit may be equal to
!             or beyond tout, but not behind it in the direction of
!             integration.  this option is useful if the problem
!             has a singularity at or beyond t = tcrit.
!          5  means take one step, without passing tcrit, and return.
!             tcrit must be input as rwork(1).
!
!          note:  if itask = 4 or 5 and the solver reaches tcrit
!          (within roundoff), it will return t = tcrit (exactly) to
!          indicate this (unless itask = 4 and tout comes before tcrit,
!          in which case answers at t = tout are returned first).
!
! istate = an index used for input and output to specify the
!          the state of the calculation.
!
!          in the input, the values of istate are as follows.
!          1  means this is the first call for the problem
!             (initializations will be done).  see note below.
!          2  means this is not the first call, and the calculation
!             is to continue normally, with no change in any input
!             parameters except possibly tout and itask.
!             (if itol, rtol, and/or atol are changed between calls
!             with istate = 2, the new values will be used but not
!             tested for legality.)
!          3  means this is not the first call, and the
!             calculation is to continue normally, but with
!             a change in input parameters other than
!             tout and itask.  changes are allowed in
!             neq, itol, rtol, atol, iopt, lrw, liw, mf, ml, mu,
!             and any of the optional input except h0.
!             (see iwork description for ml and mu.)
!          note:  a preliminary call with tout = t is not counted
!          as a first call here, as no initialization or checking of
!          input is done.  (such a call is sometimes useful to include
!          the initial conditions in the output.)
!          thus the first call for which tout .ne. t requires
!          istate = 1 in the input.
!
!          in the output, istate has the following values and meanings.
!           1  means nothing was done, as tout was equal to t with
!              istate = 1 in the input.
!           2  means the integration was performed successfully.
!          -1  means an excessive amount of work (more than mxstep
!              steps) was done on this call, before completing the
!              requested task, but the integration was otherwise
!              successful as far as t.  (mxstep is an optional input
!              and is normally 500.)  to continue, the user may
!              simply reset istate to a value .gt. 1 and call again.
!              (the excess work step counter will be reset to 0.)
!              in addition, the user may increase mxstep to avoid
!              this error return.  (see optional input below.)
!          -2  means too much accuracy was requested for the precision
!              of the machine being used.  this was detected before
!              completing the requested task, but the integration
!              was successful as far as t.  to continue, the tolerance
!              parameters must be reset, and istate must be set
!              to 3.  the optional output tolsf may be used for this
!              purpose.  (note: if this condition is detected before
!              taking any steps, then an illegal input return
!              (istate = -3) occurs instead.)
!          -3  means illegal input was detected, before taking any
!              integration steps.  see written message for details.
!              note:  if the solver detects an infinite loop of calls
!              to the solver with illegal input, it will cause
!              the run to stop.
!          -4  means there were repeated error test failures on
!              one attempted step, before completing the requested
!              task, but the integration was successful as far as t.
!              the problem may have a singularity, or the input
!              may be inappropriate.
!          -5  means there were repeated convergence test failures on
!              one attempted step, before completing the requested
!              task, but the integration was successful as far as t.
!              this may be caused by an inaccurate jacobian matrix,
!              if one is being used.
!          -6  means ewt(i) became zero for some i during the
!              integration.  pure relative error control (atol(i)=0.0)
!              was requested on a variable which has now vanished.
!              the integration was successful as far as t.
!
!          note:  since the normal output value of istate is 2,
!          it does not need to be reset for normal continuation.
!          also, since a negative input value of istate will be
!          regarded as illegal, a negative output value requires the
!          user to change it, and possibly other input, before
!          calling the solver again.
!
! iopt   = an integer flag to specify whether or not any optional
!          input is being used on this call.  input only.
!          the optional input is listed separately below.
!          iopt = 0 means no optional input is being used.
!                   default values will be used in all cases.
!          iopt = 1 means optional input is being used.
!
! rwork  = a real working array (double precision).
!          the length of rwork must be at least
!             20 + nyh*(maxord + 1) + 3*neq + lwm    where
!          nyh    = the initial value of neq,
!          maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
!                   smaller value is given as an optional input),
!          lwm = length of work space for matrix-related data:
!          lwm = 0             if miter = 0,
!          lwm = 2*neq**2 + 2  if miter = 1 or 2, and mf.gt.0,
!          lwm = neq**2 + 2    if miter = 1 or 2, and mf.lt.0,
!          lwm = neq + 2       if miter = 3,
!          lwm = (3*ml+2*mu+2)*neq + 2 if miter = 4 or 5, and mf.gt.0,
!          lwm = (2*ml+mu+1)*neq + 2   if miter = 4 or 5, and mf.lt.0.
!          (see the mf description for meth and miter.)
!          thus if maxord has its default value and neq is constant,
!          this length is:
!             20 + 16*neq                    for mf = 10,
!             22 + 16*neq + 2*neq**2         for mf = 11 or 12,
!             22 + 16*neq + neq**2           for mf = -11 or -12,
!             22 + 17*neq                    for mf = 13,
!             22 + 18*neq + (3*ml+2*mu)*neq  for mf = 14 or 15,
!             22 + 17*neq + (2*ml+mu)*neq    for mf = -14 or -15,
!             20 +  9*neq                    for mf = 20,
!             22 +  9*neq + 2*neq**2         for mf = 21 or 22,
!             22 +  9*neq + neq**2           for mf = -21 or -22,
!             22 + 10*neq                    for mf = 23,
!             22 + 11*neq + (3*ml+2*mu)*neq  for mf = 24 or 25.
!             22 + 10*neq + (2*ml+mu)*neq    for mf = -24 or -25.
!          the first 20 words of rwork are reserved for conditional
!          and optional input and optional output.
!
!          the following word in rwork is a conditional input:
!            rwork(1) = tcrit = critical value of t which the solver
!                       is not to overshoot.  required if itask is
!                       4 or 5, and ignored otherwise.  (see itask.)
!
! lrw    = the length of the array rwork, as declared by the user.
!          (this will be checked by the solver.)
!
! iwork  = an integer work array.  the length of iwork must be at least
!             30        if miter = 0 or 3 (mf = 10, 13, 20, 23), or
!             30 + neq  otherwise (abs(mf) = 11,12,14,15,21,22,24,25).
!          the first 30 words of iwork are reserved for conditional and
!          optional input and optional output.
!
!          the following 2 words in iwork are conditional input:
!            iwork(1) = ml     these are the lower and upper
!            iwork(2) = mu     half-bandwidths, respectively, of the
!                       banded jacobian, excluding the main diagonal.
!                       the band is defined by the matrix locations
!                       (i,j) with i-ml .le. j .le. i+mu.  ml and mu
!                       must satisfy  0 .le.  ml,mu  .le. neq-1.
!                       these are required if miter is 4 or 5, and
!                       ignored otherwise.  ml and mu may in fact be
!                       the band parameters for a matrix to which
!                       df/dy is only approximately equal.
!
! liw    = the length of the array iwork, as declared by the user.
!          (this will be checked by the solver.)
!
! note:  the work arrays must not be altered between calls to dvode
! for the same problem, except possibly for the conditional and
! optional input, and except for the last 3*neq words of rwork.
! the latter space is used for internal scratch space, and so is
! available for use by the user outside dvode between calls, if
! desired (but not for use by f or jac).
!
! jac    = the name of the user-supplied routine (miter = 1 or 4) to
!          compute the jacobian matrix, df/dy, as a function of
!          the scalar t and the vector y.  it is to have the form
!               subroutine jac (neq, t, y, ml, mu, pd, nrowpd,
!                               rpar, ipar)
!               double precision t, y(neq), pd(nrowpd,neq), rpar
!          where neq, t, y, ml, mu, and nrowpd are input and the array
!          pd is to be loaded with partial derivatives (elements of the
!          jacobian matrix) in the output.  pd must be given a first
!          dimension of nrowpd.  t and y have the same meaning as in
!          subroutine f.
!               in the full matrix case (miter = 1), ml and mu are
!          ignored, and the jacobian is to be loaded into pd in
!          columnwise manner, with df(i)/dy(j) loaded into pd(i,j).
!               in the band matrix case (miter = 4), the elements
!          within the band are to be loaded into pd in columnwise
!          manner, with diagonal lines of df/dy loaded into the rows
!          of pd. thus df(i)/dy(j) is to be loaded into pd(i-j+mu+1,j).
!          ml and mu are the half-bandwidth parameters. (see iwork).
!          the locations in pd in the two triangular areas which
!          correspond to nonexistent matrix elements can be ignored
!          or loaded arbitrarily, as they are overwritten by dvode.
!               jac need not provide df/dy exactly.  a crude
!          approximation (possibly with a smaller bandwidth) will do.
!               in either case, pd is preset to zero by the solver,
!          so that only the nonzero elements need be loaded by jac.
!          each call to jac is preceded by a call to f with the same
!          arguments neq, t, and y.  thus to gain some efficiency,
!          intermediate quantities shared by both calculations may be
!          saved in a user common block by f and not recomputed by jac,
!          if desired.  also, jac may alter the y array, if desired.
!          jac must be declared external in the calling program.
!               subroutine jac may access user-defined real and integer
!          work arrays, rpar and ipar, whose dimensions are set by the
!          user in the main program.
!
! mf     = the method flag.  used only for input.  the legal values of
!          mf are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, 25,
!          -11, -12, -14, -15, -21, -22, -24, -25.
!          mf is a signed two-digit integer, mf = jsv*(10*meth + miter).
!          jsv = sign(mf) indicates the jacobian-saving strategy:
!            jsv =  1 means a copy of the jacobian is saved for reuse
!                     in the corrector iteration algorithm.
!            jsv = -1 means a copy of the jacobian is not saved
!                     (valid only for miter = 1, 2, 4, or 5).
!          meth indicates the basic linear multistep method:
!            meth = 1 means the implicit adams method.
!            meth = 2 means the method based on backward
!                     differentiation formulas (bdf-s).
!          miter indicates the corrector iteration method:
!            miter = 0 means functional iteration (no jacobian matrix
!                      is involved).
!            miter = 1 means chord iteration with a user-supplied
!                      full (neq by neq) jacobian.
!            miter = 2 means chord iteration with an internally
!                      generated (difference quotient) full jacobian
!                      (using neq extra calls to f per df/dy value).
!            miter = 3 means chord iteration with an internally
!                      generated diagonal jacobian approximation
!                      (using 1 extra call to f per df/dy evaluation).
!            miter = 4 means chord iteration with a user-supplied
!                      banded jacobian.
!            miter = 5 means chord iteration with an internally
!                      generated banded jacobian (using ml+mu+1 extra
!                      calls to f per df/dy evaluation).
!          if miter = 1 or 4, the user must supply a subroutine jac
!          (the name is arbitrary) as described above under jac.
!          for other values of miter, a dummy argument can be used.
!
! rpar     user-specified array used to communicate real parameters
!          to user-supplied subroutines.  if rpar is a vector, then
!          it must be dimensioned in the user's main program.  if it
!          is unused or it is a scalar, then it need not be
!          dimensioned.
!
! ipar     user-specified array used to communicate integer parameter
!          to user-supplied subroutines.  the comments on dimensioning
!          rpar apply to ipar.
!-----------------------------------------------------------------------
! optional input.
!
! the following is a list of the optional input provided for in the
! call sequence.  (see also part ii.)  for each such input variable,
! this table lists its name as used in this documentation, its
! location in the call sequence, its meaning, and the default value.
! the use of any of this input requires iopt = 1, and in that
! case all of this input is examined.  a value of zero for any
! of these optional input variables will cause the default value to be
! used.  thus to use a subset of the optional input, simply preload
! locations 5 to 10 in rwork and iwork to 0.0 and 0 respectively, and
! then set those of interest to nonzero values.
!
! name    location      meaning and default value
!
! h0      rwork(5)  the step size to be attempted on the first step.
!                   the default value is determined by the solver.
!
! hmax    rwork(6)  the maximum absolute step size allowed.
!                   the default value is infinite.
!
! hmin    rwork(7)  the minimum absolute step size allowed.
!                   the default value is 0.  (this lower bound is not
!                   enforced on the final step before reaching tcrit
!                   when itask = 4 or 5.)
!
! maxord  iwork(5)  the maximum order to be allowed.  the default
!                   value is 12 if meth = 1, and 5 if meth = 2.
!                   if maxord exceeds the default value, it will
!                   be reduced to the default value.
!                   if maxord is changed during the problem, it may
!                   cause the current order to be reduced.
!
! mxstep  iwork(6)  maximum number of (internally defined) steps
!                   allowed during one call to the solver.
!                   the default value is 500.
!
! mxhnil  iwork(7)  maximum number of messages printed (per problem)
!                   warning that t + h = t on a step (h = step size).
!                   this must be positive to result in a non-default
!                   value.  the default value is 10.
!
!-----------------------------------------------------------------------
! optional output.
!
! as optional additional output from dvode, the variables listed
! below are quantities related to the performance of dvode
! which are available to the user.  these are communicated by way of
! the work arrays, but also have internal mnemonic names as shown.
! except where stated otherwise, all of this output is defined
! on any successful return from dvode, and on any return with
! istate = -1, -2, -4, -5, or -6.  on an illegal input return
! (istate = -3), they will be unchanged from their existing values
! (if any), except possibly for tolsf, lenrw, and leniw.
! on any error return, output relevant to the error will be defined,
! as noted below.
!
! name    location      meaning
!
! hu      rwork(11) the step size in t last used (successfully).
!
! hcur    rwork(12) the step size to be attempted on the next step.
!
! tcur    rwork(13) the current value of the independent variable
!                   which the solver has actually reached, i.e. the
!                   current internal mesh point in t.  in the output,
!                   tcur will always be at least as far from the
!                   initial value of t as the current argument t,
!                   but may be farther (if interpolation was done).
!
! tolsf   rwork(14) a tolerance scale factor, greater than 1.0,
!                   computed when a request for too much accuracy was
!                   detected (istate = -3 if detected at the start of
!                   the problem, istate = -2 otherwise).  if itol is
!                   left unaltered but rtol and atol are uniformly
!                   scaled up by a factor of tolsf for the next call,
!                   then the solver is deemed likely to succeed.
!                   (the user may also ignore tolsf and alter the
!                   tolerance parameters in any other way appropriate.)
!
! nst     iwork(11) the number of steps taken for the problem so far.
!
! nfe     iwork(12) the number of f evaluations for the problem so far.
!
! nje     iwork(13) the number of jacobian evaluations so far.
!
! nqu     iwork(14) the method order last used (successfully).
!
! nqcur   iwork(15) the order to be attempted on the next step.
!
! imxer   iwork(16) the index of the component of largest magnitude in
!                   the weighted local error vector ( e(i)/ewt(i) ),
!                   on an error return with istate = -4 or -5.
!
! lenrw   iwork(17) the length of rwork actually required.
!                   this is defined on normal returns and on an illegal
!                   input return for insufficient storage.
!
! leniw   iwork(18) the length of iwork actually required.
!                   this is defined on normal returns and on an illegal
!                   input return for insufficient storage.
!
! nlu     iwork(19) the number of matrix lu decompositions so far.
!
! nni     iwork(20) the number of nonlinear (newton) iterations so far.
!
! ncfn    iwork(21) the number of convergence failures of the nonlinear
!                   solver so far.
!
! netf    iwork(22) the number of error test failures of the integrator
!                   so far.
!
! the following two arrays are segments of the rwork array which
! may also be of interest to the user as optional output.
! for each array, the table below gives its internal name,
! its base address in rwork, and its description.
!
! name    base address      description
!
! yh      21             the nordsieck history array, of size nyh by
!                        (nqcur + 1), where nyh is the initial value
!                        of neq.  for j = 0,1,...,nqcur, column j+1
!                        of yh contains hcur**j/factorial(j) times
!                        the j-th derivative of the interpolating
!                        polynomial currently representing the
!                        solution, evaluated at t = tcur.
!
! acor     lenrw-neq+1   array of size neq used for the accumulated
!                        corrections on each step, scaled in the output
!                        to represent the estimated local error in y
!                        on the last step.  this is the vector e in
!                        the description of the error control.  it is
!                        defined only on a successful return from dvode.
!
!-----------------------------------------------------------------------
! interrupting and restarting
!
! if the integration of a given problem by dvode is to be
! interrrupted and then later continued, such as when restarting
! an interrupted run or alternating between two or more ode problems,
! the user should save, following the return from the last dvode call
! prior to the interruption, the contents of the call sequence
! variables and internal common blocks, and later restore these
! values before the next dvode call for that problem.  to save
! and restore the common blocks, use subroutine dvsrco, as
! described below in part ii.
!
! in addition, if non-default values for either lun or mflag are
! desired, an extra call to xsetun and/or xsetf should be made just
! before continuing the integration.  see part ii below for details.
!
!-----------------------------------------------------------------------
! part ii.  other routines callable.
!
! the following are optional calls which the user may make to
! gain additional capabilities in conjunction with dvode.
! (the routines xsetun and xsetf are designed to conform to the
! slatec error handling package.)
!
!     form of call                  function
!  call xsetun(lun)           set the logical unit number, lun, for
!                             output of messages from dvode, if
!                             the default is not desired.
!                             the default value of lun is 6.
!
!  call xsetf(mflag)          set a flag to control the printing of
!                             messages by dvode.
!                             mflag = 0 means do not print. (danger:
!                             this risks losing valuable information.)
!                             mflag = 1 means print (the default).
!
!                             either of the above calls may be made at
!                             any time and will take effect immediately.
!
!  call dvsrco(rsav,isav,job) saves and restores the contents of
!                             the internal common blocks used by
!                             dvode. (see part iii below.)
!                             rsav must be a real array of length 49
!                             or more, and isav must be an integer
!                             array of length 40 or more.
!                             job=1 means save common into rsav/isav.
!                             job=2 means restore common from rsav/isav.
!                                dvsrco is useful if one is
!                             interrupting a run and restarting
!                             later, or alternating between two or
!                             more problems solved with dvode.
!
!  call dvindy(,,,,,)         provide derivatives of y, of various
!        (see below.)         orders, at a specified point t, if
!                             desired.  it may be called only after
!                             a successful return from dvode.
!
! the detailed instructions for using dvindy are as follows.
! the form of the call is:
!
!  call dvindy (t, k, rwork(21), nyh, dky, iflag)
!
! the input parameters are:
!
! t         = value of independent variable where answers are desired
!             (normally the same as the t last returned by dvode).
!             for valid results, t must lie between tcur - hu and tcur.
!             (see optional output for tcur and hu.)
! k         = integer order of the derivative desired.  k must satisfy
!             0 .le. k .le. nqcur, where nqcur is the current order
!             (see optional output).  the capability corresponding
!             to k = 0, i.e. computing y(t), is already provided
!             by dvode directly.  since nqcur .ge. 1, the first
!             derivative dy/dt is always available with dvindy.
! rwork(21) = the base address of the history array yh.
! nyh       = column length of yh, equal to the initial value of neq.
!
! the output parameters are:
!
! dky       = a real array of length neq containing the computed value
!             of the k-th derivative of y(t).
! iflag     = integer flag, returned as 0 if k and t were legal,
!             -1 if k was illegal, and -2 if t was illegal.
!             on an error return, a message is also written.
!-----------------------------------------------------------------------
! part iii.  common blocks.
! if dvode is to be used in an overlay situation, the user
! must declare, in the primary overlay, the variables in:
!   (1) the call sequence to dvode,
!   (2) the two internal common blocks
!         /dvod_cmn01/  of length  81  (48 double precision words
!                         followed by 33 integer words),
!         /dvod_cmn02/  of length  9  (1 double precision word
!                         followed by 8 integer words),
!
! if dvode is used on a system in which the contents of internal
! common blocks are not preserved between calls, the user should
! declare the above two common blocks in his main program to insure
! that their contents are preserved.
!
!-----------------------------------------------------------------------
! part iv.  optionally replaceable solver routines.
!
! below are descriptions of two routines in the dvode package which
! relate to the measurement of errors.  either routine can be
! replaced by a user-supplied version, if desired.  however, since such
! a replacement may have a major impact on performance, it should be
! done only when absolutely necessary, and only with great caution.
! (note: the means by which the package version of a routine is
! superseded by the user's version may be system-dependent.)
!
! (a) dewset.
! the following subroutine is called just before each internal
! integration step, and sets the array of error weights, ewt, as
! described under itol/rtol/atol above:
!     subroutine dewset (neq, itol, rtol, atol, ycur, ewt)
! where neq, itol, rtol, and atol are as in the dvode call sequence,
! ycur contains the current dependent variable vector, and
! ewt is the array of weights set by dewset.
!
! if the user supplies this subroutine, it must return in ewt(i)
! (i = 1,...,neq) a positive quantity suitable for comparison with
! errors in y(i).  the ewt array returned by dewset is passed to the
! dvnorm routine (see below.), and also used by dvode in the computation
! of the optional output imxer, the diagonal jacobian approximation,
! and the increments for difference quotient jacobians.
!
! in the user-supplied version of dewset, it may be desirable to use
! the current values of derivatives of y.  derivatives up to order nq
! are available from the history array yh, described above under
! optional output.  in dewset, yh is identical to the ycur array,
! extended to nq + 1 columns with a column length of nyh and scale
! factors of h**j/factorial(j).  on the first call for the problem,
! given by nst = 0, nq is 1 and h is temporarily set to 1.0.
! nyh is the initial value of neq.  the quantities nq, h, and nst
! can be obtained by including in dewset the statements:
!     double precision rvod, h, hu
!     common /dvod_cmn01/ rvod(48), ivod(33)
!     common /dvod_cmn02/ hu, ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!     nq = ivod(28)
!     h = rvod(21)
! thus, for example, the current value of dy/dt can be obtained as
! ycur(nyh+i)/h  (i=1,...,neq)  (and the division by h is
! unnecessary when nst = 0).
!
! (b) dvnorm.
! the following is a real function routine which computes the weighted
! root-mean-square norm of a vector v:
!     d = dvnorm (n, v, w)
! where:
!   n = the length of the vector,
!   v = real array of length n containing the vector,
!   w = real array of length n containing weights,
!   d = sqrt( (1/n) * sum(v(i)*w(i))**2 ).
! dvnorm is called with n = neq and with w(i) = 1.0/ewt(i), where
! ewt is as set by subroutine dewset.
!
! if the user supplies this function, it should return a non-negative
! value of dvnorm suitable for use in the error control in dvode.
! none of the arguments should be altered by dvnorm.
! for example, a user-supplied dvnorm routine might:
!   -substitute a max-norm of (v(i)*w(i)) for the rms-norm, or
!   -ignore some components of v in the norm, with the effect of
!    suppressing the error control on those components of y.
!-----------------------------------------------------------------------
! revision history (yyyymmdd)
!  19890615  date written.  initial release.
!  19890922  added interrupt/restart ability, minor changes throughout.
!  19910228  minor revisions in line format,  prologue, etc.
!  19920227  modifications by d. pang:
!            (1) applied subgennam to get generic intrinsic names.
!            (2) changed intrinsic names to generic in comments.
!            (3) added *deck lines before each routine.
!  19920721  names of routines and labeled common blocks changed, so as
!            to be unique in combined single/double precision code (ach).
!  19920722  minor revisions to prologue (ach).
!  19920831  conversion to double precision done (ach).
!  19921106  fixed minor bug: etaq,etaqm1 in dvstep save statement (ach).
!  19921118  changed lunsav/mflgsv to ixsav (ach).
!  19941222  removed mf overwrite; attached sign to h in estimated second
!            deriv. in dvhin; misc. comment changes throughout (ach).
!  19970515  minor corrections to comments in prologue, dvjac (ach).
!  19981111  corrected block b by adding final line, go to 200 (ach).
!  20020430  various upgrades (ach): use odepack error handler package.
!            replaced d1mach by dumach.  various changes to main
!            prologue and other routine prologues.
!-----------------------------------------------------------------------
! other routines in the dvode package.
!
! in addition to subroutine dvode, the dvode package includes the
! following subroutines and function routines:
!  dvhin     computes an approximate step size for the initial step.
!  dvindy    computes an interpolated value of the y vector at t = tout.
!  dvstep    is the core integrator, which does one step of the
!            integration and the associated error control.
!  dvset     sets all method coefficients and test constants.
!  dvnlsd    solves the underlying nonlinear system -- the corrector.
!  dvjac     computes and preprocesses the jacobian matrix j = df/dy
!            and the newton iteration matrix p = i - (h/l1)*j.
!  dvsol     manages solution of linear system in chord iteration.
!  dvjust    adjusts the history array on a change of order.
!  dewset    sets the error weight vector ewt before each step.
!  dvnorm    computes the weighted r.m.s. norm of a vector.
!  dvsrco    is a user-callable routine to save and restore
!            the contents of the internal common blocks.
!  dacopy    is a routine to copy one two-dimensional array to another.
!  dgefa and dgesl   are routines from linpack for solving full
!            systems of linear algebraic equations.
!  dgbfa and dgbsl   are routines from linpack for solving banded
!            linear systems.
!  daxpy, dscal, and dcopy are basic linear algebra modules (blas).
!  dumach    sets the unit roundoff of the machine.
!  xerrwd, xsetun, xsetf, ixsav, and iumach handle the printing of all
!            error messages and warnings.  xerrwd is machine-dependent.
! note:  dvnorm, dumach, ixsav, and iumach are function routines.
! all the others are subroutines.
!
!-----------------------------------------------------------------------
!!
!! type declarations for labeled common block dvod_cmn01 --------------------
!!
!      double precision acnrm, ccmxj, conp, crate, drc, el,   &
!           eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!           rc, rl1, tau, tq, tn, uround
!      integer icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!              l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!              locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!              n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!              nslp, nyh
!!
!! type declarations for labeled common block dvod_cmn02 --------------------
!!
!      double precision hu
!      integer ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
! type declarations for local variables --------------------------------
!
! 27-oct-2005 rce - do not declare functions that are in the module
!     external dvnlsd
      logical ihit
      double precision atoli, big, ewti, four, h0, hmax, hmx, hun, one,   &
         pt2, rh, rtoli, size, tcrit, tnext, tolsf, tp, two, zero
      integer i, ier, iflag, imxer, jco, kgo, leniw, lenj, lenp, lenrw,   &
         lenwm, lf0, mband, mfa, ml, mord, mu, mxhnl0, mxstp0, niter,   &
         nslast
      character*80 msg
!
! type declaration for function subroutines called ---------------------
!
! 27-oct-2005 rce - do not declare functions that are in the module
!     double precision dumach, dvnorm
!
      dimension mord(2)
!-----------------------------------------------------------------------
! the following fortran-77 declaration is to cause the values of the
! listed (local) variables to be saved between calls to dvode.
!-----------------------------------------------------------------------
      save mord, mxhnl0, mxstp0
      save zero, one, two, four, pt2, hun
!-----------------------------------------------------------------------
! the following internal common blocks contain variables which are
! communicated between subroutines in the dvode package, or which are
! to be saved between calls to dvode.
! in each block, real variables precede integers.
! the block /dvod_cmn01/ appears in subroutines dvode, dvindy, dvstep,
! dvset, dvnlsd, dvjac, dvsol, dvjust and dvsrco.
! the block /dvod_cmn02/ appears in subroutines dvode, dvindy, dvstep,
! dvnlsd, dvjac, and dvsrco.
!
! the variables stored in the internal common blocks are as follows:
!
! acnrm  = weighted r.m.s. norm of accumulated correction vectors.
! ccmxj  = threshhold on drc for updating the jacobian. (see drc.)
! conp   = the saved value of tq(5).
! crate  = estimated corrector convergence rate constant.
! drc    = relative change in h*rl1 since last dvjac call.
! el     = real array of integration coefficients.  see dvset.
! eta    = saved tentative ratio of new to old h.
! etamax = saved maximum value of eta to be allowed.
! h      = the step size.
! hmin   = the minimum absolute value of the step size h to be used.
! hmxi   = inverse of the maximum absolute value of h to be used.
!          hmxi = 0.0 is allowed and corresponds to an infinite hmax.
! hnew   = the step size to be attempted on the next step.
! hscal  = stepsize in scaling of yh array.
! prl1   = the saved value of rl1.
! rc     = ratio of current h*rl1 to value on last dvjac call.
! rl1    = the reciprocal of the coefficient el(1).
! tau    = real vector of past nq step sizes, length 13.
! tq     = a real vector of length 5 in which dvset stores constants
!          used for the convergence test, the error test, and the
!          selection of h at a new order.
! tn     = the independent variable, updated on each step taken.
! uround = the machine unit roundoff.  the smallest positive real number
!          such that  1.0 + uround .ne. 1.0
! icf    = integer flag for convergence failure in dvnlsd:
!            0 means no failures.
!            1 means convergence failure with out of date jacobian
!                   (recoverable error).
!            2 means convergence failure with current jacobian or
!                   singular matrix (unrecoverable error).
! init   = saved integer flag indicating whether initialization of the
!          problem has been done (init = 1) or not.
! ipup   = saved flag to signal updating of newton matrix.
! jcur   = output flag from dvjac showing jacobian status:
!            jcur = 0 means j is not current.
!            jcur = 1 means j is current.
! jstart = integer flag used as input to dvstep:
!            0  means perform the first step.
!            1  means take a new step continuing from the last.
!            -1 means take the next step with a new value of maxord,
!                  hmin, hmxi, n, meth, miter, and/or matrix parameters.
!          on return, dvstep sets jstart = 1.
! jsv    = integer flag for jacobian saving, = sign(mf).
! kflag  = a completion code from dvstep with the following meanings:
!               0      the step was succesful.
!              -1      the requested error could not be achieved.
!              -2      corrector convergence could not be achieved.
!              -3, -4  fatal error in vnls (can not occur here).
! kuth   = input flag to dvstep showing whether h was reduced by the
!          driver.  kuth = 1 if h was reduced, = 0 otherwise.
! l      = integer variable, nq + 1, current order plus one.
! lmax   = maxord + 1 (used for dimensioning).
! locjs  = a pointer to the saved jacobian, whose storage starts at
!          wm(locjs), if jsv = 1.
! lyh, lewt, lacor, lsavf, lwm, liwm = saved integer pointers
!          to segments of rwork and iwork.
! maxord = the maximum order of integration method to be allowed.
! meth/miter = the method flags.  see mf.
! msbj   = the maximum number of steps between j evaluations, = 50.
! mxhnil = saved value of optional input mxhnil.
! mxstep = saved value of optional input mxstep.
! n      = the number of first-order odes, = neq.
! newh   = saved integer to flag change of h.
! newq   = the method order to be used on the next step.
! nhnil  = saved counter for occurrences of t + h = t.
! nq     = integer variable, the current integration method order.
! nqnyh  = saved value of nq*nyh.
! nqwait = a counter controlling the frequency of order changes.
!          an order change is about to be considered if nqwait = 1.
! nslj   = the number of steps taken as of the last jacobian update.
! nslp   = saved value of nst as of last newton matrix update.
! nyh    = saved value of the initial value of neq.
! hu     = the step size in t last used.
! ncfn   = number of nonlinear convergence failures so far.
! netf   = the number of error test failures of the integrator so far.
! nfe    = the number of f evaluations for the problem so far.
! nje    = the number of jacobian evaluations so far.
! nlu    = the number of matrix lu decompositions so far.
! nni    = number of nonlinear iterations so far.
! nqu    = the method order last used.
! nst    = the number of steps taken for the problem so far.
!-----------------------------------------------------------------------
!      common /dvod_cmn01/ acnrm, ccmxj, conp, crate, drc, el(13),   &
!                      eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!                      rc, rl1, tau(13), tq(5), tn, uround,   &
!                      icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!                      l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!                      locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!                      n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!                      nslp, nyh
!      common /dvod_cmn02/ hu, ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
! in this routine only, declare h and n as local variables
!    before/after calls where cmn12 is passed as an argument,
!    load/unload them to/from cmn12
      type( dvode_cmn_vars ) :: cmn12
      double precision h
      integer n
!
!
      data  mord(1) /12/, mord(2) /5/, mxstp0 /500/, mxhnl0 /10/
      data zero /0.0d0/, one /1.0d0/, two /2.0d0/, four /4.0d0/,   &
           pt2 /0.2d0/, hun /100.0d0/
!-----------------------------------------------------------------------
! block a.
! this code block is executed on every call.
! it tests istate and itask for legality and branches appropriately.
! if istate .gt. 1 but the flag init shows that initialization has
! not yet been done, an error return occurs.
! if istate = 1 and tout = t, return immediately.
!-----------------------------------------------------------------------

! in this modified-for-wrf version of dvode, istate must be 1
      if (istate .ne. 1) then
         msg = '*** module_cmu_dvode_solver/dvode -- ' // &
               'fatal error, istate must be 1'
         call wrf_error_fatal( msg )
      end if

      if (istate .lt. 1 .or. istate .gt. 3) go to 601
      if (itask .lt. 1 .or. itask .gt. 5) go to 602
      if (istate .eq. 1) go to 10
      if (cmn12%init .ne. 1) go to 603
      if (istate .eq. 2) go to 200
      go to 20
 10   cmn12%init = 0
      if (tout .eq. t) return
!-----------------------------------------------------------------------
! block b.
! the next code block is executed for the initial call (istate = 1),
! or for a continuation call with parameter changes (istate = 3).
! it contains checking of all input and various initializations.
!
! first check legality of the non-optional input neq, itol, iopt,
! mf, ml, and mu.
!-----------------------------------------------------------------------
 20   if (neq .le. 0) go to 604
      if (istate .eq. 1) go to 25
      if (neq .gt. n) go to 605
 25   n = neq
      if (itol .lt. 1 .or. itol .gt. 4) go to 606
      if (iopt .lt. 0 .or. iopt .gt. 1) go to 607
      cmn12%jsv = sign(1,mf)
      mfa = abs(mf)
      cmn12%meth = mfa/10
      cmn12%miter = mfa - 10*cmn12%meth
      if (cmn12%meth .lt. 1 .or. cmn12%meth .gt. 2) go to 608
      if (cmn12%miter .lt. 0 .or. cmn12%miter .gt. 5) go to 608
      if (cmn12%miter .le. 3) go to 30
      ml = iwork(1)
      mu = iwork(2)
      if (ml .lt. 0 .or. ml .ge. n) go to 609
      if (mu .lt. 0 .or. mu .ge. n) go to 610
 30   continue
! next process and check the optional input. ---------------------------
      if (iopt .eq. 1) go to 40
      cmn12%maxord = mord(cmn12%meth)
      cmn12%mxstep = mxstp0
      cmn12%mxhnil = mxhnl0
      if (istate .eq. 1) h0 = zero
      cmn12%hmxi = zero
      cmn12%hmin = zero
      go to 60
 40   cmn12%maxord = iwork(5)
      if (cmn12%maxord .lt. 0) go to 611
      if (cmn12%maxord .eq. 0) cmn12%maxord = 100
      cmn12%maxord = min(cmn12%maxord,mord(cmn12%meth))
      cmn12%mxstep = iwork(6)
      if (cmn12%mxstep .lt. 0) go to 612
      if (cmn12%mxstep .eq. 0) cmn12%mxstep = mxstp0
      cmn12%mxhnil = iwork(7)
      if (cmn12%mxhnil .lt. 0) go to 613
      if (cmn12%mxhnil .eq. 0) cmn12%mxhnil = mxhnl0
      if (istate .ne. 1) go to 50
      h0 = rwork(5)
      if ((tout - t)*h0 .lt. zero) go to 614
 50   hmax = rwork(6)
      if (hmax .lt. zero) go to 615
      cmn12%hmxi = zero
      if (hmax .gt. zero) cmn12%hmxi = one/hmax
      cmn12%hmin = rwork(7)
      if (cmn12%hmin .lt. zero) go to 616
!-----------------------------------------------------------------------
! set work array pointers and check lengths lrw and liw.
! pointers to segments of rwork and iwork are named by prefixing l to
! the name of the segment.  e.g., the segment yh starts at rwork(lyh).
! segments of rwork (in order) are denoted  yh, wm, ewt, savf, acor.
! within wm, locjs is the location of the saved jacobian (jsv .gt. 0).
!-----------------------------------------------------------------------
 60   cmn12%lyh = 21
      if (istate .eq. 1) cmn12%nyh = n
      cmn12%lwm = cmn12%lyh + (cmn12%maxord + 1)*cmn12%nyh
      jco = max(0,cmn12%jsv)
      if (cmn12%miter .eq. 0) lenwm = 0
      if (cmn12%miter .eq. 1 .or. cmn12%miter .eq. 2) then
        lenwm = 2 + (1 + jco)*n*n
        cmn12%locjs = n*n + 3
      endif
      if (cmn12%miter .eq. 3) lenwm = 2 + n
      if (cmn12%miter .eq. 4 .or. cmn12%miter .eq. 5) then
        mband = ml + mu + 1
        lenp = (mband + ml)*n
        lenj = mband*n
        lenwm = 2 + lenp + jco*lenj
        cmn12%locjs = lenp + 3
        endif
      cmn12%lewt = cmn12%lwm + lenwm
      cmn12%lsavf = cmn12%lewt + n
      cmn12%lacor = cmn12%lsavf + n
      lenrw = cmn12%lacor + n - 1
      iwork(17) = lenrw
      cmn12%liwm = 1
      leniw = 30 + n
      if (cmn12%miter .eq. 0 .or. cmn12%miter .eq. 3) leniw = 30
      iwork(18) = leniw
      if (lenrw .gt. lrw) go to 617
      if (leniw .gt. liw) go to 618
! check rtol and atol for legality. ------------------------------------
      rtoli = rtol(1)
      atoli = atol(1)
      do 70 i = 1,n
        if (itol .ge. 3) rtoli = rtol(i)
        if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
        if (rtoli .lt. zero) go to 619
        if (atoli .lt. zero) go to 620
 70     continue
      if (istate .eq. 1) go to 100
! if istate = 3, set flag to signal parameter changes to dvstep. -------
      cmn12%jstart = -1
      if (cmn12%nq .le. cmn12%maxord) go to 90
! maxord was reduced below nq.  copy yh(*,maxord+2) into savf. ---------
      call dcopy (n, rwork(cmn12%lwm), 1, rwork(cmn12%lsavf), 1)
! reload wm(1) = rwork(lwm), since lwm may have changed. ---------------
 90   if (cmn12%miter .gt. 0) rwork(cmn12%lwm) = sqrt(cmn12%uround)
      go to 200
!-----------------------------------------------------------------------
! block c.
! the next block is for the initial call only (istate = 1).
! it contains all remaining initializations, the initial call to f,
! and the calculation of the initial step size.
! the error weights in ewt are inverted after being loaded.
!-----------------------------------------------------------------------
 100  cmn12%uround = dumach()
      cmn12%tn = t
      if (itask .ne. 4 .and. itask .ne. 5) go to 110
      tcrit = rwork(1)
      if ((tcrit - tout)*(tout - t) .lt. zero) go to 625
      if (h0 .ne. zero .and. (t + h0 - tcrit)*h0 .gt. zero)   &
         h0 = tcrit - t
 110  cmn12%jstart = 0
      if (cmn12%miter .gt. 0) rwork(cmn12%lwm) = sqrt(cmn12%uround)
      cmn12%ccmxj = pt2
      cmn12%msbj = 50
      cmn12%nhnil = 0
      cmn12%nst = 0
      cmn12%nje = 0
      cmn12%nni = 0
      cmn12%ncfn = 0
      cmn12%netf = 0
      cmn12%nlu = 0
      cmn12%nslj = 0
      nslast = 0
      cmn12%hu = zero
      cmn12%nqu = 0
! initial call to f.  (lf0 points to yh(*,2).) -------------------------
      lf0 = cmn12%lyh + cmn12%nyh
      call f (n, t, y, rwork(lf0), rpar, ipar)
      cmn12%nfe = 1
! load the initial value vector in yh. ---------------------------------
      call dcopy (n, y, 1, rwork(cmn12%lyh), 1)
! load and invert the ewt array.  (h is temporarily set to 1.0.) -------
      cmn12%nq = 1
      h = one
      call dewset (n, itol, rtol, atol, rwork(cmn12%lyh), rwork(cmn12%lewt))
      do 120 i = 1,n
        if (rwork(i+cmn12%lewt-1) .le. zero) go to 621
 120    rwork(i+cmn12%lewt-1) = one/rwork(i+cmn12%lewt-1)
      if (h0 .ne. zero) go to 180
! call dvhin to set initial step size h0 to be attempted. --------------
      call dvhin (n, t, rwork(cmn12%lyh), rwork(lf0), f, rpar, ipar, tout,   &
         cmn12%uround, rwork(cmn12%lewt), itol, atol, y, rwork(cmn12%lacor), h0,   &
         niter, ier)
      cmn12%nfe = cmn12%nfe + niter
      if (ier .ne. 0) go to 622
! adjust h0 if necessary to meet hmax bound. ---------------------------
 180  rh = abs(h0)*cmn12%hmxi
      if (rh .gt. one) h0 = h0/rh
! load h with h0 and scale yh(*,2) by h0. ------------------------------
      h = h0
      call dscal (n, h0, rwork(lf0), 1)
      go to 270
!-----------------------------------------------------------------------
! block d.
! the next code block is for continuation calls only (istate = 2 or 3)
! and is to check stop conditions before taking a step.
!-----------------------------------------------------------------------
 200  nslast = cmn12%nst
      cmn12%kuth = 0
      go to (210, 250, 220, 230, 240), itask
 210  if ((cmn12%tn - tout)*h .lt. zero) go to 250
      cmn12%h = h ; cmn12%n = n   ! load h, n into cmn12
      call dvindy (tout, 0, rwork(cmn12%lyh), cmn12%nyh, y, iflag, cmn12)
      h = cmn12%h ; n = cmn12%n   ! unload h, n from cmn12
      if (iflag .ne. 0) go to 627
      t = tout
      go to 420
 220  tp = cmn12%tn - cmn12%hu*(one + hun*cmn12%uround)
      if ((tp - tout)*h .gt. zero) go to 623
      if ((cmn12%tn - tout)*h .lt. zero) go to 250
      go to 400
 230  tcrit = rwork(1)
      if ((cmn12%tn - tcrit)*h .gt. zero) go to 624
      if ((tcrit - tout)*h .lt. zero) go to 625
      if ((cmn12%tn - tout)*h .lt. zero) go to 245
      cmn12%h = h ; cmn12%n = n   ! load h, n into cmn12
      call dvindy (tout, 0, rwork(cmn12%lyh), cmn12%nyh, y, iflag, cmn12)
      h = cmn12%h ; n = cmn12%n   ! unload h, n from cmn12
      if (iflag .ne. 0) go to 627
      t = tout
      go to 420
 240  tcrit = rwork(1)
      if ((cmn12%tn - tcrit)*h .gt. zero) go to 624
 245  hmx = abs(cmn12%tn) + abs(h)
      ihit = abs(cmn12%tn - tcrit) .le. hun*cmn12%uround*hmx
      if (ihit) go to 400
      tnext = cmn12%tn + cmn12%hnew*(one + four*cmn12%uround)
      if ((tnext - tcrit)*h .le. zero) go to 250
      h = (tcrit - cmn12%tn)*(one - four*cmn12%uround)
      cmn12%kuth = 1
!-----------------------------------------------------------------------
! block e.
! the next block is normally executed for all calls and contains
! the call to the one-step core integrator dvstep.
!
! this is a looping point for the integration steps.
!
! first check for too many steps being taken, update ewt (if not at
! start of problem), check for too much accuracy being requested, and
! check for h below the roundoff level in t.
!-----------------------------------------------------------------------
 250  continue
      if ((cmn12%nst-nslast) .ge. cmn12%mxstep) go to 500
      call dewset (n, itol, rtol, atol, rwork(cmn12%lyh), rwork(cmn12%lewt))
      do 260 i = 1,n
        if (rwork(i+cmn12%lewt-1) .le. zero) go to 510
 260    rwork(i+cmn12%lewt-1) = one/rwork(i+cmn12%lewt-1)
 270  tolsf = cmn12%uround*dvnorm (n, rwork(cmn12%lyh), rwork(cmn12%lewt))
      if (tolsf .le. one) go to 280
      tolsf = tolsf*two
      if (cmn12%nst .eq. 0) go to 626
      go to 520
 280  if ((cmn12%tn + h) .ne. cmn12%tn) go to 290
      cmn12%nhnil = cmn12%nhnil + 1
      if (cmn12%nhnil .gt. cmn12%mxhnil) go to 290
      msg = 'dvode--  warning: internal t (=r1) and h (=r2) are'
      call xerrwd (msg, 50, 101, 1, 0, 0, 0, 0, zero, zero)
      msg='      such that in the machine, t + h = t on the next step  '
      call xerrwd (msg, 60, 101, 1, 0, 0, 0, 0, zero, zero)
      msg = '      (h = step size). solver will continue anyway'
      call xerrwd (msg, 50, 101, 1, 0, 0, 0, 2, cmn12%tn, h)
      if (cmn12%nhnil .lt. cmn12%mxhnil) go to 290
      msg = 'dvode--  above warning has been issued i1 times.  '
      call xerrwd (msg, 50, 102, 1, 0, 0, 0, 0, zero, zero)
      msg = '      it will not be issued again for this problem'
      call xerrwd (msg, 50, 102, 1, 1, cmn12%mxhnil, 0, 0, zero, zero)
 290  continue
!-----------------------------------------------------------------------
! call dvstep (y, yh, nyh, yh, ewt, savf, vsav, acor,
!              wm, iwm, f, jac, f, dvnlsd, rpar, ipar)
!-----------------------------------------------------------------------
      cmn12%h = h ; cmn12%n = n   ! load h, n into cmn12
      call dvstep (y, rwork(cmn12%lyh), cmn12%nyh, rwork(cmn12%lyh), rwork(cmn12%lewt),   &
         rwork(cmn12%lsavf), y, rwork(cmn12%lacor), rwork(cmn12%lwm), iwork(cmn12%liwm),   &
         f, jac, f, dvnlsd, rpar, ipar, cmn12)
      h = cmn12%h ; n = cmn12%n   ! unload h, n from cmn12
      kgo = 1 - cmn12%kflag
! branch on kflag.  note: in this version, kflag can not be set to -3.
!  kflag .eq. 0,   -1,  -2
      go to (300, 530, 540), kgo
!-----------------------------------------------------------------------
! block f.
! the following block handles the case of a successful return from the
! core integrator (kflag = 0).  test for stop conditions.
!-----------------------------------------------------------------------
 300  cmn12%init = 1
      cmn12%kuth = 0
      go to (310, 400, 330, 340, 350), itask
! itask = 1.  if tout has been reached, interpolate. -------------------
 310  if ((cmn12%tn - tout)*h .lt. zero) go to 250
      cmn12%h = h ; cmn12%n = n   ! load h, n into cmn12
      call dvindy (tout, 0, rwork(cmn12%lyh), cmn12%nyh, y, iflag, cmn12)
      h = cmn12%h ; n = cmn12%n   ! unload h, n from cmn12
      t = tout
      go to 420
! itask = 3.  jump to exit if tout was reached. ------------------------
 330  if ((cmn12%tn - tout)*h .ge. zero) go to 400
      go to 250
! itask = 4.  see if tout or tcrit was reached.  adjust h if necessary.
 340  if ((cmn12%tn - tout)*h .lt. zero) go to 345
      cmn12%h = h ; cmn12%n = n   ! load h, n into cmn12
      call dvindy (tout, 0, rwork(cmn12%lyh), cmn12%nyh, y, iflag, cmn12)
      h = cmn12%h ; n = cmn12%n   ! unload h, n from cmn12
      t = tout
      go to 420
 345  hmx = abs(cmn12%tn) + abs(h)
      ihit = abs(cmn12%tn - tcrit) .le. hun*cmn12%uround*hmx
      if (ihit) go to 400
      tnext = cmn12%tn + cmn12%hnew*(one + four*cmn12%uround)
      if ((tnext - tcrit)*h .le. zero) go to 250
      h = (tcrit - cmn12%tn)*(one - four*cmn12%uround)
      cmn12%kuth = 1
      go to 250
! itask = 5.  see if tcrit was reached and jump to exit. ---------------
 350  hmx = abs(cmn12%tn) + abs(h)
      ihit = abs(cmn12%tn - tcrit) .le. hun*cmn12%uround*hmx
!-----------------------------------------------------------------------
! block g.
! the following block handles all successful returns from dvode.
! if itask .ne. 1, y is loaded from yh and t is set accordingly.
! istate is set to 2, and the optional output is loaded into the work
! arrays before returning.
!-----------------------------------------------------------------------
 400  continue
      call dcopy (n, rwork(cmn12%lyh), 1, y, 1)
      t = cmn12%tn
      if (itask .ne. 4 .and. itask .ne. 5) go to 420
      if (ihit) t = tcrit
 420  istate = 2
      rwork(11) = cmn12%hu
      rwork(12) = cmn12%hnew
      rwork(13) = cmn12%tn
      iwork(11) = cmn12%nst
      iwork(12) = cmn12%nfe
      iwork(13) = cmn12%nje
      iwork(14) = cmn12%nqu
      iwork(15) = cmn12%newq
      iwork(19) = cmn12%nlu
      iwork(20) = cmn12%nni
      iwork(21) = cmn12%ncfn
      iwork(22) = cmn12%netf
      return
!-----------------------------------------------------------------------
! block h.
! the following block handles all unsuccessful returns other than
! those for illegal input.  first the error message routine is called.
! if there was an error test or convergence test failure, imxer is set.
! then y is loaded from yh, and t is set to tn.
! the optional output is loaded into the work arrays before returning.
!-----------------------------------------------------------------------
! the maximum number of steps was taken before reaching tout. ----------
 500  msg = 'dvode--  at current t (=r1), mxstep (=i1) steps   '
      call xerrwd (msg, 50, 201, 1, 0, 0, 0, 0, zero, zero)
      msg = '      taken on this call before reaching tout     '
      call xerrwd (msg, 50, 201, 1, 1, cmn12%mxstep, 0, 1, cmn12%tn, zero)
      istate = -1
      go to 580
! ewt(i) .le. 0.0 for some i (not at start of problem). ----------------
 510  ewti = rwork(cmn12%lewt+i-1)
      msg = 'dvode--  at t (=r1), ewt(i1) has become r2 .le. 0.'
      call xerrwd (msg, 50, 202, 1, 1, i, 0, 2, cmn12%tn, ewti)
      istate = -6
      go to 580
! too much accuracy requested for machine precision. -------------------
 520  msg = 'dvode--  at t (=r1), too much accuracy requested  '
      call xerrwd (msg, 50, 203, 1, 0, 0, 0, 0, zero, zero)
      msg = '      for precision of machine:   see tolsf (=r2) '
      call xerrwd (msg, 50, 203, 1, 0, 0, 0, 2, cmn12%tn, tolsf)
      rwork(14) = tolsf
      istate = -2
      go to 580
! kflag = -1.  error test failed repeatedly or with abs(h) = hmin. -----
 530  msg = 'dvode--  at t(=r1) and step size h(=r2), the error'
      call xerrwd (msg, 50, 204, 1, 0, 0, 0, 0, zero, zero)
      msg = '      test failed repeatedly or with abs(h) = hmin'
      call xerrwd (msg, 50, 204, 1, 0, 0, 0, 2, cmn12%tn, h)
      istate = -4
      go to 560
! kflag = -2.  convergence failed repeatedly or with abs(h) = hmin. ----
 540  msg = 'dvode--  at t (=r1) and step size h (=r2), the    '
      call xerrwd (msg, 50, 205, 1, 0, 0, 0, 0, zero, zero)
      msg = '      corrector convergence failed repeatedly     '
      call xerrwd (msg, 50, 205, 1, 0, 0, 0, 0, zero, zero)
      msg = '      or with abs(h) = hmin   '
      call xerrwd (msg, 30, 205, 1, 0, 0, 0, 2, cmn12%tn, h)
      istate = -5
! compute imxer if relevant. -------------------------------------------
 560  big = zero
      imxer = 1
      do 570 i = 1,n
        size = abs(rwork(i+cmn12%lacor-1)*rwork(i+cmn12%lewt-1))
        if (big .ge. size) go to 570
        big = size
        imxer = i
 570    continue
      iwork(16) = imxer
! set y vector, t, and optional output. --------------------------------
 580  continue
      call dcopy (n, rwork(cmn12%lyh), 1, y, 1)
      t = cmn12%tn
      rwork(11) = cmn12%hu
      rwork(12) = h
      rwork(13) = cmn12%tn
      iwork(11) = cmn12%nst
      iwork(12) = cmn12%nfe
      iwork(13) = cmn12%nje
      iwork(14) = cmn12%nqu
      iwork(15) = cmn12%nq
      iwork(19) = cmn12%nlu
      iwork(20) = cmn12%nni
      iwork(21) = cmn12%ncfn
      iwork(22) = cmn12%netf
      return
!-----------------------------------------------------------------------
! block i.
! the following block handles all error returns due to illegal input
! (istate = -3), as detected before calling the core integrator.
! first the error message routine is called.   if the illegal input
! is a negative istate, the run is aborted (apparent infinite loop).
!-----------------------------------------------------------------------
 601  msg = 'dvode--  istate (=i1) illegal '
      call xerrwd (msg, 30, 1, 1, 1, istate, 0, 0, zero, zero)
      if (istate .lt. 0) go to 800
      go to 700
 602  msg = 'dvode--  itask (=i1) illegal  '
      call xerrwd (msg, 30, 2, 1, 1, itask, 0, 0, zero, zero)
      go to 700
 603  msg='dvode--  istate (=i1) .gt. 1 but dvode not initialized      '
      call xerrwd (msg, 60, 3, 1, 1, istate, 0, 0, zero, zero)
      go to 700
 604  msg = 'dvode--  neq (=i1) .lt. 1     '
      call xerrwd (msg, 30, 4, 1, 1, neq, 0, 0, zero, zero)
      go to 700
 605  msg = 'dvode--  istate = 3 and neq increased (i1 to i2)  '
      call xerrwd (msg, 50, 5, 1, 2, n, neq, 0, zero, zero)
      go to 700
 606  msg = 'dvode--  itol (=i1) illegal   '
      call xerrwd (msg, 30, 6, 1, 1, itol, 0, 0, zero, zero)
      go to 700
 607  msg = 'dvode--  iopt (=i1) illegal   '
      call xerrwd (msg, 30, 7, 1, 1, iopt, 0, 0, zero, zero)
      go to 700
 608  msg = 'dvode--  mf (=i1) illegal     '
      call xerrwd (msg, 30, 8, 1, 1, mf, 0, 0, zero, zero)
      go to 700
 609  msg = 'dvode--  ml (=i1) illegal:  .lt.0 or .ge.neq (=i2)'
      call xerrwd (msg, 50, 9, 1, 2, ml, neq, 0, zero, zero)
      go to 700
 610  msg = 'dvode--  mu (=i1) illegal:  .lt.0 or .ge.neq (=i2)'
      call xerrwd (msg, 50, 10, 1, 2, mu, neq, 0, zero, zero)
      go to 700
 611  msg = 'dvode--  maxord (=i1) .lt. 0  '
      call xerrwd (msg, 30, 11, 1, 1, cmn12%maxord, 0, 0, zero, zero)
      go to 700
 612  msg = 'dvode--  mxstep (=i1) .lt. 0  '
      call xerrwd (msg, 30, 12, 1, 1, cmn12%mxstep, 0, 0, zero, zero)
      go to 700
 613  msg = 'dvode--  mxhnil (=i1) .lt. 0  '
      call xerrwd (msg, 30, 13, 1, 1, cmn12%mxhnil, 0, 0, zero, zero)
      go to 700
 614  msg = 'dvode--  tout (=r1) behind t (=r2)      '
      call xerrwd (msg, 40, 14, 1, 0, 0, 0, 2, tout, t)
      msg = '      integration direction is given by h0 (=r1)  '
      call xerrwd (msg, 50, 14, 1, 0, 0, 0, 1, h0, zero)
      go to 700
 615  msg = 'dvode--  hmax (=r1) .lt. 0.0  '
      call xerrwd (msg, 30, 15, 1, 0, 0, 0, 1, hmax, zero)
      go to 700
 616  msg = 'dvode--  hmin (=r1) .lt. 0.0  '
      call xerrwd (msg, 30, 16, 1, 0, 0, 0, 1, cmn12%hmin, zero)
      go to 700
 617  continue
      msg='dvode--  rwork length needed, lenrw (=i1), exceeds lrw (=i2)'
      call xerrwd (msg, 60, 17, 1, 2, lenrw, lrw, 0, zero, zero)
      go to 700
 618  continue
      msg='dvode--  iwork length needed, leniw (=i1), exceeds liw (=i2)'
      call xerrwd (msg, 60, 18, 1, 2, leniw, liw, 0, zero, zero)
      go to 700
 619  msg = 'dvode--  rtol(i1) is r1 .lt. 0.0        '
      call xerrwd (msg, 40, 19, 1, 1, i, 0, 1, rtoli, zero)
      go to 700
 620  msg = 'dvode--  atol(i1) is r1 .lt. 0.0        '
      call xerrwd (msg, 40, 20, 1, 1, i, 0, 1, atoli, zero)
      go to 700
 621  ewti = rwork(cmn12%lewt+i-1)
      msg = 'dvode--  ewt(i1) is r1 .le. 0.0         '
      call xerrwd (msg, 40, 21, 1, 1, i, 0, 1, ewti, zero)
      go to 700
 622  continue
      msg='dvode--  tout (=r1) too close to t(=r2) to start integration'
      call xerrwd (msg, 60, 22, 1, 0, 0, 0, 2, tout, t)
      go to 700
 623  continue
      msg='dvode--  itask = i1 and tout (=r1) behind tcur - hu (= r2)  '
      call xerrwd (msg, 60, 23, 1, 1, itask, 0, 2, tout, tp)
      go to 700
 624  continue
      msg='dvode--  itask = 4 or 5 and tcrit (=r1) behind tcur (=r2)   '
      call xerrwd (msg, 60, 24, 1, 0, 0, 0, 2, tcrit, cmn12%tn)
      go to 700
 625  continue
      msg='dvode--  itask = 4 or 5 and tcrit (=r1) behind tout (=r2)   '
      call xerrwd (msg, 60, 25, 1, 0, 0, 0, 2, tcrit, tout)
      go to 700
 626  msg = 'dvode--  at start of problem, too much accuracy   '
      call xerrwd (msg, 50, 26, 1, 0, 0, 0, 0, zero, zero)
      msg='      requested for precision of machine:   see tolsf (=r1) '
      call xerrwd (msg, 60, 26, 1, 0, 0, 0, 1, tolsf, zero)
      rwork(14) = tolsf
      go to 700
 627  msg='dvode--  trouble from dvindy.  itask = i1, tout = r1.       '
      call xerrwd (msg, 60, 27, 1, 1, itask, 0, 1, tout, zero)
!
 700  continue
      istate = -3
      return
!
 800  msg = 'dvode--  run aborted:  apparent infinite loop     '
      call xerrwd (msg, 50, 303, 2, 0, 0, 0, 0, zero, zero)
      return
!----------------------- end of subroutine dvode -----------------------
      end subroutine dvode
!deck dvhin
      subroutine dvhin (n, t0, y0, ydot, f, rpar, ipar, tout, uround,   &
         ewt, itol, atol, y, temp, h0, niter, ier)
      external f
      double precision t0, y0, ydot, rpar, tout, uround, ewt, atol, y,   &
         temp, h0
      integer n, ipar, itol, niter, ier
      dimension y0(*), ydot(*), ewt(*), atol(*), y(*),   &
         temp(*), rpar(*), ipar(*)
!-----------------------------------------------------------------------
! call sequence input -- n, t0, y0, ydot, f, rpar, ipar, tout, uround,
!                        ewt, itol, atol, y, temp
! call sequence output -- h0, niter, ier
! common block variables accessed -- none
!
! subroutines called by dvhin:  f
! function routines called by dvhi: dvnorm
!-----------------------------------------------------------------------
! this routine computes the step size, h0, to be attempted on the
! first step, when the user has not supplied a value for this.
!
! first we check that tout - t0 differs significantly from zero.  then
! an iteration is done to approximate the initial second derivative
! and this is used to define h from w.r.m.s.norm(h**2 * yddot / 2) = 1.
! a bias factor of 1/2 is applied to the resulting h.
! the sign of h0 is inferred from the initial values of tout and t0.
!
! communication with dvhin is done with the following variables:
!
! n      = size of ode system, input.
! t0     = initial value of independent variable, input.
! y0     = vector of initial conditions, input.
! ydot   = vector of initial first derivatives, input.
! f      = name of subroutine for right-hand side f(t,y), input.
! rpar, ipar = dummy names for user's real and integer work arrays.
! tout   = first output value of independent variable
! uround = machine unit roundoff
! ewt, itol, atol = error weights and tolerance parameters
!                   as described in the driver routine, input.
! y, temp = work arrays of length n.
! h0     = step size to be attempted, output.
! niter  = number of iterations (and of f evaluations) to compute h0,
!          output.
! ier    = the error flag, returned with the value
!          ier = 0  if no trouble occurred, or
!          ier = -1 if tout and t0 are considered too close to proceed.
!-----------------------------------------------------------------------
!
! type declarations for local variables --------------------------------
!
      double precision afi, atoli, delyi, h, half, hg, hlb, hnew, hrat,   &
           hub, hun, pt1, t1, tdist, tround, two, yddnrm
      integer i, iter
!
! type declaration for function subroutines called ---------------------
!
! 27-oct-2005 rce - do not declare functions that are in the module
!     double precision dvnorm
!-----------------------------------------------------------------------
! the following fortran-77 declaration is to cause the values of the
! listed (local) variables to be saved between calls to this integrator.
!-----------------------------------------------------------------------
      save half, hun, pt1, two
      data half /0.5d0/, hun /100.0d0/, pt1 /0.1d0/, two /2.0d0/
!
      niter = 0
      tdist = abs(tout - t0)
      tround = uround*max(abs(t0),abs(tout))
      if (tdist .lt. two*tround) go to 100
!
! set a lower bound on h based on the roundoff level in t0 and tout. ---
      hlb = hun*tround
! set an upper bound on h based on tout-t0 and the initial y and ydot. -
      hub = pt1*tdist
      atoli = atol(1)
      do 10 i = 1, n
        if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
        delyi = pt1*abs(y0(i)) + atoli
        afi = abs(ydot(i))
        if (afi*hub .gt. delyi) hub = delyi/afi
 10     continue
!
! set initial guess for h as geometric mean of upper and lower bounds. -
      iter = 0
      hg = sqrt(hlb*hub)
! if the bounds have crossed, exit with the mean value. ----------------
      if (hub .lt. hlb) then
        h0 = hg
        go to 90
      endif
!
! looping point for iteration. -----------------------------------------
 50   continue
! estimate the second derivative as a difference quotient in f. --------
      h = sign (hg, tout - t0)
      t1 = t0 + h
      do 60 i = 1, n
 60     y(i) = y0(i) + h*ydot(i)
      call f (n, t1, y, temp, rpar, ipar)
      do 70 i = 1, n
 70     temp(i) = (temp(i) - ydot(i))/h
      yddnrm = dvnorm (n, temp, ewt)
! get the corresponding new value of h. --------------------------------
      if (yddnrm*hub*hub .gt. two) then
        hnew = sqrt(two/yddnrm)
      else
        hnew = sqrt(hg*hub)
      endif
      iter = iter + 1
!-----------------------------------------------------------------------
! test the stopping conditions.
! stop if the new and previous h values differ by a factor of .lt. 2.
! stop if four iterations have been done.  also, stop with previous h
! if hnew/hg .gt. 2 after first iteration, as this probably means that
! the second derivative value is bad because of cancellation error.
!-----------------------------------------------------------------------
      if (iter .ge. 4) go to 80
      hrat = hnew/hg
      if ( (hrat .gt. half) .and. (hrat .lt. two) ) go to 80
      if ( (iter .ge. 2) .and. (hnew .gt. two*hg) ) then
        hnew = hg
        go to 80
      endif
      hg = hnew
      go to 50
!
! iteration done.  apply bounds, bias factor, and sign.  then exit. ----
 80   h0 = hnew*half
      if (h0 .lt. hlb) h0 = hlb
      if (h0 .gt. hub) h0 = hub
 90   h0 = sign(h0, tout - t0)
      niter = iter
      ier = 0
      return
! error return for tout - t0 too small. --------------------------------
 100  ier = -1
      return
!----------------------- end of subroutine dvhin -----------------------
      end subroutine dvhin
!deck dvindy
      subroutine dvindy (t, k, yh, ldyh, dky, iflag, cmn12)
      implicit none
      double precision t, yh, dky
      integer k, ldyh, iflag
      dimension yh(ldyh,*), dky(*)
      type( dvode_cmn_vars ) :: cmn12
!-----------------------------------------------------------------------
! call sequence input -- t, k, yh, ldyh
! call sequence output -- dky, iflag
! common block variables accessed:
!     /dvod_cmn01/ --  h, tn, uround, l, n, nq
!     /dvod_cmn02/ --  hu
!
! subroutines called by dvindy: dscal, xerrwd
! function routines called by dvindy: none
!-----------------------------------------------------------------------
! dvindy computes interpolated values of the k-th derivative of the
! dependent variable vector y, and stores it in dky.  this routine
! is called within the package with k = 0 and t = tout, but may
! also be called by the user for any k up to the current order.
! (see detailed instructions in the usage documentation.)
!-----------------------------------------------------------------------
! the computed values in dky are gotten by interpolation using the
! nordsieck history array yh.  this array corresponds uniquely to a
! vector-valued polynomial of degree nqcur or less, and dky is set
! to the k-th derivative of this polynomial at t.
! the formula for dky is:
!              q
!  dky(i)  =  sum  c(j,k) * (t - tn)**(j-k) * h**(-j) * yh(i,j+1)
!             j=k
! where  c(j,k) = j*(j-1)*...*(j-k+1), q = nqcur, tn = tcur, h = hcur.
! the quantities  nq = nqcur, l = nq+1, n, tn, and h are
! communicated by common.  the above sum is done in reverse order.
! iflag is returned negative if either k or t is out of bounds.
!
! discussion above and comments in driver explain all variables.
!-----------------------------------------------------------------------
!!
!! type declarations for labeled common block dvod_cmn01 --------------------
!!
!      double precision acnrm, ccmxj, conp, crate, drc, el,   &
!           eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!           rc, rl1, tau, tq, tn, uround
!      integer icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!              l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!              locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!              n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!              nslp, nyh
!!
!! type declarations for labeled common block dvod_cmn02 --------------------
!!
!      double precision hu
!      integer ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
! type declarations for local variables --------------------------------
!
      double precision c, hun, r, s, tfuzz, tn1, tp, zero
      integer i, ic, j, jb, jb2, jj, jj1, jp1
      character*80 msg
!-----------------------------------------------------------------------
! the following fortran-77 declaration is to cause the values of the
! listed (local) variables to be saved between calls to this integrator.
!-----------------------------------------------------------------------
      save hun, zero
!
!      common /dvod_cmn01/ acnrm, ccmxj, conp, crate, drc, el(13),   &
!                      eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!                      rc, rl1, tau(13), tq(5), tn, uround,   &
!                      icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!                      l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!                      locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!                      n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!                      nslp, nyh
!      common /dvod_cmn02/ hu, ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
      data hun /100.0d0/, zero /0.0d0/
!
      iflag = 0
      if (k .lt. 0 .or. k .gt. cmn12%nq) go to 80
      tfuzz = hun*cmn12%uround*(cmn12%tn + cmn12%hu)
      tp = cmn12%tn - cmn12%hu - tfuzz
      tn1 = cmn12%tn + tfuzz
      if ((t-tp)*(t-tn1) .gt. zero) go to 90
!
      s = (t - cmn12%tn)/cmn12%h
      ic = 1
      if (k .eq. 0) go to 15
      jj1 = cmn12%l - k
      do 10 jj = jj1, cmn12%nq
 10     ic = ic*jj
 15   c = real(ic)
      do 20 i = 1, cmn12%n
 20     dky(i) = c*yh(i,cmn12%l)
      if (k .eq. cmn12%nq) go to 55
      jb2 = cmn12%nq - k
      do 50 jb = 1, jb2
        j = cmn12%nq - jb
        jp1 = j + 1
        ic = 1
        if (k .eq. 0) go to 35
        jj1 = jp1 - k
        do 30 jj = jj1, j
 30       ic = ic*jj
 35     c = real(ic)
        do 40 i = 1, cmn12%n
 40       dky(i) = c*yh(i,jp1) + s*dky(i)
 50     continue
      if (k .eq. 0) return
 55   r = cmn12%h**(-k)
      call dscal (cmn12%n, r, dky, 1)
      return
!
 80   msg = 'dvindy-- k (=i1) illegal      '
      call xerrwd (msg, 30, 51, 1, 1, k, 0, 0, zero, zero)
      iflag = -1
      return
 90   msg = 'dvindy-- t (=r1) illegal      '
      call xerrwd (msg, 30, 52, 1, 0, 0, 0, 1, t, zero)
      msg='      t not in interval tcur - hu (= r1) to tcur (=r2)      '
      call xerrwd (msg, 60, 52, 1, 0, 0, 0, 2, tp, cmn12%tn)
      iflag = -2
      return
!----------------------- end of subroutine dvindy ----------------------
      end subroutine dvindy
!deck dvstep
      subroutine dvstep (y, yh, ldyh, yh1, ewt, savf, vsav, acor,   &
                        wm, iwm, f, jac, psol, vnls, rpar, ipar, cmn12)
      implicit none
      external f, jac, psol, vnls
      double precision y, yh, yh1, ewt, savf, vsav, acor, wm, rpar
      integer ldyh, iwm, ipar
      dimension y(*), yh(ldyh,*), yh1(*), ewt(*), savf(*), vsav(*),   &
         acor(*), wm(*), iwm(*), rpar(*), ipar(*)
      type( dvode_cmn_vars ) :: cmn12
!-----------------------------------------------------------------------
! call sequence input -- y, yh, ldyh, yh1, ewt, savf, vsav,
!                        acor, wm, iwm, f, jac, psol, vnls, rpar, ipar
! call sequence output -- yh, acor, wm, iwm
! common block variables accessed:
!     /dvod_cmn01/  acnrm, el(13), h, hmin, hmxi, hnew, hscal, rc, tau(13),
!               tq(5), tn, jcur, jstart, kflag, kuth,
!               l, lmax, maxord, n, newq, nq, nqwait
!     /dvod_cmn02/  hu, ncfn, netf, nfe, nqu, nst
!
! subroutines called by dvstep: f, daxpy, dcopy, dscal,
!                               dvjust, vnls, dvset
! function routines called by dvstep: dvnorm
!-----------------------------------------------------------------------
! dvstep performs one step of the integration of an initial value
! problem for a system of ordinary differential equations.
! dvstep calls subroutine vnls for the solution of the nonlinear system
! arising in the time step.  thus it is independent of the problem
! jacobian structure and the type of nonlinear system solution method.
! dvstep returns a completion flag kflag (in common).
! a return with kflag = -1 or -2 means either abs(h) = hmin or 10
! consecutive failures occurred.  on a return with kflag negative,
! the values of tn and the yh array are as of the beginning of the last
! step, and h is the last step size attempted.
!
! communication with dvstep is done with the following variables:
!
! y      = an array of length n used for the dependent variable vector.
! yh     = an ldyh by lmax array containing the dependent variables
!          and their approximate scaled derivatives, where
!          lmax = maxord + 1.  yh(i,j+1) contains the approximate
!          j-th derivative of y(i), scaled by h**j/factorial(j)
!          (j = 0,1,...,nq).  on entry for the first step, the first
!          two columns of yh must be set from the initial values.
! ldyh   = a constant integer .ge. n, the first dimension of yh.
!          n is the number of odes in the system.
! yh1    = a one-dimensional array occupying the same space as yh.
! ewt    = an array of length n containing multiplicative weights
!          for local error measurements.  local errors in y(i) are
!          compared to 1.0/ewt(i) in various error tests.
! savf   = an array of working storage, of length n.
!          also used for input of yh(*,maxord+2) when jstart = -1
!          and maxord .lt. the current order nq.
! vsav   = a work array of length n passed to subroutine vnls.
! acor   = a work array of length n, used for the accumulated
!          corrections.  on a successful return, acor(i) contains
!          the estimated one-step local error in y(i).
! wm,iwm = real and integer work arrays associated with matrix
!          operations in vnls.
! f      = dummy name for the user supplied subroutine for f.
! jac    = dummy name for the user supplied jacobian subroutine.
! psol   = dummy name for the subroutine passed to vnls, for
!          possible use there.
! vnls   = dummy name for the nonlinear system solving subroutine,
!          whose real name is dependent on the method used.
! rpar, ipar = dummy names for user's real and integer work arrays.
!-----------------------------------------------------------------------
!!
!! type declarations for labeled common block dvod_cmn01 --------------------
!!
!      double precision acnrm, ccmxj, conp, crate, drc, el,   &
!           eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!           rc, rl1, tau, tq, tn, uround
!      integer icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!              l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!              locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!              n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!              nslp, nyh
!!
!! type declarations for labeled common block dvod_cmn02 --------------------
!!
!      double precision hu
!      integer ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
! type declarations for local variables --------------------------------
!
! note -- etaq & etaqm1 are saved variables that are not fixed constants,
!         so put them into cmn12
!     double precision addon, bias1,bias2,bias3, cnquot, ddn, dsm, dup,   &
!          etacf, etamin, etamx1, etamx2, etamx3, etamxf,   &
!          etaq, etaqm1, etaqp1, flotl, one, onepsm,   &
!          r, thresh, told, zero
      double precision addon, bias1,bias2,bias3, cnquot, ddn, dsm, dup,   &
           etacf, etamin, etamx1, etamx2, etamx3, etamxf,   &
           etaqp1, flotl, one, onepsm,   &
           r, thresh, told, zero
      integer i, i1, i2, iback, j, jb, kfc, kfh, mxncf, ncf, nflag
!
! type declaration for function subroutines called ---------------------
!
! 27-oct-2005 rce - do not declare functions that are in the module
!     double precision dvnorm
!-----------------------------------------------------------------------
! the following fortran-77 declaration is to cause the values of the
! listed (local) variables to be saved between calls to this integrator.
!-----------------------------------------------------------------------
! note -- etaq & etaqm1 are saved variables that are not fixed constants,
!         so put them into cmn12
!     save addon, bias1, bias2, bias3,   &
!          etacf, etamin, etamx1, etamx2, etamx3, etamxf, etaq, etaqm1,   &
!          kfc, kfh, mxncf, onepsm, thresh, one, zero
      save addon, bias1, bias2, bias3,   &
           etacf, etamin, etamx1, etamx2, etamx3, etamxf,   &
           kfc, kfh, mxncf, onepsm, thresh, one, zero
!-----------------------------------------------------------------------
!      common /dvod_cmn01/ acnrm, ccmxj, conp, crate, drc, el(13),   &
!                      eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!                      rc, rl1, tau(13), tq(5), tn, uround,   &
!                      icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!                      l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!                      locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!                      n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!                      nslp, nyh
!      common /dvod_cmn02/ hu, ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
      data kfc/-3/, kfh/-7/, mxncf/10/
      data addon  /1.0d-6/,    bias1  /6.0d0/,     bias2  /6.0d0/,   &
           bias3  /10.0d0/,    etacf  /0.25d0/,    etamin /0.1d0/,   &
           etamxf /0.2d0/,     etamx1 /1.0d4/,     etamx2 /10.0d0/,   &
           etamx3 /10.0d0/,    onepsm /1.00001d0/, thresh /1.5d0/
      data one/1.0d0/, zero/0.0d0/
!
      cmn12%kflag = 0
      told = cmn12%tn
      ncf = 0
      cmn12%jcur = 0
      nflag = 0
      if (cmn12%jstart .gt. 0) go to 20
      if (cmn12%jstart .eq. -1) go to 100
!-----------------------------------------------------------------------
! on the first call, the order is set to 1, and other variables are
! initialized.  etamax is the maximum ratio by which h can be increased
! in a single step.  it is normally 10, but is larger during the
! first step to compensate for the small initial h.  if a failure
! occurs (in corrector convergence or error test), etamax is set to 1
! for the next increase.
!-----------------------------------------------------------------------
      cmn12%lmax = cmn12%maxord + 1
      cmn12%nq = 1
      cmn12%l = 2
      cmn12%nqnyh = cmn12%nq*ldyh
      cmn12%tau(1) = cmn12%h
      cmn12%prl1 = one
      cmn12%rc = zero
      cmn12%etamax = etamx1
      cmn12%nqwait = 2
      cmn12%hscal = cmn12%h
      go to 200
!-----------------------------------------------------------------------
! take preliminary actions on a normal continuation step (jstart.gt.0).
! if the driver changed h, then eta must be reset and newh set to 1.
! if a change of order was dictated on the previous step, then
! it is done here and appropriate adjustments in the history are made.
! on an order decrease, the history array is adjusted by dvjust.
! on an order increase, the history array is augmented by a column.
! on a change of step size h, the history array yh is rescaled.
!-----------------------------------------------------------------------
 20   continue
      if (cmn12%kuth .eq. 1) then
        cmn12%eta = min(cmn12%eta,cmn12%h/cmn12%hscal)
        cmn12%newh = 1
        endif
 50   if (cmn12%newh .eq. 0) go to 200
      if (cmn12%newq .eq. cmn12%nq) go to 150
      if (cmn12%newq .lt. cmn12%nq) then
        call dvjust (yh, ldyh, -1, cmn12)
        cmn12%nq = cmn12%newq
        cmn12%l = cmn12%nq + 1
        cmn12%nqwait = cmn12%l
        go to 150
        endif
      if (cmn12%newq .gt. cmn12%nq) then
        call dvjust (yh, ldyh, 1, cmn12)
        cmn12%nq = cmn12%newq
        cmn12%l = cmn12%nq + 1
        cmn12%nqwait = cmn12%l
        go to 150
      endif
!-----------------------------------------------------------------------
! the following block handles preliminaries needed when jstart = -1.
! if n was reduced, zero out part of yh to avoid undefined references.
! if maxord was reduced to a value less than the tentative order newq,
! then nq is set to maxord, and a new h ratio eta is chosen.
! otherwise, we take the same preliminary actions as for jstart .gt. 0.
! in any case, nqwait is reset to l = nq + 1 to prevent further
! changes in order for that many steps.
! the new h ratio eta is limited by the input h if kuth = 1,
! by hmin if kuth = 0, and by hmxi in any case.
! finally, the history array yh is rescaled.
!-----------------------------------------------------------------------
 100  continue
      cmn12%lmax = cmn12%maxord + 1
      if (cmn12%n .eq. ldyh) go to 120
      i1 = 1 + (cmn12%newq + 1)*ldyh
      i2 = (cmn12%maxord + 1)*ldyh
      if (i1 .gt. i2) go to 120
      do 110 i = i1, i2
 110    yh1(i) = zero
 120  if (cmn12%newq .le. cmn12%maxord) go to 140
      flotl = real(cmn12%lmax)
      if (cmn12%maxord .lt. cmn12%nq-1) then
        ddn = dvnorm (cmn12%n, savf, ewt)/cmn12%tq(1)
        cmn12%eta = one/((bias1*ddn)**(one/flotl) + addon)
        endif
      if (cmn12%maxord .eq. cmn12%nq .and. cmn12%newq .eq. cmn12%nq+1) cmn12%eta = cmn12%etaq
      if (cmn12%maxord .eq. cmn12%nq-1 .and. cmn12%newq .eq. cmn12%nq+1) then
        cmn12%eta = cmn12%etaqm1
        call dvjust (yh, ldyh, -1, cmn12)
        endif
      if (cmn12%maxord .eq. cmn12%nq-1 .and. cmn12%newq .eq. cmn12%nq) then
        ddn = dvnorm (cmn12%n, savf, ewt)/cmn12%tq(1)
        cmn12%eta = one/((bias1*ddn)**(one/flotl) + addon)
        call dvjust (yh, ldyh, -1, cmn12)
        endif
      cmn12%eta = min(cmn12%eta,one)
      cmn12%nq = cmn12%maxord
      cmn12%l = cmn12%lmax
 140  if (cmn12%kuth .eq. 1) cmn12%eta = min(cmn12%eta,abs(cmn12%h/cmn12%hscal))
      if (cmn12%kuth .eq. 0) cmn12%eta = max(cmn12%eta,cmn12%hmin/abs(cmn12%hscal))
      cmn12%eta = cmn12%eta/max(one,abs(cmn12%hscal)*cmn12%hmxi*cmn12%eta)
      cmn12%newh = 1
      cmn12%nqwait = cmn12%l
      if (cmn12%newq .le. cmn12%maxord) go to 50
! rescale the history array for a change in h by a factor of eta. ------
 150  r = one
      do 180 j = 2, cmn12%l
        r = r*cmn12%eta
        call dscal (cmn12%n, r, yh(1,j), 1 )
 180    continue
      cmn12%h = cmn12%hscal*cmn12%eta
      cmn12%hscal = cmn12%h
      cmn12%rc = cmn12%rc*cmn12%eta
      cmn12%nqnyh = cmn12%nq*ldyh
!-----------------------------------------------------------------------
! this section computes the predicted values by effectively
! multiplying the yh array by the pascal triangle matrix.
! dvset is called to calculate all integration coefficients.
! rc is the ratio of new to old values of the coefficient h/el(2)=h/l1.
!-----------------------------------------------------------------------
 200  cmn12%tn = cmn12%tn + cmn12%h
      i1 = cmn12%nqnyh + 1
      do 220 jb = 1, cmn12%nq
        i1 = i1 - ldyh
        do 210 i = i1, cmn12%nqnyh
 210      yh1(i) = yh1(i) + yh1(i+ldyh)
 220  continue
      call dvset( cmn12 )
      cmn12%rl1 = one/cmn12%el(2)
      cmn12%rc = cmn12%rc*(cmn12%rl1/cmn12%prl1)
      cmn12%prl1 = cmn12%rl1
!
! call the nonlinear system solver. ------------------------------------
!
      call vnls (y, yh, ldyh, vsav, savf, ewt, acor, iwm, wm,   &
                 f, jac, psol, nflag, rpar, ipar, cmn12)
!
      if (nflag .eq. 0) go to 450
!-----------------------------------------------------------------------
! the vnls routine failed to achieve convergence (nflag .ne. 0).
! the yh array is retracted to its values before prediction.
! the step size h is reduced and the step is retried, if possible.
! otherwise, an error exit is taken.
!-----------------------------------------------------------------------
        ncf = ncf + 1
        cmn12%ncfn = cmn12%ncfn + 1
        cmn12%etamax = one
        cmn12%tn = told
        i1 = cmn12%nqnyh + 1
        do 430 jb = 1, cmn12%nq
          i1 = i1 - ldyh
          do 420 i = i1, cmn12%nqnyh
 420        yh1(i) = yh1(i) - yh1(i+ldyh)
 430      continue
        if (nflag .lt. -1) go to 680
        if (abs(cmn12%h) .le. cmn12%hmin*onepsm) go to 670
        if (ncf .eq. mxncf) go to 670
        cmn12%eta = etacf
        cmn12%eta = max(cmn12%eta,cmn12%hmin/abs(cmn12%h))
        nflag = -1
        go to 150
!-----------------------------------------------------------------------
! the corrector has converged (nflag = 0).  the local error test is
! made and control passes to statement 500 if it fails.
!-----------------------------------------------------------------------
 450  continue
      dsm = cmn12%acnrm/cmn12%tq(2)
      if (dsm .gt. one) go to 500
!-----------------------------------------------------------------------
! after a successful step, update the yh and tau arrays and decrement
! nqwait.  if nqwait is then 1 and nq .lt. maxord, then acor is saved
! for use in a possible order increase on the next step.
! if etamax = 1 (a failure occurred this step), keep nqwait .ge. 2.
!-----------------------------------------------------------------------
      cmn12%kflag = 0
      cmn12%nst = cmn12%nst + 1
      cmn12%hu = cmn12%h
      cmn12%nqu = cmn12%nq
      do 470 iback = 1, cmn12%nq
        i = cmn12%l - iback
 470    cmn12%tau(i+1) = cmn12%tau(i)
      cmn12%tau(1) = cmn12%h
      do 480 j = 1, cmn12%l
        call daxpy (cmn12%n, cmn12%el(j), acor, 1, yh(1,j), 1 )
 480    continue
      cmn12%nqwait = cmn12%nqwait - 1
      if ((cmn12%l .eq. cmn12%lmax) .or. (cmn12%nqwait .ne. 1)) go to 490
      call dcopy (cmn12%n, acor, 1, yh(1,cmn12%lmax), 1 )
      cmn12%conp = cmn12%tq(5)
 490  if (cmn12%etamax .ne. one) go to 560
      if (cmn12%nqwait .lt. 2) cmn12%nqwait = 2
      cmn12%newq = cmn12%nq
      cmn12%newh = 0
      cmn12%eta = one
      cmn12%hnew = cmn12%h
      go to 690
!-----------------------------------------------------------------------
! the error test failed.  kflag keeps track of multiple failures.
! restore tn and the yh array to their previous values, and prepare
! to try the step again.  compute the optimum step size for the
! same order.  after repeated failures, h is forced to decrease
! more rapidly.
!-----------------------------------------------------------------------
 500  cmn12%kflag = cmn12%kflag - 1
      cmn12%netf = cmn12%netf + 1
      nflag = -2
      cmn12%tn = told
      i1 = cmn12%nqnyh + 1
      do 520 jb = 1, cmn12%nq
        i1 = i1 - ldyh
        do 510 i = i1, cmn12%nqnyh
 510      yh1(i) = yh1(i) - yh1(i+ldyh)
 520  continue
      if (abs(cmn12%h) .le. cmn12%hmin*onepsm) go to 660
      cmn12%etamax = one
      if (cmn12%kflag .le. kfc) go to 530
! compute ratio of new h to current h at the current order. ------------
      flotl = real(cmn12%l)
      cmn12%eta = one/((bias2*dsm)**(one/flotl) + addon)
      cmn12%eta = max(cmn12%eta,cmn12%hmin/abs(cmn12%h),etamin)
      if ((cmn12%kflag .le. -2) .and. (cmn12%eta .gt. etamxf)) cmn12%eta = etamxf
      go to 150
!-----------------------------------------------------------------------
! control reaches this section if 3 or more consecutive failures
! have occurred.  it is assumed that the elements of the yh array
! have accumulated errors of the wrong order.  the order is reduced
! by one, if possible.  then h is reduced by a factor of 0.1 and
! the step is retried.  after a total of 7 consecutive failures,
! an exit is taken with kflag = -1.
!-----------------------------------------------------------------------
 530  if (cmn12%kflag .eq. kfh) go to 660
      if (cmn12%nq .eq. 1) go to 540
      cmn12%eta = max(etamin,cmn12%hmin/abs(cmn12%h))
      call dvjust (yh, ldyh, -1, cmn12)
      cmn12%l = cmn12%nq
      cmn12%nq = cmn12%nq - 1
      cmn12%nqwait = cmn12%l
      go to 150
 540  cmn12%eta = max(etamin,cmn12%hmin/abs(cmn12%h))
      cmn12%h = cmn12%h*cmn12%eta
      cmn12%hscal = cmn12%h
      cmn12%tau(1) = cmn12%h
      call f (cmn12%n, cmn12%tn, y, savf, rpar, ipar)
      cmn12%nfe = cmn12%nfe + 1
      do 550 i = 1, cmn12%n
 550    yh(i,2) = cmn12%h*savf(i)
      cmn12%nqwait = 10
      go to 200
!-----------------------------------------------------------------------
! if nqwait = 0, an increase or decrease in order by one is considered.
! factors etaq, etaqm1, etaqp1 are computed by which h could
! be multiplied at order q, q-1, or q+1, respectively.
! the largest of these is determined, and the new order and
! step size set accordingly.
! a change of h or nq is made only if h increases by at least a
! factor of thresh.  if an order change is considered and rejected,
! then nqwait is set to 2 (reconsider it after 2 steps).
!-----------------------------------------------------------------------
! compute ratio of new h to current h at the current order. ------------
 560  flotl = real(cmn12%l)
      cmn12%etaq = one/((bias2*dsm)**(one/flotl) + addon)
      if (cmn12%nqwait .ne. 0) go to 600
      cmn12%nqwait = 2
      cmn12%etaqm1 = zero
      if (cmn12%nq .eq. 1) go to 570
! compute ratio of new h to current h at the current order less one. ---
      ddn = dvnorm (cmn12%n, yh(1,cmn12%l), ewt)/cmn12%tq(1)
      cmn12%etaqm1 = one/((bias1*ddn)**(one/(flotl - one)) + addon)
 570  etaqp1 = zero
      if (cmn12%l .eq. cmn12%lmax) go to 580
! compute ratio of new h to current h at current order plus one. -------
      cnquot = (cmn12%tq(5)/cmn12%conp)*(cmn12%h/cmn12%tau(2))**cmn12%l
      do 575 i = 1, cmn12%n
 575    savf(i) = acor(i) - cnquot*yh(i,cmn12%lmax)
      dup = dvnorm (cmn12%n, savf, ewt)/cmn12%tq(3)
      etaqp1 = one/((bias3*dup)**(one/(flotl + one)) + addon)
 580  if (cmn12%etaq .ge. etaqp1) go to 590
      if (etaqp1 .gt. cmn12%etaqm1) go to 620
      go to 610
 590  if (cmn12%etaq .lt. cmn12%etaqm1) go to 610
 600  cmn12%eta = cmn12%etaq
      cmn12%newq = cmn12%nq
      go to 630
 610  cmn12%eta = cmn12%etaqm1
      cmn12%newq = cmn12%nq - 1
      go to 630
 620  cmn12%eta = etaqp1
      cmn12%newq = cmn12%nq + 1
      call dcopy (cmn12%n, acor, 1, yh(1,cmn12%lmax), 1)
! test tentative new h against thresh, etamax, and hmxi, then exit. ----
 630  if (cmn12%eta .lt. thresh .or. cmn12%etamax .eq. one) go to 640
      cmn12%eta = min(cmn12%eta,cmn12%etamax)
      cmn12%eta = cmn12%eta/max(one,abs(cmn12%h)*cmn12%hmxi*cmn12%eta)
      cmn12%newh = 1
      cmn12%hnew = cmn12%h*cmn12%eta
      go to 690
 640  cmn12%newq = cmn12%nq
      cmn12%newh = 0
      cmn12%eta = one
      cmn12%hnew = cmn12%h
      go to 690
!-----------------------------------------------------------------------
! all returns are made through this section.
! on a successful return, etamax is reset and acor is scaled.
!-----------------------------------------------------------------------
 660  cmn12%kflag = -1
      go to 720
 670  cmn12%kflag = -2
      go to 720
 680  if (nflag .eq. -2) cmn12%kflag = -3
      if (nflag .eq. -3) cmn12%kflag = -4
      go to 720
 690  cmn12%etamax = etamx3
      if (cmn12%nst .le. 10) cmn12%etamax = etamx2
 700  r = one/cmn12%tq(2)
      call dscal (cmn12%n, r, acor, 1)
 720  cmn12%jstart = 1
      return
!----------------------- end of subroutine dvstep ----------------------
      end subroutine dvstep
!deck dvset
      subroutine dvset( cmn12 )
      implicit none
      type( dvode_cmn_vars ) :: cmn12
!-----------------------------------------------------------------------
! call sequence communication: none
! common block variables accessed:
!     /dvod_cmn01/ -- el(13), h, tau(13), tq(5), l(= nq + 1),
!                 meth, nq, nqwait
!
! subroutines called by dvset: none
! function routines called by dvset: none
!-----------------------------------------------------------------------
! dvset is called by dvstep and sets coefficients for use there.
!
! for each order nq, the coefficients in el are calculated by use of
!  the generating polynomial lambda(x), with coefficients el(i).
!      lambda(x) = el(1) + el(2)*x + ... + el(nq+1)*(x**nq).
! for the backward differentiation formulas,
!                                     nq-1
!      lambda(x) = (1 + x/xi*(nq)) * product (1 + x/xi(i) ) .
!                                     i = 1
! for the adams formulas,
!                              nq-1
!      (d/dx) lambda(x) = c * product (1 + x/xi(i) ) ,
!                              i = 1
!      lambda(-1) = 0,    lambda(0) = 1,
! where c is a normalization constant.
! in both cases, xi(i) is defined by
!      h*xi(i) = t sub n  -  t sub (n-i)
!              = h + tau(1) + tau(2) + ... tau(i-1).
!
!
! in addition to variables described previously, communication
! with dvset uses the following:
!   tau    = a vector of length 13 containing the past nq values
!            of h.
!   el     = a vector of length 13 in which vset stores the
!            coefficients for the corrector formula.
!   tq     = a vector of length 5 in which vset stores constants
!            used for the convergence test, the error test, and the
!            selection of h at a new order.
!   meth   = the basic method indicator.
!   nq     = the current order.
!   l      = nq + 1, the length of the vector stored in el, and
!            the number of columns of the yh array being used.
!   nqwait = a counter controlling the frequency of order changes.
!            an order change is about to be considered if nqwait = 1.
!-----------------------------------------------------------------------
!
! type declarations for labeled common block dvod_cmn01 --------------------
!
!      double precision acnrm, ccmxj, conp, crate, drc, el,   &
!           eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!           rc, rl1, tau, tq, tn, uround
!      integer icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!              l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!              locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!              n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!              nslp, nyh
!
! type declarations for local variables --------------------------------
!
      double precision ahatn0, alph0, cnqm1, cortes, csum, elp, em,   &
           em0, floti, flotl, flotnq, hsum, one, rxi, rxis, s, six,   &
           t1, t2, t3, t4, t5, t6, two, xi, zero
      integer i, iback, j, jp1, nqm1, nqm2
!
      dimension em(13)
!-----------------------------------------------------------------------
! the following fortran-77 declaration is to cause the values of the
! listed (local) variables to be saved between calls to this integrator.
!-----------------------------------------------------------------------
      save cortes, one, six, two, zero
!
!      common /dvod_cmn01/ acnrm, ccmxj, conp, crate, drc, el(13),   &
!                      eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!                      rc, rl1, tau(13), tq(5), tn, uround,   &
!                      icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!                      l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!                      locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!                      n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!                      nslp, nyh
!
      data cortes /0.1d0/
      data one  /1.0d0/, six /6.0d0/, two /2.0d0/, zero /0.0d0/
!
      flotl = real(cmn12%l)
      nqm1 = cmn12%nq - 1
      nqm2 = cmn12%nq - 2
      go to (100, 200), cmn12%meth
!
! set coefficients for adams methods. ----------------------------------
 100  if (cmn12%nq .ne. 1) go to 110
      cmn12%el(1) = one
      cmn12%el(2) = one
      cmn12%tq(1) = one
      cmn12%tq(2) = two
      cmn12%tq(3) = six*cmn12%tq(2)
      cmn12%tq(5) = one
      go to 300
 110  hsum = cmn12%h
      em(1) = one
      flotnq = flotl - one
      do 115 i = 2, cmn12%l
 115    em(i) = zero
      do 150 j = 1, nqm1
        if ((j .ne. nqm1) .or. (cmn12%nqwait .ne. 1)) go to 130
        s = one
        csum = zero
        do 120 i = 1, nqm1
          csum = csum + s*em(i)/real(i+1)
 120      s = -s
        cmn12%tq(1) = em(nqm1)/(flotnq*csum)
 130    rxi = cmn12%h/hsum
        do 140 iback = 1, j
          i = (j + 2) - iback
 140      em(i) = em(i) + em(i-1)*rxi
        hsum = hsum + cmn12%tau(j)
 150    continue
! compute integral from -1 to 0 of polynomial and of x times it. -------
      s = one
      em0 = zero
      csum = zero
      do 160 i = 1, cmn12%nq
        floti = real(i)
        em0 = em0 + s*em(i)/floti
        csum = csum + s*em(i)/(floti+one)
 160    s = -s
! in el, form coefficients of normalized integrated polynomial. --------
      s = one/em0
      cmn12%el(1) = one
      do 170 i = 1, cmn12%nq
 170    cmn12%el(i+1) = s*em(i)/real(i)
      xi = hsum/cmn12%h
      cmn12%tq(2) = xi*em0/csum
      cmn12%tq(5) = xi/cmn12%el(cmn12%l)
      if (cmn12%nqwait .ne. 1) go to 300
! for higher order control constant, multiply polynomial by 1+x/xi(q). -
      rxi = one/xi
      do 180 iback = 1, cmn12%nq
        i = (cmn12%l + 1) - iback
 180    em(i) = em(i) + em(i-1)*rxi
! compute integral of polynomial. --------------------------------------
      s = one
      csum = zero
      do 190 i = 1, cmn12%l
        csum = csum + s*em(i)/real(i+1)
 190    s = -s
      cmn12%tq(3) = flotl*em0/csum
      go to 300
!
! set coefficients for bdf methods. ------------------------------------
 200  do 210 i = 3, cmn12%l
 210    cmn12%el(i) = zero
      cmn12%el(1) = one
      cmn12%el(2) = one
      alph0 = -one
      ahatn0 = -one
      hsum = cmn12%h
      rxi = one
      rxis = one
      if (cmn12%nq .eq. 1) go to 240
      do 230 j = 1, nqm2
! in el, construct coefficients of (1+x/xi(1))*...*(1+x/xi(j+1)). ------
        hsum = hsum + cmn12%tau(j)
        rxi = cmn12%h/hsum
        jp1 = j + 1
        alph0 = alph0 - one/real(jp1)
        do 220 iback = 1, jp1
          i = (j + 3) - iback
 220      cmn12%el(i) = cmn12%el(i) + cmn12%el(i-1)*rxi
 230    continue
      alph0 = alph0 - one/real(cmn12%nq)
      rxis = -cmn12%el(2) - alph0
      hsum = hsum + cmn12%tau(nqm1)
      rxi = cmn12%h/hsum
      ahatn0 = -cmn12%el(2) - rxi
      do 235 iback = 1, cmn12%nq
        i = (cmn12%nq + 2) - iback
 235    cmn12%el(i) = cmn12%el(i) + cmn12%el(i-1)*rxis
 240  t1 = one - ahatn0 + alph0
      t2 = one + real(cmn12%nq)*t1
      cmn12%tq(2) = abs(alph0*t2/t1)
      cmn12%tq(5) = abs(t2/(cmn12%el(cmn12%l)*rxi/rxis))
      if (cmn12%nqwait .ne. 1) go to 300
      cnqm1 = rxis/cmn12%el(cmn12%l)
      t3 = alph0 + one/real(cmn12%nq)
      t4 = ahatn0 + rxi
      elp = t3/(one - t4 + t3)
      cmn12%tq(1) = abs(elp/cnqm1)
      hsum = hsum + cmn12%tau(cmn12%nq)
      rxi = cmn12%h/hsum
      t5 = alph0 - one/real(cmn12%nq+1)
      t6 = ahatn0 - rxi
      elp = t2/(one - t6 + t5)
      cmn12%tq(3) = abs(elp*rxi*(flotl + one)*t5)
 300  cmn12%tq(4) = cortes*cmn12%tq(2)
      return
!----------------------- end of subroutine dvset -----------------------
      end subroutine dvset
!deck dvjust
      subroutine dvjust (yh, ldyh, iord, cmn12)
      implicit none
      double precision yh
      integer ldyh, iord
      dimension yh(ldyh,*)
      type( dvode_cmn_vars ) :: cmn12
!-----------------------------------------------------------------------
! call sequence input -- yh, ldyh, iord
! call sequence output -- yh
! common block input -- nq, meth, lmax, hscal, tau(13), n
! common block variables accessed:
!     /dvod_cmn01/ -- hscal, tau(13), lmax, meth, n, nq,
!
! subroutines called by dvjust: daxpy
! function routines called by dvjust: none
!-----------------------------------------------------------------------
! this subroutine adjusts the yh array on reduction of order,
! and also when the order is increased for the stiff option (meth = 2).
! communication with dvjust uses the following:
! iord  = an integer flag used when meth = 2 to indicate an order
!         increase (iord = +1) or an order decrease (iord = -1).
! hscal = step size h used in scaling of nordsieck array yh.
!         (if iord = +1, dvjust assumes that hscal = tau(1).)
! see references 1 and 2 for details.
!-----------------------------------------------------------------------
!
! type declarations for labeled common block dvod_cmn01 --------------------
!
!      double precision acnrm, ccmxj, conp, crate, drc, el,   &
!           eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!           rc, rl1, tau, tq, tn, uround
!      integer icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!              l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!              locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!              n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!              nslp, nyh
!
! type declarations for local variables --------------------------------
!
      double precision alph0, alph1, hsum, one, prod, t1, xi,xiold, zero
      integer i, iback, j, jp1, lp1, nqm1, nqm2, nqp1
!-----------------------------------------------------------------------
! the following fortran-77 declaration is to cause the values of the
! listed (local) variables to be saved between calls to this integrator.
!-----------------------------------------------------------------------
      save one, zero
!
!      common /dvod_cmn01/ acnrm, ccmxj, conp, crate, drc, el(13),   &
!                      eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!                      rc, rl1, tau(13), tq(5), tn, uround,   &
!                      icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!                      l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!                      locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!                      n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!                      nslp, nyh
!
      data one /1.0d0/, zero /0.0d0/
!
      if ((cmn12%nq .eq. 2) .and. (iord .ne. 1)) return
      nqm1 = cmn12%nq - 1
      nqm2 = cmn12%nq - 2
      go to (100, 200), cmn12%meth
!-----------------------------------------------------------------------
! nonstiff option...
! check to see if the order is being increased or decreased.
!-----------------------------------------------------------------------
 100  continue
      if (iord .eq. 1) go to 180
! order decrease. ------------------------------------------------------
      do 110 j = 1, cmn12%lmax
 110    cmn12%el(j) = zero
      cmn12%el(2) = one
      hsum = zero
      do 130 j = 1, nqm2
! construct coefficients of x*(x+xi(1))*...*(x+xi(j)). -----------------
        hsum = hsum + cmn12%tau(j)
        xi = hsum/cmn12%hscal
        jp1 = j + 1
        do 120 iback = 1, jp1
          i = (j + 3) - iback
 120      cmn12%el(i) = cmn12%el(i)*xi + cmn12%el(i-1)
 130    continue
! construct coefficients of integrated polynomial. ---------------------
      do 140 j = 2, nqm1
 140    cmn12%el(j+1) = real(cmn12%nq)*cmn12%el(j)/real(j)
! subtract correction terms from yh array. -----------------------------
      do 170 j = 3, cmn12%nq
        do 160 i = 1, cmn12%n
 160      yh(i,j) = yh(i,j) - yh(i,cmn12%l)*cmn12%el(j)
 170    continue
      return
! order increase. ------------------------------------------------------
! zero out next column in yh array. ------------------------------------
 180  continue
      lp1 = cmn12%l + 1
      do 190 i = 1, cmn12%n
 190    yh(i,lp1) = zero
      return
!-----------------------------------------------------------------------
! stiff option...
! check to see if the order is being increased or decreased.
!-----------------------------------------------------------------------
 200  continue
      if (iord .eq. 1) go to 300
! order decrease. ------------------------------------------------------
      do 210 j = 1, cmn12%lmax
 210    cmn12%el(j) = zero
      cmn12%el(3) = one
      hsum = zero
      do 230 j = 1,nqm2
! construct coefficients of x*x*(x+xi(1))*...*(x+xi(j)). ---------------
        hsum = hsum + cmn12%tau(j)
        xi = hsum/cmn12%hscal
        jp1 = j + 1
        do 220 iback = 1, jp1
          i = (j + 4) - iback
 220      cmn12%el(i) = cmn12%el(i)*xi + cmn12%el(i-1)
 230    continue
! subtract correction terms from yh array. -----------------------------
      do 250 j = 3,cmn12%nq
        do 240 i = 1, cmn12%n
 240      yh(i,j) = yh(i,j) - yh(i,cmn12%l)*cmn12%el(j)
 250    continue
      return
! order increase. ------------------------------------------------------
 300  do 310 j = 1, cmn12%lmax
 310    cmn12%el(j) = zero
      cmn12%el(3) = one
      alph0 = -one
      alph1 = one
      prod = one
      xiold = one
      hsum = cmn12%hscal
      if (cmn12%nq .eq. 1) go to 340
      do 330 j = 1, nqm1
! construct coefficients of x*x*(x+xi(1))*...*(x+xi(j)). ---------------
        jp1 = j + 1
        hsum = hsum + cmn12%tau(jp1)
        xi = hsum/cmn12%hscal
        prod = prod*xi
        alph0 = alph0 - one/real(jp1)
        alph1 = alph1 + one/xi
        do 320 iback = 1, jp1
          i = (j + 4) - iback
 320      cmn12%el(i) = cmn12%el(i)*xiold + cmn12%el(i-1)
        xiold = xi
 330    continue
 340  continue
      t1 = (-alph0 - alph1)/prod
! load column l + 1 in yh array. ---------------------------------------
      lp1 = cmn12%l + 1
      do 350 i = 1, cmn12%n
 350    yh(i,lp1) = t1*yh(i,cmn12%lmax)
! add correction terms to yh array. ------------------------------------
      nqp1 = cmn12%nq + 1
      do 370 j = 3, nqp1
        call daxpy (cmn12%n, cmn12%el(j), yh(1,lp1), 1, yh(1,j), 1 )
 370  continue
      return
!----------------------- end of subroutine dvjust ----------------------
      end subroutine dvjust
!deck dvnlsd
      subroutine dvnlsd (y, yh, ldyh, vsav, savf, ewt, acor, iwm, wm,   &
                       f, jac, pdum, nflag, rpar, ipar, cmn12)
      implicit none
      external f, jac, pdum
      double precision y, yh, vsav, savf, ewt, acor, wm, rpar
      integer ldyh, iwm, nflag, ipar
      dimension y(*), yh(ldyh,*), vsav(*), savf(*), ewt(*), acor(*),   &
                iwm(*), wm(*), rpar(*), ipar(*)
      type( dvode_cmn_vars ) :: cmn12
!-----------------------------------------------------------------------
! call sequence input -- y, yh, ldyh, savf, ewt, acor, iwm, wm,
!                        f, jac, nflag, rpar, ipar
! call sequence output -- yh, acor, wm, iwm, nflag
! common block variables accessed:
!     /dvod_cmn01/ acnrm, crate, drc, h, rc, rl1, tq(5), tn, icf,
!                jcur, meth, miter, n, nslp
!     /dvod_cmn02/ hu, ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
! subroutines called by dvnlsd: f, daxpy, dcopy, dscal, dvjac, dvsol
! function routines called by dvnlsd: dvnorm
!-----------------------------------------------------------------------
! subroutine dvnlsd is a nonlinear system solver, which uses functional
! iteration or a chord (modified newton) method.  for the chord method
! direct linear algebraic system solvers are used.  subroutine dvnlsd
! then handles the corrector phase of this integration package.
!
! communication with dvnlsd is done with the following variables. (for
! more details, please see the comments in the driver subroutine.)
!
! y          = the dependent variable, a vector of length n, input.
! yh         = the nordsieck (taylor) array, ldyh by lmax, input
!              and output.  on input, it contains predicted values.
! ldyh       = a constant .ge. n, the first dimension of yh, input.
! vsav       = unused work array.
! savf       = a work array of length n.
! ewt        = an error weight vector of length n, input.
! acor       = a work array of length n, used for the accumulated
!              corrections to the predicted y vector.
! wm,iwm     = real and integer work arrays associated with matrix
!              operations in chord iteration (miter .ne. 0).
! f          = dummy name for user supplied routine for f.
! jac        = dummy name for user supplied jacobian routine.
! pdum       = unused dummy subroutine name.  included for uniformity
!              over collection of integrators.
! nflag      = input/output flag, with values and meanings as follows:
!              input
!                  0 first call for this time step.
!                 -1 convergence failure in previous call to dvnlsd.
!                 -2 error test failure in dvstep.
!              output
!                  0 successful completion of nonlinear solver.
!                 -1 convergence failure or singular matrix.
!                 -2 unrecoverable error in matrix preprocessing
!                    (cannot occur here).
!                 -3 unrecoverable error in solution (cannot occur
!                    here).
! rpar, ipar = dummy names for user's real and integer work arrays.
!
! ipup       = own variable flag with values and meanings as follows:
!              0,            do not update the newton matrix.
!              miter .ne. 0, update newton matrix, because it is the
!                            initial step, order was changed, the error
!                            test failed, or an update is indicated by
!                            the scalar rc or step counter nst.
!
! for more details, see comments in driver subroutine.
!-----------------------------------------------------------------------
!! type declarations for labeled common block dvod_cmn01 --------------------
!!
!      double precision acnrm, ccmxj, conp, crate, drc, el,   &
!           eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!           rc, rl1, tau, tq, tn, uround
!      integer icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!              l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!              locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!              n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!              nslp, nyh
!!
!! type declarations for labeled common block dvod_cmn02 --------------------
!!
!      double precision hu
!      integer ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
! type declarations for local variables --------------------------------
!
      double precision ccmax, crdown, cscale, dcon, del, delp, one,   &
           rdiv, two, zero
      integer i, ierpj, iersl, m, maxcor, msbp
!
! type declaration for function subroutines called ---------------------
!
! 27-oct-2005 rce - do not declare functions that are in the module
!     double precision dvnorm
!-----------------------------------------------------------------------
! the following fortran-77 declaration is to cause the values of the
! listed (local) variables to be saved between calls to this integrator.
!-----------------------------------------------------------------------
      save ccmax, crdown, maxcor, msbp, rdiv, one, two, zero
!
!      common /dvod_cmn01/ acnrm, ccmxj, conp, crate, drc, el(13),   &
!                      eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!                      rc, rl1, tau(13), tq(5), tn, uround,   &
!                      icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!                      l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!                      locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!                      n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!                      nslp, nyh
!      common /dvod_cmn02/ hu, ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
      data ccmax /0.3d0/, crdown /0.3d0/, maxcor /3/, msbp /20/,   &
           rdiv  /2.0d0/
      data one /1.0d0/, two /2.0d0/, zero /0.0d0/
!-----------------------------------------------------------------------
! on the first step, on a change of method order, or after a
! nonlinear convergence failure with nflag = -2, set ipup = miter
! to force a jacobian update when miter .ne. 0.
!-----------------------------------------------------------------------
      if (cmn12%jstart .eq. 0) cmn12%nslp = 0
      if (nflag .eq. 0) cmn12%icf = 0
      if (nflag .eq. -2) cmn12%ipup = cmn12%miter
      if ( (cmn12%jstart .eq. 0) .or. (cmn12%jstart .eq. -1) ) cmn12%ipup = cmn12%miter
! if this is functional iteration, set crate .eq. 1 and drop to 220
      if (cmn12%miter .eq. 0) then
        cmn12%crate = one
        go to 220
      endif
!-----------------------------------------------------------------------
! rc is the ratio of new to old values of the coefficient h/el(2)=h/l1.
! when rc differs from 1 by more than ccmax, ipup is set to miter
! to force dvjac to be called, if a jacobian is involved.
! in any case, dvjac is called at least every msbp steps.
!-----------------------------------------------------------------------
      cmn12%drc = abs(cmn12%rc-one)
      if (cmn12%drc .gt. ccmax .or. cmn12%nst .ge. cmn12%nslp+msbp) cmn12%ipup = cmn12%miter
!-----------------------------------------------------------------------
! up to maxcor corrector iterations are taken.  a convergence test is
! made on the r.m.s. norm of each correction, weighted by the error
! weight vector ewt.  the sum of the corrections is accumulated in the
! vector acor(i).  the yh array is not altered in the corrector loop.
!-----------------------------------------------------------------------
 220  m = 0
      delp = zero
      call dcopy (cmn12%n, yh(1,1), 1, y, 1 )
      call f (cmn12%n, cmn12%tn, y, savf, rpar, ipar)
      cmn12%nfe = cmn12%nfe + 1
      if (cmn12%ipup .le. 0) go to 250
!-----------------------------------------------------------------------
! if indicated, the matrix p = i - h*rl1*j is reevaluated and
! preprocessed before starting the corrector iteration.  ipup is set
! to 0 as an indicator that this has been done.
!-----------------------------------------------------------------------
      call dvjac (y, yh, ldyh, ewt, acor, savf, wm, iwm, f, jac, ierpj,   &
                 rpar, ipar, cmn12)
      cmn12%ipup = 0
      cmn12%rc = one
      cmn12%drc = zero
      cmn12%crate = one
      cmn12%nslp = cmn12%nst
! if matrix is singular, take error return to force cut in step size. --
      if (ierpj .ne. 0) go to 430
 250  do 260 i = 1,cmn12%n
 260    acor(i) = zero
! this is a looping point for the corrector iteration. -----------------
 270  if (cmn12%miter .ne. 0) go to 350
!-----------------------------------------------------------------------
! in the case of functional iteration, update y directly from
! the result of the last function evaluation.
!-----------------------------------------------------------------------
      do 280 i = 1,cmn12%n
 280    savf(i) = cmn12%rl1*(cmn12%h*savf(i) - yh(i,2))
      do 290 i = 1,cmn12%n
 290    y(i) = savf(i) - acor(i)
      del = dvnorm (cmn12%n, y, ewt)
      do 300 i = 1,cmn12%n
 300    y(i) = yh(i,1) + savf(i)
      call dcopy (cmn12%n, savf, 1, acor, 1)
      go to 400
!-----------------------------------------------------------------------
! in the case of the chord method, compute the corrector error,
! and solve the linear system with that as right-hand side and
! p as coefficient matrix.  the correction is scaled by the factor
! 2/(1+rc) to account for changes in h*rl1 since the last dvjac call.
!-----------------------------------------------------------------------
 350  do 360 i = 1,cmn12%n
 360    y(i) = (cmn12%rl1*cmn12%h)*savf(i) - (cmn12%rl1*yh(i,2) + acor(i))
      call dvsol (wm, iwm, y, iersl, cmn12)
      cmn12%nni = cmn12%nni + 1
      if (iersl .gt. 0) go to 410
      if (cmn12%meth .eq. 2 .and. cmn12%rc .ne. one) then
        cscale = two/(one + cmn12%rc)
        call dscal (cmn12%n, cscale, y, 1)
      endif
      del = dvnorm (cmn12%n, y, ewt)
      call daxpy (cmn12%n, one, y, 1, acor, 1)
      do 380 i = 1,cmn12%n
 380    y(i) = yh(i,1) + acor(i)
!-----------------------------------------------------------------------
! test for convergence.  if m .gt. 0, an estimate of the convergence
! rate constant is stored in crate, and this is used in the test.
!-----------------------------------------------------------------------
 400  if (m .ne. 0) cmn12%crate = max(crdown*cmn12%crate,del/delp)
      dcon = del*min(one,cmn12%crate)/cmn12%tq(4)
      if (dcon .le. one) go to 450
      m = m + 1
      if (m .eq. maxcor) go to 410
      if (m .ge. 2 .and. del .gt. rdiv*delp) go to 410
      delp = del
      call f (cmn12%n, cmn12%tn, y, savf, rpar, ipar)
      cmn12%nfe = cmn12%nfe + 1
      go to 270
!
 410  if (cmn12%miter .eq. 0 .or. cmn12%jcur .eq. 1) go to 430
      cmn12%icf = 1
      cmn12%ipup = cmn12%miter
      go to 220
!
 430  continue
      nflag = -1
      cmn12%icf = 2
      cmn12%ipup = cmn12%miter
      return
!
! return for successful step. ------------------------------------------
 450  nflag = 0
      cmn12%jcur = 0
      cmn12%icf = 0
      if (m .eq. 0) cmn12%acnrm = del
      if (m .gt. 0) cmn12%acnrm = dvnorm (cmn12%n, acor, ewt)
      return
!----------------------- end of subroutine dvnlsd ----------------------
      end subroutine dvnlsd
!deck dvjac
      subroutine dvjac (y, yh, ldyh, ewt, ftem, savf, wm, iwm, f, jac,   &
                       ierpj, rpar, ipar, cmn12)
      implicit none
      external f, jac
      double precision y, yh, ewt, ftem, savf, wm, rpar
      integer ldyh, iwm, ierpj, ipar
      dimension y(*), yh(ldyh,*), ewt(*), ftem(*), savf(*),   &
         wm(*), iwm(*), rpar(*), ipar(*)
      type( dvode_cmn_vars ) :: cmn12
!-----------------------------------------------------------------------
! call sequence input -- y, yh, ldyh, ewt, ftem, savf, wm, iwm,
!                        f, jac, rpar, ipar
! call sequence output -- wm, iwm, ierpj
! common block variables accessed:
!     /dvod_cmn01/  ccmxj, drc, h, rl1, tn, uround, icf, jcur, locjs,
!               miter, msbj, n, nslj
!     /dvod_cmn02/  nfe, nst, nje, nlu
!
! subroutines called by dvjac: f, jac, dacopy, dcopy, dgbfa, dgefa,
!                              dscal
! function routines called by dvjac: dvnorm
!-----------------------------------------------------------------------
! dvjac is called by dvnlsd to compute and process the matrix
! p = i - h*rl1*j , where j is an approximation to the jacobian.
! here j is computed by the user-supplied routine jac if
! miter = 1 or 4, or by finite differencing if miter = 2, 3, or 5.
! if miter = 3, a diagonal approximation to j is used.
! if jsv = -1, j is computed from scratch in all cases.
! if jsv = 1 and miter = 1, 2, 4, or 5, and if the saved value of j is
! considered acceptable, then p is constructed from the saved j.
! j is stored in wm and replaced by p.  if miter .ne. 3, p is then
! subjected to lu decomposition in preparation for later solution
! of linear systems with p as coefficient matrix. this is done
! by dgefa if miter = 1 or 2, and by dgbfa if miter = 4 or 5.
!
! communication with dvjac is done with the following variables.  (for
! more details, please see the comments in the driver subroutine.)
! y          = vector containing predicted values on entry.
! yh         = the nordsieck array, an ldyh by lmax array, input.
! ldyh       = a constant .ge. n, the first dimension of yh, input.
! ewt        = an error weight vector of length n.
! savf       = array containing f evaluated at predicted y, input.
! wm         = real work space for matrices.  in the output, it contains
!              the inverse diagonal matrix if miter = 3 and the lu
!              decomposition of p if miter is 1, 2 , 4, or 5.
!              storage of matrix elements starts at wm(3).
!              storage of the saved jacobian starts at wm(locjs).
!              wm also contains the following matrix-related data:
!              wm(1) = sqrt(uround), used in numerical jacobian step.
!              wm(2) = h*rl1, saved for later use if miter = 3.
! iwm        = integer work space containing pivot information,
!              starting at iwm(31), if miter is 1, 2, 4, or 5.
!              iwm also contains band parameters ml = iwm(1) and
!              mu = iwm(2) if miter is 4 or 5.
! f          = dummy name for the user supplied subroutine for f.
! jac        = dummy name for the user supplied jacobian subroutine.
! rpar, ipar = dummy names for user's real and integer work arrays.
! rl1        = 1/el(2) (input).
! ierpj      = output error flag,  = 0 if no trouble, 1 if the p
!              matrix is found to be singular.
! jcur       = output flag to indicate whether the jacobian matrix
!              (or approximation) is now current.
!              jcur = 0 means j is not current.
!              jcur = 1 means j is current.
!-----------------------------------------------------------------------
!!
!! type declarations for labeled common block dvod_cmn01 --------------------
!!
!      double precision acnrm, ccmxj, conp, crate, drc, el,   &
!           eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!           rc, rl1, tau, tq, tn, uround
!      integer icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!              l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!              locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!              n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!              nslp, nyh
!!
!! type declarations for labeled common block dvod_cmn02 --------------------
!!
!      double precision hu
!      integer ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
! type declarations for local variables --------------------------------
!
      double precision con, di, fac, hrl1, one, pt1, r, r0, srur, thou,   &
           yi, yj, yjj, zero
      integer i, i1, i2, ier, ii, j, j1, jj, jok, lenp, mba, mband,   &
              meb1, meband, ml, ml3, mu, np1
!
! type declaration for function subroutines called ---------------------
!
! 27-oct-2005 rce - do not declare functions that are in the module
!     double precision dvnorm
!-----------------------------------------------------------------------
! the following fortran-77 declaration is to cause the values of the
! listed (local) variables to be saved between calls to this subroutine.
!-----------------------------------------------------------------------
      save one, pt1, thou, zero
!-----------------------------------------------------------------------
!      common /dvod_cmn01/ acnrm, ccmxj, conp, crate, drc, el(13),   &
!                      eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!                      rc, rl1, tau(13), tq(5), tn, uround,   &
!                      icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!                      l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!                      locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!                      n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!                      nslp, nyh
!      common /dvod_cmn02/ hu, ncfn, netf, nfe, nje, nlu, nni, nqu, nst
!
      data one /1.0d0/, thou /1000.0d0/, zero /0.0d0/, pt1 /0.1d0/
!
      ierpj = 0
      hrl1 = cmn12%h*cmn12%rl1
! see whether j should be evaluated (jok = -1) or not (jok = 1). -------
      jok = cmn12%jsv
      if (cmn12%jsv .eq. 1) then
        if (cmn12%nst .eq. 0 .or. cmn12%nst .gt. cmn12%nslj+cmn12%msbj) jok = -1
        if (cmn12%icf .eq. 1 .and. cmn12%drc .lt. cmn12%ccmxj) jok = -1
        if (cmn12%icf .eq. 2) jok = -1
      endif
! end of setting jok. --------------------------------------------------
!
      if (jok .eq. -1 .and. cmn12%miter .eq. 1) then
! if jok = -1 and miter = 1, call jac to evaluate jacobian. ------------
      cmn12%nje = cmn12%nje + 1
      cmn12%nslj = cmn12%nst
      cmn12%jcur = 1
      lenp = cmn12%n*cmn12%n
      do 110 i = 1,lenp
 110    wm(i+2) = zero
      call jac (cmn12%n, cmn12%tn, y, 0, 0, wm(3), cmn12%n, rpar, ipar)
      if (cmn12%jsv .eq. 1) call dcopy (lenp, wm(3), 1, wm(cmn12%locjs), 1)
      endif
!
      if (jok .eq. -1 .and. cmn12%miter .eq. 2) then
! if miter = 2, make n calls to f to approximate the jacobian. ---------
      cmn12%nje = cmn12%nje + 1
      cmn12%nslj = cmn12%nst
      cmn12%jcur = 1
      fac = dvnorm (cmn12%n, savf, ewt)
      r0 = thou*abs(cmn12%h)*cmn12%uround*real(cmn12%n)*fac
      if (r0 .eq. zero) r0 = one
      srur = wm(1)
      j1 = 2
      do 230 j = 1,cmn12%n
        yj = y(j)
        r = max(srur*abs(yj),r0/ewt(j))
        y(j) = y(j) + r
        fac = one/r
        call f (cmn12%n, cmn12%tn, y, ftem, rpar, ipar)
        do 220 i = 1,cmn12%n
 220      wm(i+j1) = (ftem(i) - savf(i))*fac
        y(j) = yj
        j1 = j1 + cmn12%n
 230    continue
      cmn12%nfe = cmn12%nfe + cmn12%n
      lenp = cmn12%n*cmn12%n
      if (cmn12%jsv .eq. 1) call dcopy (lenp, wm(3), 1, wm(cmn12%locjs), 1)
      endif
!
      if (jok .eq. 1 .and. (cmn12%miter .eq. 1 .or. cmn12%miter .eq. 2)) then
      cmn12%jcur = 0
      lenp = cmn12%n*cmn12%n
      call dcopy (lenp, wm(cmn12%locjs), 1, wm(3), 1)
      endif
!
      if (cmn12%miter .eq. 1 .or. cmn12%miter .eq. 2) then
! multiply jacobian by scalar, add identity, and do lu decomposition. --
      con = -hrl1
      call dscal (lenp, con, wm(3), 1)
      j = 3
      np1 = cmn12%n + 1
      do 250 i = 1,cmn12%n
        wm(j) = wm(j) + one
 250    j = j + np1
      cmn12%nlu = cmn12%nlu + 1
      call dgefa (wm(3), cmn12%n, cmn12%n, iwm(31), ier)
      if (ier .ne. 0) ierpj = 1
      return
      endif
! end of code block for miter = 1 or 2. --------------------------------
!
      if (cmn12%miter .eq. 3) then
! if miter = 3, construct a diagonal approximation to j and p. ---------
      cmn12%nje = cmn12%nje + 1
      cmn12%jcur = 1
      wm(2) = hrl1
      r = cmn12%rl1*pt1
      do 310 i = 1,cmn12%n
 310    y(i) = y(i) + r*(cmn12%h*savf(i) - yh(i,2))
      call f (cmn12%n, cmn12%tn, y, wm(3), rpar, ipar)
      cmn12%nfe = cmn12%nfe + 1
      do 320 i = 1,cmn12%n
        r0 = cmn12%h*savf(i) - yh(i,2)
        di = pt1*r0 - cmn12%h*(wm(i+2) - savf(i))
        wm(i+2) = one
        if (abs(r0) .lt. cmn12%uround/ewt(i)) go to 320
        if (abs(di) .eq. zero) go to 330
        wm(i+2) = pt1*r0/di
 320    continue
      return
 330  ierpj = 1
      return
      endif
! end of code block for miter = 3. -------------------------------------
!
! set constants for miter = 4 or 5. ------------------------------------
      ml = iwm(1)
      mu = iwm(2)
      ml3 = ml + 3
      mband = ml + mu + 1
      meband = mband + ml
      lenp = meband*cmn12%n
!
      if (jok .eq. -1 .and. cmn12%miter .eq. 4) then
! if jok = -1 and miter = 4, call jac to evaluate jacobian. ------------
      cmn12%nje = cmn12%nje + 1
      cmn12%nslj = cmn12%nst
      cmn12%jcur = 1
      do 410 i = 1,lenp
 410    wm(i+2) = zero
      call jac (cmn12%n, cmn12%tn, y, ml, mu, wm(ml3), meband, rpar, ipar)
      if (cmn12%jsv .eq. 1)   &
         call dacopy (mband, cmn12%n, wm(ml3), meband, wm(cmn12%locjs), mband)
      endif
!
      if (jok .eq. -1 .and. cmn12%miter .eq. 5) then
! if miter = 5, make ml+mu+1 calls to f to approximate the jacobian. ---
      cmn12%nje = cmn12%nje + 1
      cmn12%nslj = cmn12%nst
      cmn12%jcur = 1
      mba = min(mband,cmn12%n)
      meb1 = meband - 1
      srur = wm(1)
      fac = dvnorm (cmn12%n, savf, ewt)
      r0 = thou*abs(cmn12%h)*cmn12%uround*real(cmn12%n)*fac
      if (r0 .eq. zero) r0 = one
      do 560 j = 1,mba
        do 530 i = j,cmn12%n,mband
          yi = y(i)
          r = max(srur*abs(yi),r0/ewt(i))
 530      y(i) = y(i) + r
        call f (cmn12%n, cmn12%tn, y, ftem, rpar, ipar)
        do 550 jj = j,cmn12%n,mband
          y(jj) = yh(jj,1)
          yjj = y(jj)
          r = max(srur*abs(yjj),r0/ewt(jj))
          fac = one/r
          i1 = max(jj-mu,1)
          i2 = min(jj+ml,cmn12%n)
          ii = jj*meb1 - ml + 2
          do 540 i = i1,i2
 540        wm(ii+i) = (ftem(i) - savf(i))*fac
 550      continue
 560    continue
      cmn12%nfe = cmn12%nfe + mba
      if (cmn12%jsv .eq. 1)   &
         call dacopy (mband, cmn12%n, wm(ml3), meband, wm(cmn12%locjs), mband)
      endif
!
      if (jok .eq. 1) then
      cmn12%jcur = 0
      call dacopy (mband, cmn12%n, wm(cmn12%locjs), mband, wm(ml3), meband)
      endif
!
! multiply jacobian by scalar, add identity, and do lu decomposition.
      con = -hrl1
      call dscal (lenp, con, wm(3), 1 )
      ii = mband + 2
      do 580 i = 1,cmn12%n
        wm(ii) = wm(ii) + one
 580    ii = ii + meband
      cmn12%nlu = cmn12%nlu + 1
      call dgbfa (wm(3), meband, cmn12%n, ml, mu, iwm(31), ier)
      if (ier .ne. 0) ierpj = 1
      return
! end of code block for miter = 4 or 5. --------------------------------
!
!----------------------- end of subroutine dvjac -----------------------
      end subroutine dvjac
!deck dacopy
      subroutine dacopy (nrow, ncol, a, nrowa, b, nrowb)
      double precision a, b
      integer nrow, ncol, nrowa, nrowb
      dimension a(nrowa,ncol), b(nrowb,ncol)
!-----------------------------------------------------------------------
! call sequence input -- nrow, ncol, a, nrowa, nrowb
! call sequence output -- b
! common block variables accessed -- none
!
! subroutines called by dacopy: dcopy
! function routines called by dacopy: none
!-----------------------------------------------------------------------
! this routine copies one rectangular array, a, to another, b,
! where a and b may have different row dimensions, nrowa and nrowb.
! the data copied consists of nrow rows and ncol columns.
!-----------------------------------------------------------------------
      integer ic
!
      do 20 ic = 1,ncol
        call dcopy (nrow, a(1,ic), 1, b(1,ic), 1)
 20     continue
!
      return
!----------------------- end of subroutine dacopy ----------------------
      end subroutine dacopy
!deck dvsol
      subroutine dvsol (wm, iwm, x, iersl, cmn12)
      implicit none
      double precision wm, x
      integer iwm, iersl
      dimension wm(*), iwm(*), x(*)
      type( dvode_cmn_vars ) :: cmn12
!-----------------------------------------------------------------------
! call sequence input -- wm, iwm, x
! call sequence output -- x, iersl
! common block variables accessed:
!     /dvod_cmn01/ -- h, rl1, miter, n
!
! subroutines called by dvsol: dgesl, dgbsl
! function routines called by dvsol: none
!-----------------------------------------------------------------------
! this routine manages the solution of the linear system arising from
! a chord iteration.  it is called if miter .ne. 0.
! if miter is 1 or 2, it calls dgesl to accomplish this.
! if miter = 3 it updates the coefficient h*rl1 in the diagonal
! matrix, and then computes the solution.
! if miter is 4 or 5, it calls dgbsl.
! communication with dvsol uses the following variables:
! wm    = real work space containing the inverse diagonal matrix if
!         miter = 3 and the lu decomposition of the matrix otherwise.
!         storage of matrix elements starts at wm(3).
!         wm also contains the following matrix-related data:
!         wm(1) = sqrt(uround) (not used here),
!         wm(2) = hrl1, the previous value of h*rl1, used if miter = 3.
! iwm   = integer work space containing pivot information, starting at
!         iwm(31), if miter is 1, 2, 4, or 5.  iwm also contains band
!         parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
! x     = the right-hand side vector on input, and the solution vector
!         on output, of length n.
! iersl = output flag.  iersl = 0 if no trouble occurred.
!         iersl = 1 if a singular matrix arose with miter = 3.
!-----------------------------------------------------------------------
!!
!! type declarations for labeled common block dvod_cmn01 --------------------
!!
!      double precision acnrm, ccmxj, conp, crate, drc, el,   &
!           eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!           rc, rl1, tau, tq, tn, uround
!      integer icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!              l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!              locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!              n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!              nslp, nyh
!
! type declarations for local variables --------------------------------
!
      integer i, meband, ml, mu
      double precision di, hrl1, one, phrl1, r, zero
!-----------------------------------------------------------------------
! the following fortran-77 declaration is to cause the values of the
! listed (local) variables to be saved between calls to this integrator.
!-----------------------------------------------------------------------
      save one, zero
!
!      common /dvod_cmn01/ acnrm, ccmxj, conp, crate, drc, el(13),   &
!                      eta, etamax, h, hmin, hmxi, hnew, hscal, prl1,   &
!                      rc, rl1, tau(13), tq(5), tn, uround,   &
!                      icf, init, ipup, jcur, jstart, jsv, kflag, kuth,   &
!                      l, lmax, lyh, lewt, lacor, lsavf, lwm, liwm,   &
!                      locjs, maxord, meth, miter, msbj, mxhnil, mxstep,   &
!                      n, newh, newq, nhnil, nq, nqnyh, nqwait, nslj,   &
!                      nslp, nyh
!
      data one /1.0d0/, zero /0.0d0/
!
      iersl = 0
      go to (100, 100, 300, 400, 400), cmn12%miter
 100  call dgesl (wm(3), cmn12%n, cmn12%n, iwm(31), x, 0)
      return
!
 300  phrl1 = wm(2)
      hrl1 = cmn12%h*cmn12%rl1
      wm(2) = hrl1
      if (hrl1 .eq. phrl1) go to 330
      r = hrl1/phrl1
      do 320 i = 1,cmn12%n
        di = one - r*(one - one/wm(i+2))
        if (abs(di) .eq. zero) go to 390
 320    wm(i+2) = one/di
!
 330  do 340 i = 1,cmn12%n
 340    x(i) = wm(i+2)*x(i)
      return
 390  iersl = 1
      return
!
 400  ml = iwm(1)
      mu = iwm(2)
      meband = 2*ml + mu + 1
      call dgbsl (wm(3), meband, cmn12%n, ml, mu, iwm(31), x, 0)
      return
!----------------------- end of subroutine dvsol -----------------------
      end subroutine dvsol
!!deck dvsrco
!      subroutine dvsrco (rsav, isav, job)
!      double precision rsav
!      integer isav, job
!      dimension rsav(*), isav(*)
!!-----------------------------------------------------------------------
!! call sequence input -- rsav, isav, job
!! call sequence output -- rsav, isav
!! common block variables accessed -- all of /dvod_cmn01/ and /dvod_cmn02/
!!
!! subroutines/functions called by dvsrco: none
!!-----------------------------------------------------------------------
!! this routine saves or restores (depending on job) the contents of the
!! common blocks dvod_cmn01 and dvod_cmn02, which are used internally by dvode.
!!
!! rsav = real array of length 49 or more.
!! isav = integer array of length 41 or more.
!! job  = flag indicating to save or restore the common blocks:
!!        job  = 1 if common is to be saved (written to rsav/isav).
!!        job  = 2 if common is to be restored (read from rsav/isav).
!!        a call with job = 2 presumes a prior call with job = 1.
!!-----------------------------------------------------------------------
!      double precision rvod1, rvod2
!      integer ivod1, ivod2
!      integer i, leniv1, leniv2, lenrv1, lenrv2
!!-----------------------------------------------------------------------
!! the following fortran-77 declaration is to cause the values of the
!! listed (local) variables to be saved between calls to this integrator.
!!-----------------------------------------------------------------------
!      save lenrv1, leniv1, lenrv2, leniv2
!!
!      common /dvod_cmn01/ rvod1(48), ivod1(33)
!      common /dvod_cmn02/ rvod2(1), ivod2(8)
!      data lenrv1/48/, leniv1/33/, lenrv2/1/, leniv2/8/
!!
!      if (job .eq. 2) go to 100
!      do 10 i = 1,lenrv1
! 10     rsav(i) = rvod1(i)
!      do 15 i = 1,lenrv2
! 15     rsav(lenrv1+i) = rvod2(i)
!!
!      do 20 i = 1,leniv1
! 20     isav(i) = ivod1(i)
!      do 25 i = 1,leniv2
! 25     isav(leniv1+i) = ivod2(i)
!!
!      return
!!
! 100  continue
!      do 110 i = 1,lenrv1
! 110     rvod1(i) = rsav(i)
!      do 115 i = 1,lenrv2
! 115     rvod2(i) = rsav(lenrv1+i)
!!
!      do 120 i = 1,leniv1
! 120     ivod1(i) = isav(i)
!      do 125 i = 1,leniv2
! 125     ivod2(i) = isav(leniv1+i)
!!
!      return
!!----------------------- end of subroutine dvsrco ----------------------
!      end subroutine dvsrco
!deck dewset
      subroutine dewset (n, itol, rtol, atol, ycur, ewt)
!***begin prologue  dewset
!***subsidiary
!***purpose  set error weight vector.
!***type      double precision (sewset-s, dewset-d)
!***author  hindmarsh, alan c., (llnl)
!***description
!
!  this subroutine sets the error weight vector ewt according to
!      ewt(i) = rtol(i)*abs(ycur(i)) + atol(i),  i = 1,...,n,
!  with the subscript on rtol and/or atol possibly replaced by 1 above,
!  depending on the value of itol.
!
!***see also  dlsode
!***routines called  (none)
!***revision history  (yymmdd)
!   791129  date written
!   890501  modified prologue to slatec/ldoc format.  (fnf)
!   890503  minor cosmetic changes.  (fnf)
!   930809  renamed to allow single/double precision versions. (ach)
!***end prologue  dewset
!**end
      integer n, itol
      integer i
      double precision rtol, atol, ycur, ewt
      dimension rtol(*), atol(*), ycur(n), ewt(n)
!
!***first executable statement  dewset
      go to (10, 20, 30, 40), itol
 10   continue
      do 15 i = 1,n
 15     ewt(i) = rtol(1)*abs(ycur(i)) + atol(1)
      return
 20   continue
      do 25 i = 1,n
 25     ewt(i) = rtol(1)*abs(ycur(i)) + atol(i)
      return
 30   continue
      do 35 i = 1,n
 35     ewt(i) = rtol(i)*abs(ycur(i)) + atol(1)
      return
 40   continue
      do 45 i = 1,n
 45     ewt(i) = rtol(i)*abs(ycur(i)) + atol(i)
      return
!----------------------- end of subroutine dewset ----------------------
      end subroutine dewset
!deck dvnorm
      double precision function dvnorm (n, v, w)
!***begin prologue  dvnorm
!***subsidiary
!***purpose  weighted root-mean-square vector norm.
!***type      double precision (svnorm-s, dvnorm-d)
!***author  hindmarsh, alan c., (llnl)
!***description
!
!  this function routine computes the weighted root-mean-square norm
!  of the vector of length n contained in the array v, with weights
!  contained in the array w of length n:
!    dvnorm = sqrt( (1/n) * sum( v(i)*w(i) )**2 )
!
!***see also  dlsode
!***routines called  (none)
!***revision history  (yymmdd)
!   791129  date written
!   890501  modified prologue to slatec/ldoc format.  (fnf)
!   890503  minor cosmetic changes.  (fnf)
!   930809  renamed to allow single/double precision versions. (ach)
!***end prologue  dvnorm
!**end
      integer n,   i
      double precision v, w,   sum
      dimension v(n), w(n)
!
!***first executable statement  dvnorm
      sum = 0.0d0
      do 10 i = 1,n
 10     sum = sum + (v(i)*w(i))**2
      dvnorm = sqrt(sum/n)
      return
!----------------------- end of function dvnorm ------------------------
      end function dvnorm
!deck xerrwd
      subroutine xerrwd (msg, nmes, nerr, level, ni, i1, i2, nr, r1, r2)
!***begin prologue  xerrwd
!***subsidiary
!***purpose  write error message with values.
!***category  r3c
!***type      double precision (xerrwv-s, xerrwd-d)
!***author  hindmarsh, alan c., (llnl)
!***description
!
!  subroutines xerrwd, xsetf, xsetun, and the function routine ixsav,
!  as given here, constitute a simplified version of the slatec error
!  handling package.
!
!  all arguments are input arguments.
!
!  msg    = the message (character array).
!  nmes   = the length of msg (number of characters).
!  nerr   = the error number (not used).
!  level  = the error level..
!           0 or 1 means recoverable (control returns to caller).
!           2 means fatal (run is aborted--see note below).
!  ni     = number of integers (0, 1, or 2) to be printed with message.
!  i1,i2  = integers to be printed, depending on ni.
!  nr     = number of reals (0, 1, or 2) to be printed with message.
!  r1,r2  = reals to be printed, depending on nr.
!
!  note..  this routine is machine-dependent and specialized for use
!  in limited context, in the following ways..
!  1. the argument msg is assumed to be of type character, and
!     the message is printed with a format of (1x,a).
!  2. the message is assumed to take only one line.
!     multi-line messages are generated by repeated calls.
!  3. if level = 2, control passes to the statement   stop
!     to abort the run.  this statement may be machine-dependent.
!  4. r1 and r2 are assumed to be in double precision and are printed
!     in d21.13 format.
!
!***routines called  ixsav
!***revision history  (yymmdd)
!   920831  date written
!   921118  replaced mflgsv/lunsav by ixsav. (ach)
!   930329  modified prologue to slatec format. (fnf)
!   930407  changed msg from character*1 array to variable. (fnf)
!   930922  minor cosmetic change. (fnf)
!***end prologue  xerrwd
!
!*internal notes:
!
! for a different default logical unit number, ixsav (or a subsidiary
! routine that it calls) will need to be modified.
! for a different run-abort command, change the statement following
! statement 100 at the end.
!-----------------------------------------------------------------------
! subroutines called by xerrwd.. none
! function routine called by xerrwd.. ixsav
!-----------------------------------------------------------------------
!**end
!
!  declare arguments.
!
      double precision r1, r2
      integer nmes, nerr, level, ni, i1, i2, nr
      character*(*) msg
!
!  declare local variables.
!
! 27-oct-2005 rce - do not declare functions that are in the module
!     integer lunit, ixsav, mesflg
      integer lunit,        mesflg
!
!  get logical unit number and message print flag.
!
!***first executable statement  xerrwd
      lunit = ixsav (1, 0, .false.)
      mesflg = ixsav (2, 0, .false.)
      if (mesflg .eq. 0) go to 100
!
!  write the message.
!
      write (lunit,10)  msg
 10   format(1x,a)
      if (ni .eq. 1) write (lunit, 20) i1
 20   format(6x,'in above message,  i1 =',i10)
      if (ni .eq. 2) write (lunit, 30) i1,i2
 30   format(6x,'in above message,  i1 =',i10,3x,'i2 =',i10)
      if (nr .eq. 1) write (lunit, 40) r1
 40   format(6x,'in above message,  r1 =',d21.13)
      if (nr .eq. 2) write (lunit, 50) r1,r2
 50   format(6x,'in above,  r1 =',d21.13,3x,'r2 =',d21.13)
!
!  abort the run if level = 2.
!
 100  if (level .ne. 2) return
         call wrf_error_fatal(" ")

!----------------------- end of subroutine xerrwd ----------------------
      end subroutine xerrwd
!deck xsetf
      subroutine xsetf (mflag)
!***begin prologue  xsetf
!***purpose  reset the error print control flag.
!***category  r3a
!***type      all (xsetf-a)
!***keywords  error control
!***author  hindmarsh, alan c., (llnl)
!***description
!
!   xsetf sets the error print control flag to mflag:
!      mflag=1 means print all messages (the default).
!      mflag=0 means no printing.
!
!***see also  xerrwd, xerrwv
!***references  (none)
!***routines called  ixsav
!***revision history  (yymmdd)
!   921118  date written
!   930329  added slatec format prologue. (fnf)
!   930407  corrected see also section. (fnf)
!   930922  made user-callable, and other cosmetic changes. (fnf)
!***end prologue  xsetf
!
! subroutines called by xsetf.. none
! function routine called by xsetf.. ixsav
!-----------------------------------------------------------------------
!**end
! 27-oct-2005 rce - do not declare functions that are in the module
!     integer mflag, junk, ixsav
      integer mflag, junk
!
!***first executable statement  xsetf
      if (mflag .eq. 0 .or. mflag .eq. 1) junk = ixsav (2,mflag,.true.)
      return
!----------------------- end of subroutine xsetf -----------------------
      end subroutine xsetf
!deck xsetun
      subroutine xsetun (lun)
!***begin prologue  xsetun
!***purpose  reset the logical unit number for error messages.
!***category  r3b
!***type      all (xsetun-a)
!***keywords  error control
!***description
!
!   xsetun sets the logical unit number for error messages to lun.
!
!***author  hindmarsh, alan c., (llnl)
!***see also  xerrwd, xerrwv
!***references  (none)
!***routines called  ixsav
!***revision history  (yymmdd)
!   921118  date written
!   930329  added slatec format prologue. (fnf)
!   930407  corrected see also section. (fnf)
!   930922  made user-callable, and other cosmetic changes. (fnf)
!***end prologue  xsetun
!
! subroutines called by xsetun.. none
! function routine called by xsetun.. ixsav
!-----------------------------------------------------------------------
!**end
! 27-oct-2005 rce - do not declare functions that are in the module
!     integer lun, junk, ixsav
      integer lun, junk
!
!***first executable statement  xsetun
      if (lun .gt. 0) junk = ixsav (1,lun,.true.)
      return
!----------------------- end of subroutine xsetun ----------------------
      end subroutine xsetun
!deck ixsav
      integer function ixsav (ipar, ivalue, iset)
!***begin prologue  ixsav
!***subsidiary
!***purpose  save and recall error message control parameters.
!***category  r3c
!***type      all (ixsav-a)
!***author  hindmarsh, alan c., (llnl)
!***description
!
!  ixsav saves and recalls one of two error message parameters:
!    lunit, the logical unit number to which messages are printed, and
!    mesflg, the message print flag.
!  this is a modification of the slatec library routine j4save.
!
!  saved local variables..
!   lunit  = logical unit number for messages.  the default is obtained
!            by a call to iumach (may be machine-dependent).
!   mesflg = print control flag..
!            1 means print all messages (the default).
!            0 means no printing.
!
!  on input..
!    ipar   = parameter indicator (1 for lunit, 2 for mesflg).
!    ivalue = the value to be set for the parameter, if iset = .true.
!    iset   = logical flag to indicate whether to read or write.
!             if iset = .true., the parameter will be given
!             the value ivalue.  if iset = .false., the parameter
!             will be unchanged, and ivalue is a dummy argument.
!
!  on return..
!    ixsav = the (old) value of the parameter.
!
!***see also  xerrwd, xerrwv
!***routines called  iumach
!***revision history  (yymmdd)
!   921118  date written
!   930329  modified prologue to slatec format. (fnf)
!   930915  added iumach call to get default output unit.  (ach)
!   930922  minor cosmetic changes. (fnf)
!   010425  type declaration for iumach added. (ach)
!***end prologue  ixsav
!
! subroutines called by ixsav.. none
! function routine called by ixsav.. iumach
!-----------------------------------------------------------------------
!**end
      logical iset
      integer ipar, ivalue
!-----------------------------------------------------------------------
! 27-oct-2005 rce - do not declare functions that are in the module
!     integer iumach, lunit, mesflg
      integer         lunit, mesflg
!-----------------------------------------------------------------------
! the following fortran-77 declaration is to cause the values of the
! listed (local) variables to be saved between calls to this routine.
!-----------------------------------------------------------------------
      save lunit, mesflg
      data lunit/-1/, mesflg/1/
!
!***first executable statement  ixsav
      if (ipar .eq. 1) then
        if (lunit .eq. -1) lunit = iumach()
        ixsav = lunit
        if (iset) lunit = ivalue
        endif
!
      if (ipar .eq. 2) then
        ixsav = mesflg
        if (iset) mesflg = ivalue
        endif
!
      return
!----------------------- end of function ixsav -------------------------
      end function ixsav
!deck iumach
      integer function iumach()
!***begin prologue  iumach
!***purpose  provide standard output unit number.
!***category  r1
!***type      integer (iumach-i)
!***keywords  machine constants
!***author  hindmarsh, alan c., (llnl)
!***description
! *usage:
!        integer  lout, iumach
!        lout = iumach()
!
! *function return values:
!     lout : the standard logical unit for fortran output.
!
!***references  (none)
!***routines called  (none)
!***revision history  (yymmdd)
!   930915  date written
!   930922  made user-callable, and other cosmetic changes. (fnf)
!***end prologue  iumach
!
!*internal notes:
!  the built-in value of 6 is standard on a wide range of fortran
!  systems.  this may be machine-dependent.
!**end
!***first executable statement  iumach
      iumach = 6
!
      return
!----------------------- end of function iumach ------------------------
      end function iumach
!deck dumach
      double precision function dumach ()
!***begin prologue  dumach
!***purpose  compute the unit roundoff of the machine.
!***category  r1
!***type      double precision (rumach-s, dumach-d)
!***keywords  machine constants
!***author  hindmarsh, alan c., (llnl)
!***description
! *usage:
!        double precision  a, dumach
!        a = dumach()
!
! *function return values:
!     a : the unit roundoff of the machine.
!
! *description:
!     the unit roundoff is defined as the smallest positive machine
!     number u such that  1.0 + u .ne. 1.0.  this is computed by dumach
!     in a machine-independent manner.
!
!***references  (none)
!***routines called  dumsum
!***revision history  (yyyymmdd)
!   19930216  date written
!   19930818  added slatec-format prologue.  (fnf)
!   20030707  added dumsum to force normal storage of comp.  (ach)
!***end prologue  dumach
!
      double precision u, comp
!***first executable statement  dumach
      u = 1.0d0
 10   u = u*0.5d0
      call dumsum(1.0d0, u, comp)
      if (comp .ne. 1.0d0) go to 10
      dumach = u*2.0d0
      return
!----------------------- end of function dumach ------------------------
      end function dumach
      subroutine dumsum(a,b,c)
!     routine to force normal storing of a + b, for dumach.
      double precision a, b, c
      c = a + b
      return
      end subroutine dumsum

!-----------------------------------------------------------------------
!   vode_subs.f - created on 28-jul-2004
!       by downloading following from www.netlib.org
!       1.  daxpy, dcopy, ddot, dnrm2, dscal, idamax from blas
!       2.  dgbfa, dbgsl, dgefa, dgesl from linpack
!
!   27-oct-2005 rce - change '1' dimensions in subr dgefa & dgesl
!-----------------------------------------------------------------------


!-----------------------------------------------------------------------
      subroutine daxpy(n,da,dx,incx,dy,incy)
!
!     constant times a vector plus a vector.
!     uses unrolled loops for increments equal to one.
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
      double precision dx(*),dy(*),da
      integer i,incx,incy,ix,iy,m,mp1,n
!
      if(n.le.0)return
      if (da .eq. 0.0d0) return
      if(incx.eq.1.and.incy.eq.1)go to 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dy(iy) + da*dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 m = mod(n,4)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dy(i) + da*dx(i)
   30 continue
      if( n .lt. 4 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,4
        dy(i) = dy(i) + da*dx(i)
        dy(i + 1) = dy(i + 1) + da*dx(i + 1)
        dy(i + 2) = dy(i + 2) + da*dx(i + 2)
        dy(i + 3) = dy(i + 3) + da*dx(i + 3)
   50 continue
      return
      end subroutine daxpy


!-----------------------------------------------------------------------
      subroutine  dcopy(n,dx,incx,dy,incy)
!
!     copies a vector, x, to a vector, y.
!     uses unrolled loops for increments equal to one.
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
      double precision dx(*),dy(*)
      integer i,incx,incy,ix,iy,m,mp1,n
!
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dy(iy) = dx(ix)
        ix = ix + incx
        iy = iy + incy
   10 continue
      return
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 m = mod(n,7)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dy(i) = dx(i)
   30 continue
      if( n .lt. 7 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,7
        dy(i) = dx(i)
        dy(i + 1) = dx(i + 1)
        dy(i + 2) = dx(i + 2)
        dy(i + 3) = dx(i + 3)
        dy(i + 4) = dx(i + 4)
        dy(i + 5) = dx(i + 5)
        dy(i + 6) = dx(i + 6)
   50 continue
      return
      end subroutine  dcopy


      double precision function ddot(n,dx,incx,dy,incy)
!
!     forms the dot product of two vectors.
!     uses unrolled loops for increments equal to one.
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
      double precision dx(*),dy(*),dtemp
      integer i,incx,incy,ix,iy,m,mp1,n
!
      ddot = 0.0d0
      dtemp = 0.0d0
      if(n.le.0)return
      if(incx.eq.1.and.incy.eq.1)go to 20
!
!        code for unequal increments or equal increments
!          not equal to 1
!
      ix = 1
      iy = 1
      if(incx.lt.0)ix = (-n+1)*incx + 1
      if(incy.lt.0)iy = (-n+1)*incy + 1
      do 10 i = 1,n
        dtemp = dtemp + dx(ix)*dy(iy)
        ix = ix + incx
        iy = iy + incy
   10 continue
      ddot = dtemp
      return
!
!        code for both increments equal to 1
!
!
!        clean-up loop
!
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dtemp = dtemp + dx(i)*dy(i)
   30 continue
      if( n .lt. 5 ) go to 60
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +   &
         dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
   50 continue
   60 ddot = dtemp
      return
      end function ddot


!-----------------------------------------------------------------------
      double precision function dnrm2 ( n, x, incx )
!     .. scalar arguments ..
      integer                           incx, n
!     .. array arguments ..
      double precision                  x( * )
!     ..
!
!  dnrm2 returns the euclidean norm of a vector via the function
!  name, so that
!
!     dnrm2 := sqrt( x'*x )
!
!
!
!  -- this version written on 25-october-1982.
!     modified on 14-october-1993 to inline the call to dlassq.
!     sven hammarling, nag ltd.
!
!
!     .. parameters ..
      double precision      one         , zero
      parameter           ( one = 1.0d+0, zero = 0.0d+0 )
!     .. local scalars ..
      integer               ix
      double precision      absxi, norm, scale, ssq
!     .. intrinsic functions ..
      intrinsic             abs, sqrt
!     ..
!     .. executable statements ..
      if( n.lt.1 .or. incx.lt.1 )then
         norm  = zero
      else if( n.eq.1 )then
         norm  = abs( x( 1 ) )
      else
         scale = zero
         ssq   = one
!        the following loop is equivalent to this call to the lapack
!        auxiliary routine:
!        call dlassq( n, x, incx, scale, ssq )
!
         do 10, ix = 1, 1 + ( n - 1 )*incx, incx
            if( x( ix ).ne.zero )then
               absxi = abs( x( ix ) )
               if( scale.lt.absxi )then
                  ssq   = one   + ssq*( scale/absxi )**2
                  scale = absxi
               else
                  ssq   = ssq   +     ( absxi/scale )**2
               end if
            end if
   10    continue
         norm  = scale * sqrt( ssq )
      end if
!
      dnrm2 = norm
      return
!
!     end of dnrm2.
!
      end function dnrm2


!-----------------------------------------------------------------------
      subroutine  dscal(n,da,dx,incx)
!
!     scales a vector by a constant.
!     uses unrolled loops for increment equal to one.
!     jack dongarra, linpack, 3/11/78.
!     modified 3/93 to return if incx .le. 0.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
      double precision da,dx(*)
      integer i,incx,m,mp1,n,nincx
!
      if( n.le.0 .or. incx.le.0 )return
      if(incx.eq.1)go to 20
!
!        code for increment not equal to 1
!
      nincx = n*incx
      do 10 i = 1,nincx,incx
        dx(i) = da*dx(i)
   10 continue
      return
!
!        code for increment equal to 1
!
!
!        clean-up loop
!
   20 m = mod(n,5)
      if( m .eq. 0 ) go to 40
      do 30 i = 1,m
        dx(i) = da*dx(i)
   30 continue
      if( n .lt. 5 ) return
   40 mp1 = m + 1
      do 50 i = mp1,n,5
        dx(i) = da*dx(i)
        dx(i + 1) = da*dx(i + 1)
        dx(i + 2) = da*dx(i + 2)
        dx(i + 3) = da*dx(i + 3)
        dx(i + 4) = da*dx(i + 4)
   50 continue
      return
      end subroutine  dscal


!-----------------------------------------------------------------------
      integer function idamax(n,dx,incx)
!
!     finds the index of element having max. absolute value.
!     jack dongarra, linpack, 3/11/78.
!     modified 3/93 to return if incx .le. 0.
!     modified 12/3/93, array(1) declarations changed to array(*)
!
      double precision dx(*),dmax
      integer i,incx,ix,n
!
      idamax = 0
      if( n.lt.1 .or. incx.le.0 ) return
      idamax = 1
      if(n.eq.1)return
      if(incx.eq.1)go to 20
!
!        code for increment not equal to 1
!
      ix = 1
      dmax = dabs(dx(1))
      ix = ix + incx
      do 10 i = 2,n
         if(dabs(dx(ix)).le.dmax) go to 5
         idamax = i
         dmax = dabs(dx(ix))
    5    ix = ix + incx
   10 continue
      return
!
!        code for increment equal to 1
!
   20 dmax = dabs(dx(1))
      do 30 i = 2,n
         if(dabs(dx(i)).le.dmax) go to 30
         idamax = i
         dmax = dabs(dx(i))
   30 continue
      return
      end function idamax


!-----------------------------------------------------------------------
      subroutine dgbfa(abd,lda,n,ml,mu,ipvt,info)
      integer lda,n,ml,mu,ipvt(1),info
      double precision abd(lda,1)
!
!     dgbfa factors a double precision band matrix by elimination.
!
!     dgbfa is usually called by dgbco, but it can be called
!     directly with a saving in time if  rcond  is not needed.
!
!     on entry
!
!        abd     double precision(lda, n)
!                contains the matrix in band storage.  the columns
!                of the matrix are stored in the columns of  abd  and
!                the diagonals of the matrix are stored in rows
!                ml+1 through 2*ml+mu+1 of  abd .
!                see the comments below for details.
!
!        lda     integer
!                the leading dimension of the array  abd .
!                lda must be .ge. 2*ml + mu + 1 .
!
!        n       integer
!                the order of the original matrix.
!
!        ml      integer
!                number of diagonals below the main diagonal.
!                0 .le. ml .lt. n .
!
!        mu      integer
!                number of diagonals above the main diagonal.
!                0 .le. mu .lt. n .
!                more efficient if  ml .le. mu .
!     on return
!
!        abd     an upper triangular matrix in band storage and
!                the multipliers which were used to obtain it.
!                the factorization can be written  a = l*u  where
!                l  is a product of permutation and unit lower
!                triangular matrices and  u  is upper triangular.
!
!        ipvt    integer(n)
!                an integer vector of pivot indices.
!
!        info    integer
!                = 0  normal value.
!                = k  if  u(k,k) .eq. 0.0 .  this is not an error
!                     condition for this subroutine, but it does
!                     indicate that dgbsl will divide by zero if
!                     called.  use  rcond  in dgbco for a reliable
!                     indication of singularity.
!
!     band storage
!
!           if  a  is a band matrix, the following program segment
!           will set up the input.
!
!                   ml = (band width below the diagonal)
!                   mu = (band width above the diagonal)
!                   m = ml + mu + 1
!                   do 20 j = 1, n
!                      i1 = max0(1, j-mu)
!                      i2 = min0(n, j+ml)
!                      do 10 i = i1, i2
!                         k = i - j + m
!                         abd(k,j) = a(i,j)
!                10    continue
!                20 continue
!
!           this uses rows  ml+1  through  2*ml+mu+1  of  abd .
!           in addition, the first  ml  rows in  abd  are used for
!           elements generated during the triangularization.
!           the total number of rows needed in  abd  is  2*ml+mu+1 .
!           the  ml+mu by ml+mu  upper left triangle and the
!           ml by ml  lower right triangle are not referenced.
!
!     linpack. this version dated 08/14/78 .
!     cleve moler, university of new mexico, argonne national lab.
!
!     subroutines and functions
!
!     blas daxpy,dscal,idamax
!     fortran max0,min0
!
!     internal variables
!
      double precision t
! 27-oct-2005 rce - do not declare functions that are in the module
!     integer i,idamax,i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
      integer i,       i0,j,ju,jz,j0,j1,k,kp1,l,lm,m,mm,nm1
!
!
      m = ml + mu + 1
      info = 0
!
!     zero initial fill-in columns
!
      j0 = mu + 2
      j1 = min0(n,m) - 1
      if (j1 .lt. j0) go to 30
      do 20 jz = j0, j1
         i0 = m + 1 - jz
         do 10 i = i0, ml
            abd(i,jz) = 0.0d0
   10    continue
   20 continue
   30 continue
      jz = j1
      ju = 0
!
!     gaussian elimination with partial pivoting
!
      nm1 = n - 1
      if (nm1 .lt. 1) go to 130
      do 120 k = 1, nm1
         kp1 = k + 1
!
!        zero next fill-in column
!
         jz = jz + 1
         if (jz .gt. n) go to 50
         if (ml .lt. 1) go to 50
            do 40 i = 1, ml
               abd(i,jz) = 0.0d0
   40       continue
   50    continue
!
!        find l = pivot index
!
         lm = min0(ml,n-k)
         l = idamax(lm+1,abd(m,k),1) + m - 1
         ipvt(k) = l + k - m
!
!        zero pivot implies this column already triangularized
!
         if (abd(l,k) .eq. 0.0d0) go to 100
!
!           interchange if necessary
!
            if (l .eq. m) go to 60
               t = abd(l,k)
               abd(l,k) = abd(m,k)
               abd(m,k) = t
   60       continue
!
!           compute multipliers
!
            t = -1.0d0/abd(m,k)
            call dscal(lm,t,abd(m+1,k),1)
!
!           row elimination with column indexing
!
            ju = min0(max0(ju,mu+ipvt(k)),n)
            mm = m
            if (ju .lt. kp1) go to 90
            do 80 j = kp1, ju
               l = l - 1
               mm = mm - 1
               t = abd(l,j)
               if (l .eq. mm) go to 70
                  abd(l,j) = abd(mm,j)
                  abd(mm,j) = t
   70          continue
               call daxpy(lm,t,abd(m+1,k),1,abd(mm+1,j),1)
   80       continue
   90       continue
         go to 110
  100    continue
            info = k
  110    continue
  120 continue
  130 continue
      ipvt(n) = n
      if (abd(m,n) .eq. 0.0d0) info = n
      return
      end subroutine dgbfa


!-----------------------------------------------------------------------
      subroutine dgbsl(abd,lda,n,ml,mu,ipvt,b,job)
      integer lda,n,ml,mu,ipvt(1),job
      double precision abd(lda,1),b(1)
!
!     dgbsl solves the double precision band system
!     a * x = b  or  trans(a) * x = b
!     using the factors computed by dgbco or dgbfa.
!
!     on entry
!
!        abd     double precision(lda, n)
!                the output from dgbco or dgbfa.
!
!        lda     integer
!                the leading dimension of the array  abd .
!
!        n       integer
!                the order of the original matrix.
!
!        ml      integer
!                number of diagonals below the main diagonal.
!
!        mu      integer
!                number of diagonals above the main diagonal.
!
!        ipvt    integer(n)
!                the pivot vector from dgbco or dgbfa.
!
!        b       double precision(n)
!                the right hand side vector.
!
!        job     integer
!                = 0         to solve  a*x = b ,
!                = nonzero   to solve  trans(a)*x = b , where
!                            trans(a)  is the transpose.
!
!     on return
!
!        b       the solution vector  x .
!
!     error condition
!
!        a division by zero will occur if the input factor contains a
!        zero on the diagonal.  technically this indicates singularity
!        but it is often caused by improper arguments or improper
!        setting of lda .  it will not occur if the subroutines are
!        called correctly and if dgbco has set rcond .gt. 0.0
!        or dgbfa has set info .eq. 0 .
!
!     to compute  inverse(a) * c  where  c  is a matrix
!     with  p  columns
!           call dgbco(abd,lda,n,ml,mu,ipvt,rcond,z)
!           if (rcond is too small) go to ...
!           do 10 j = 1, p
!              call dgbsl(abd,lda,n,ml,mu,ipvt,c(1,j),0)
!        10 continue
!
!     linpack. this version dated 08/14/78 .
!     cleve moler, university of new mexico, argonne national lab.
!
!     subroutines and functions
!
!     blas daxpy,ddot
!     fortran min0
!
!     internal variables
!
! 27-oct-2005 rce - do not declare functions that are in the module
!     double precision ddot,t
      double precision      t
      integer k,kb,l,la,lb,lm,m,nm1
!
      m = mu + ml + 1
      nm1 = n - 1
      if (job .ne. 0) go to 50
!
!        job = 0 , solve  a * x = b
!        first solve l*y = b
!
         if (ml .eq. 0) go to 30
         if (nm1 .lt. 1) go to 30
            do 20 k = 1, nm1
               lm = min0(ml,n-k)
               l = ipvt(k)
               t = b(l)
               if (l .eq. k) go to 10
                  b(l) = b(k)
                  b(k) = t
   10          continue
               call daxpy(lm,t,abd(m+1,k),1,b(k+1),1)
   20       continue
   30    continue
!
!        now solve  u*x = y
!
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/abd(m,k)
            lm = min0(k,m) - 1
            la = m - lm
            lb = k - lm
            t = -b(k)
            call daxpy(lm,t,abd(la,k),1,b(lb),1)
   40    continue
      go to 100
   50 continue
!
!        job = nonzero, solve  trans(a) * x = b
!        first solve  trans(u)*y = b
!
         do 60 k = 1, n
            lm = min0(k,m) - 1
            la = m - lm
            lb = k - lm
            t = ddot(lm,abd(la,k),1,b(lb),1)
            b(k) = (b(k) - t)/abd(m,k)
   60    continue
!
!        now solve trans(l)*x = y
!
         if (ml .eq. 0) go to 90
         if (nm1 .lt. 1) go to 90
            do 80 kb = 1, nm1
               k = n - kb
               lm = min0(ml,n-k)
               b(k) = b(k) + ddot(lm,abd(m+1,k),1,b(k+1),1)
               l = ipvt(k)
               if (l .eq. k) go to 70
                  t = b(l)
                  b(l) = b(k)
                  b(k) = t
   70          continue
   80       continue
   90    continue
  100 continue
      return
      end subroutine dgbsl


!-----------------------------------------------------------------------
      subroutine dgefa(a,lda,n,ipvt,info)
! 27-oct-2005 rce - change '1' dimensions
!     integer lda,n,ipvt(1),info
      integer lda,n,ipvt(n),info
!     double precision a(lda,1)
      double precision a(lda,n)
!
!     dgefa factors a double precision matrix by gaussian elimination.
!
!     dgefa is usually called by dgeco, but it can be called
!     directly with a saving in time if  rcond  is not needed.
!     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
!
!     on entry
!
!        a       double precision(lda, n)
!                the matrix to be factored.
!
!        lda     integer
!                the leading dimension of the array  a .
!
!        n       integer
!                the order of the matrix  a .
!
!     on return
!
!        a       an upper triangular matrix and the multipliers
!                which were used to obtain it.
!                the factorization can be written  a = l*u  where
!                l  is a product of permutation and unit lower
!                triangular matrices and  u  is upper triangular.
!
!        ipvt    integer(n)
!                an integer vector of pivot indices.
!
!        info    integer
!                = 0  normal value.
!                = k  if  u(k,k) .eq. 0.0 .  this is not an error
!                     condition for this subroutine, but it does
!                     indicate that dgesl or dgedi will divide by zero
!                     if called.  use  rcond  in dgeco for a reliable
!                     indication of singularity.
!
!     linpack. this version dated 08/14/78 .
!     cleve moler, university of new mexico, argonne national lab.
!
!     subroutines and functions
!
!     blas daxpy,dscal,idamax
!
!     internal variables
!
      double precision t
! 27-oct-2005 rce - do not declare functions that are in the module
!     integer idamax,j,k,kp1,l,nm1
      integer        j,k,kp1,l,nm1
!
!
!     gaussian elimination with partial pivoting
!
      info = 0
      nm1 = n - 1
      if (nm1 .lt. 1) go to 70
      do 60 k = 1, nm1
         kp1 = k + 1
!
!        find l = pivot index
!
         l = idamax(n-k+1,a(k,k),1) + k - 1
         ipvt(k) = l
!
!        zero pivot implies this column already triangularized
!
         if (a(l,k) .eq. 0.0d0) go to 40
!
!           interchange if necessary
!
            if (l .eq. k) go to 10
               t = a(l,k)
               a(l,k) = a(k,k)
               a(k,k) = t
   10       continue
!
!           compute multipliers
!
            t = -1.0d0/a(k,k)
            call dscal(n-k,t,a(k+1,k),1)
!
!           row elimination with column indexing
!
            do 30 j = kp1, n
               t = a(l,j)
               if (l .eq. k) go to 20
                  a(l,j) = a(k,j)
                  a(k,j) = t
   20          continue
               call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
   30       continue
         go to 50
   40    continue
            info = k
   50    continue
   60 continue
   70 continue
      ipvt(n) = n
      if (a(n,n) .eq. 0.0d0) info = n
      return
      end subroutine dgefa


!-----------------------------------------------------------------------
      subroutine dgesl(a,lda,n,ipvt,b,job)
! 27-oct-2005 rce - change '1' dimensions
!     integer lda,n,ipvt(1),job
      integer lda,n,ipvt(n),job
!     double precision a(lda,1),b(1)
      double precision a(lda,n),b(n)
!
!     dgesl solves the double precision system
!     a * x = b  or  trans(a) * x = b
!     using the factors computed by dgeco or dgefa.
!
!     on entry
!
!        a       double precision(lda, n)
!                the output from dgeco or dgefa.
!
!        lda     integer
!                the leading dimension of the array  a .
!
!        n       integer
!                the order of the matrix  a .
!
!        ipvt    integer(n)
!                the pivot vector from dgeco or dgefa.
!
!        b       double precision(n)
!                the right hand side vector.
!
!        job     integer
!                = 0         to solve  a*x = b ,
!                = nonzero   to solve  trans(a)*x = b  where
!                            trans(a)  is the transpose.
!
!     on return
!
!        b       the solution vector  x .
!
!     error condition
!
!        a division by zero will occur if the input factor contains a
!        zero on the diagonal.  technically this indicates singularity
!        but it is often caused by improper arguments or improper
!        setting of lda .  it will not occur if the subroutines are
!        called correctly and if dgeco has set rcond .gt. 0.0
!        or dgefa has set info .eq. 0 .
!
!     to compute  inverse(a) * c  where  c  is a matrix
!     with  p  columns
!           call dgeco(a,lda,n,ipvt,rcond,z)
!           if (rcond is too small) go to ...
!           do 10 j = 1, p
!              call dgesl(a,lda,n,ipvt,c(1,j),0)
!        10 continue
!
!     linpack. this version dated 08/14/78 .
!     cleve moler, university of new mexico, argonne national lab.
!
!     subroutines and functions
!
!     blas daxpy,ddot
!
!     internal variables
!
! 27-oct-2005 rce - do not declare functions that are in the module
!     double precision ddot,t
      double precision      t
      integer k,kb,l,nm1
!
      nm1 = n - 1
      if (job .ne. 0) go to 50
!
!        job = 0 , solve  a * x = b
!        first solve  l*y = b
!
         if (nm1 .lt. 1) go to 30
         do 20 k = 1, nm1
            l = ipvt(k)
            t = b(l)
            if (l .eq. k) go to 10
               b(l) = b(k)
               b(k) = t
   10       continue
            call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
   20    continue
   30    continue
!
!        now solve  u*x = y
!
         do 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            call daxpy(k-1,t,a(1,k),1,b(1),1)
   40    continue
      go to 100
   50 continue
!
!        job = nonzero, solve  trans(a) * x = b
!        first solve  trans(u)*y = b
!
         do 60 k = 1, n
            t = ddot(k-1,a(1,k),1,b(1),1)
            b(k) = (b(k) - t)/a(k,k)
   60    continue
!
!        now solve trans(l)*x = y
!
         if (nm1 .lt. 1) go to 90
         do 80 kb = 1, nm1
            k = n - kb
            b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
            l = ipvt(k)
            if (l .eq. k) go to 70
               t = b(l)
               b(l) = b(k)
               b(k) = t
   70       continue
   80    continue
   90    continue
  100 continue
      return
      end subroutine dgesl



      end module module_cmu_dvode_solver