      subroutine mirkdc(method,tol,neqns,Nsub,MxNsub,
     +            mesh_input,mesh,sol_input,Y,ldY,output_control,info,
     +            iwork,work,init_de,init_Y,fsub,gsub,dfsub,dgsub)
c
c MIRKDC_4, Dec. 1999.
c
c Authors: Wayne Enright (enright@cs.toronto.edu),
c          Paul Muir (muir@stmarys.ca).
c Additional programming assistance: Mark Adams, Nicole DeYoung.
c Thanks to Beta-Testers: Ian Gladwell, Larry Shampine, Richard Pancer.
c
c Note: This code makes use of the Level 1 BLAS routines. These must
c be linked with the MIRKDC package in order for it to work.
c
c***declaration of constants***
      integer Mxs
      parameter(Mxs=10)
c
c     `Mxs' Maximum number of stages of the Runge-Kutta method.
c
c***declaration of parameters***
c   imports:
      integer               method
      double precision      tol
      integer               neqns, Nsub, MxNsub
      integer               mesh_input
      double precision      mesh(0:MxNsub)
      integer               sol_input
      double precision      Y(ldY*(Nsub+1))
      integer               ldY, output_control
c
c     `method' Defines the MIRK method to be used. Possible values
c              for `method' and the corresponding MIRK schemes are:
c                    method  | MIRK formula
c                    ---------------------------------------------------------
c                    Externally documented methods.
c                     2      | MIRK221 scheme (trapezoidal rule, second order)
c                     4      | MIRK343 scheme (Lobatto scheme, fourth order)
c                     6      | MIRK563 scheme (6th order)
c                    Internally available methods.
c                    121     | MIRK121 scheme (midpoint rule, second order)
c                    221     | MIRK221 scheme (trapezoidal rule, second order)
c                    232     | MIRK232 scheme (one-sided, abscissa 0, 2/3)
c                    2321    | MIRK232 scheme (one-sided, abscissa 1, 1/3)
c                    343     | MIRK343 scheme (Lobatto scheme, fourth order)
c                    453     | MIRK453 scheme (one-sided, 5th order)
c                    563     | MIRK563 scheme (6th order)
c
c     `tol' Tolerance applied to an estimate of the defect of the approximate
c           solution; the defect is the amount by which the continuous
c           approximate solution, u(t), fails to satisfy the ODE system.
c           We require |def(t)|/(|f(t,u(t))|+1) to be less than `tol', where
c           y'(t)=f(t,y(t)) is the ODE and def(t)= u'(t)-f(t,u(t)).
c           The same tolerance is applied to all defect components.
c
c     `neqns' Number of differential equations. Also number of boundary
c                                                                  conditions.
c     `Nsub' On input, number of subintervals into which the problem
c            interval is initially partitioned.
c     `MxNsub' Maximum number of subintervals.
c     `mesh_input' a) A value of 0 indicates that a uniform initial mesh
c                     should be set up by `mirkdc'.
c                  b) A value of 1 indicates that an initial mesh has
c                     already been provided by the user in the array `mesh'.
c                     This option should be chosen if the user has specific
c                     knowledge of the solution behavior.
c                  c) A value of 2 indicates that the user is employing
c                     continuation to solve the problem and that the mesh
c                     from a previous run will be stored in the array `mesh'.
c     `mesh' On input, the set of points that partition the problem
c            interval; initially empty if `mesh_input' = 0.
c     `sol_input' a) Normal use of `mirkdc' requires that `sol_input'=1,
c                    indicating that an initial solution estimate is to be
c                    provided by the user through the `init_Y' routine. It
c                    is recommended that `Y' not be set to a constant.
c                 b) The option `sol_input' = 2 should only be chosen when
c                    using continuation to solve a problem, where `Y' takes
c                    its value from the solution of a previous run. A value
c                    of 2 indicates that an initial solution estimate at each
c                    mesh point is available in the array `Y'.
c     `Y' On input, if `sol_input'=2, this is the initial discrete
c         approximation to the solution, evaluated at each meshpoint.
c         Externally, in the interfaces with the user main program and
c         subroutines, `Y' is a two-dimensional array; the i-th column of `Y'
c         stores the solution approximation for the i-th mesh point. Each
c         solution approximation is of length `neqns'; the resulting
c         dimensions of `Y' are `neqns x (Nsub+1)'. Inside MIRKDC it is
c         convenient for `Y' to be stored as a column vector; in this
c         one-dimensional version the solutions at each meshpoint come one
c         after the other rather than side by side as they are in the
c         two-dimensional array.
c     `ldY' The leading dimension of Y. (`ldY' >= `neqns').
c     `output_control' Controls the amount of information output. Provides
c                      profiling of code execution for activities such as the
c                      Newton iteration, mesh selection, or defect estimation.
c                             0 - None. 1 - Intermediate. 2 - Full.
c   exports:
c     double precision      mesh(0:Nsub)
c     double precision      Y(ldY*(Nsub+1))
      integer               info
c
c     `mesh' On output, the set of points which define the final mesh.
c     `Y' On output, if `info'=0, the converged discrete solution, evaluated
c         on the final mesh.
c     `info' Communication flag:
c                0 - implies a successful termination - solution obtained such
c                    that the defect is within user specified tolerance.
c               -1 - implies an unsuccessful termination - the required size for
c                    the new mesh is be too large.
c            Internal flag:
c              The Newton iteration can fail for one of several reasons, each
c              signaled by a non-zero `info' value. Even if the Newton
c              iteration converges, the resultant solution may be
c              unsuitable because the defect is too large. During the
c              execution of `mirkdc', `info' may take on the values given
c              below. None of these conditions lead directly to an
c              unsuccessful termination. Rather, the code will react by
c              attempting to double the size of the mesh and try again.
c                1 - too many iterations were taken without convergence being
c                    obtained.
c                2 - a singular coefficient matrix was encountered during the
c                    attempted solution of the Newton system.
c                3 - (a) it was impossible to obtain a suitable damping factor
c                    for the Newton update (indicative of an "effectively
c                    singular" Jacobian) or (b) evaluation of natural criterion
c                    function overflowed, indicative of a divergent iteration.
c                4 - the defect was too large and the solution from the Newton
c                    system could not be trusted.
c   work space:
      integer               iwork(neqns*(MxNsub+1))
      double precision      work(2*MxNsub*Mxs*neqns +
     +                            5*neqns*MxNsub + 2*neqns**2*MxNsub +
     +                            6*neqns + 4*neqns**2 + 2*neqns**2*Mxs)
c
c     `iwork' Integer work array provided to the `newiter' routine.
c     `work' Double precision work array used within `mirkdc' and provided to
c            `NewIter', `defect_estimate', `mesh_selector', and `interp_eval'.
c            The largest amount of workspace used is during the call to
c            `Newiter'. The space required then is at most
c                    2 * MxNsub*Mxs*neqns + 5*neqns*MxNsub +
c                       2*neqns^2*MxNsub + 6*neqns + 4*neqns^2 + 2*neqns^2*Mxs.
c
c   user-supplied subroutines:
      external    init_de,init_Y,fsub,gsub,dfsub,dgsub
c
c     `init_de' User-supplied subroutine to initialize the problem dependent
c               variables, `leftbc', `a', and `b'.
c
c                  call init_de(leftbc, a, b)
c                    exports:
c                      integer leftbc
c                      double precision a,b
c                      `leftbc' Number of boundary conditions applied at the
c                               left endpoint of the problem interval.
c                      `a' Left endpoint of the problem interval.
c                      `b' Right endpoint of the problem interval.
c
c     `init_Y' User-supplied subroutine to initialize the discrete solution
c              approximation, `Y'. (Needed if `sol_input'=1).
c
c                  call init_Y(Nsub, neqns, mesh, Y, ldY)
c                    imports:
c                      integer Nsub, neqns
c                      double precision mesh(0:Nsub)
c                      integer ldY
c                      `Nsub' Number of subintervals in the initial mesh.
c                      `neqns' Number of differential equations in the
c                              first order system.
c                      `mesh' Points making up the partition of the problem
c                             interval.
c                      `ldY' Leading dimension of Y. (`ldY' >= `neqns').
c                    exports:
c                      double precision Y(ldY,Nsub+1)
c                      `Y' Initial approximation to the discrete
c                          solution at the points in the initial mesh.
c                          The i-th column of `Y' is the solution
c                          approximation vector at mesh point, `mesh(i-1)'.
c
c     `fsub' User-supplied subroutine which defines f(t,y) for the first order
c            system of differential equations, y' = f(t,y).
c
c                  call fsub(neqns, t, y, f)
c                    imports:
c                      integer neqns
c                      double precision t, y(neqns)
c                      `neqns' Number of differential equations in the
c                              first order system.
c                      `t' A point in the problem interval.
c                      `y' Current solution approximation at t.
c                    exports:
c                      double precision f(neqns)
c                      `f' Value of f(t,y).
c
c     `gsub' User-supplied subroutine which defines the boundary condition
c            equations.
c
c                  call gsub(neqns, ya, yb, bc)
c                    imports:
c                      integer neqns
c                      double precision ya(neqns), yb(neqns)
c                      `neqns' Number of differential equations in the
c                              first order system. Also the total number
c                              of boundary conditions.
c                      `ya' Current solution approximation at left endpoint.
c                      `yb' Current solution approximation at right endpoint.
c                    exports:
c                      double precision bc(neqns)
c                      `bc' Boundary condition equations evaluated at
c                           `ya' and `yb'. The first `leftbc' components
c                           of `bc' are `g_a(ya)'; the remaining
c                           `neqns-leftbc' components of `bc' are `g_b(yb)'.
c
c     `dfsub' User-supplied subroutine which defines the Jacobian, df/dy,
c             of the system of differential equations.  Since the array
c             `JJ' has already been set to zero, it is only necessary to
c             set the non-zero elements.
c
c                  call dfsub(neqns, t, y, JJ)
c                    imports:
c                      integer neqns
c                      double precision t, y(neqns)
c                      `neqns' Number of differential equations in the
c                              first order system.
c                      `t' A point in the problem interval.
c                      `y' Current solution approximation at t.
c                    exports:
c                      double precision JJ(neqns,neqns)
c                      `JJ' Contains the Jacobian, df/dy.
c
c     `dgsub' User-supplied subroutine which defines the Jacobian of
c             the boundary conditions, that is d(bc)/dy. The array `dbc'
c             has been set to zero prior to the call to `dgsub', so it
c             is only necessary for the user to set the non-zero elements.
c
c                  call dgsub(neqns, ya, yb, dbc)
c                    imports:
c                      integer neqns
c                      double precision ya(neqns), yb(neqns)
c                      `neqns' Number of differential equations in the
c                              first order system. Also the total number
c                              of boundary conditions.
c                      `ya' Current solution approximation at left endpoint.
c                      `yb' Current solution approximation at right endpoint.
c                    exports:
c                      double precision dbc(neqns,neqns)
c                      `dbc' Contains the Jacobian of the boundary conditions.
c
c***declaration of local variables***
      integer              leftbc
      double precision     a,b,h
      integer              maxiter
      double precision     newton_tol
      double precision     defect_norm
      integer              Nsub_star
      integer              i,j
c
c     `leftbc' The number of boundary conditions at the left end of the
c              problem interval.
c     `a' Left endpoint.
c     `b' Right endpoint.
c     `h' Initial subinterval size, when a uniform mesh is employed.
c     `maxiter' Maximum number of Newton iterations.
c     `newton_tol' Tolerance for each Newton iteration.
c     `defect_norm' Max norm of the defect estimate.
c     `Nsub_star' Number of subintervals of the new mesh.
c     `i,j' Loop indices.
c
c***declaration of index related info for the work array***
c
c   Arrays local to `mirkdc' to be represented within the `work' array:
c          `k_discrete', `k_interp', `defect', `mesh_new', `Y_new'.
c   Sizes:
      integer s_k_discrete, s_k_interp, s_defect
      integer s_mesh_new, s_Y_new
c
c     The actual array declarations will be replaced by assignments giving
c     the amount of space need within `work'.
c     double precision     k_discrete(Mxs*neqns*Nsub) =>
c                                       `s_k_discrete' = `Mxs*neqns*Nsub'
c     double precision     k_interp(Mxs*neqns*Nsub) =>
c                                        `s_k_interp' = `Mxs*neqns*Nsub'
c     double precision     defect(Nsub*neqns) => `s_defect' = `Nsub*neqns'
c     double precision     mesh_new(0:MxNsub) => `s_mesh_new' = `MxNsub+1'
c     double precision     Y_new((Nsub_star+1)*neqns) =>
c                                             `s_Y_new' = `(Nsub_star+1)*neqns'
c
c     `k_discrete' Stage information associated with the discrete
c                  Runge-Kutta method: the i-th set of `s*neqns' locations
c                  contain the `s' vectors of length `neqns' associated
c                  with the Runge-Kutta method on the i-th subinterval.
c                  (`Mxs' is the maximum number of stages; `s' is the
c                  actual number of stages.)
c     `k_interp' Stage information associated with the continuous
c                Runge-Kutta method; format is the same as `k_discrete'.
c     `defect'   Estimate of the (scaled) defect
c                                (u'(t) - f(t,u(t)))/(|f(t,u(t))|+1)
c                on each subinterval (where u(t) is the continuous
c                approximate solution). The defect estimate is a vector with
c                `neqns' components and there is one such vector for each
c                subinterval.
c     `mesh_new' Points defining the new mesh.
c     `Y_new'    Discrete approximation to the solution on `mesh_new'
c
c   Locations within `work':
      integer i_k_discrete, i_k_interp
      integer i_defect, i_mesh_new, i_Y_new
c
c     `i_k_discrete' = 1;
c     `i_k_interp' = `s_k_discrete' + 1
c     `i_defect' = `s_k_discrete' + `s_k_interp' + 1
c     `i_mesh_new' = `s_k_discrete' + `s_k_interp' + `s_defect' + 1
c     `i_Y_new' = `s_k_discrete'+`s_k_interp'+`s_defect'+`s_mesh_new'+1
c
c   Index related info for the representation within `work' of arrays local
c   to `newiter'. See `newiter' for explanations of the use of the arrays.
c   Sizes:
      integer s_PHI, s_delta, s_top, s_bot, s_blocks
c
c     `s_PHI' = `neqns * (Nsub+1)'
c     `s_delta' = `neqns * (Nsub+1)'
c     `s_top' = `neqns**2'
c     `s_bot' = `neqns**2'
c     `s_blocks' = `2 * neqns**2 * Nsub'
c
c   Locations within `work' during call to `newiter'.
      integer i_PHI, i_delta, i_top, i_bot, i_blocks
c
c     `i_PHI' = `s_k_discrete + s_k_interp + 1'
c     `i_delta' = `s_k_discrete + s_k_interp + s_PHI + 1'
c     `i_top' = `s_k_discrete + s_k_interp + s_PHI + s_delta + 1'
c     `i_bot' = `s_k_discrete + s_k_interp + s_PHI +  s_delta + s_top + 1'
c     `i_blocks' = `s_k_discrete+s_k_interp+s_PHI+s_delta+s_top+s_bot+1'
c
c   Index related info for the work array for arrays local to `defect_estimate'.c   See `defect_estimate' for explanations of the use of these arrays.
c   Sizes:
c     (All these arrays are of length `neqns'.)
c
c   Locations within `work' during the call to `defect_estimate'.
      integer i_z, i_z_prime, i_def_1, i_def_2
      integer i_f_1, i_f_2, i_temp_1, i_temp_2
c
c     `i_z = i_shift + 1; i_z_prime = i_shift + neqns + 1'
c     `i_def_1 = i_shift + 2*neqns + 1; i_def_2 = i_shift + 3*neqns + 1'
c     `i_f_1 = i_shift + 4*neqns + 1; i_f_2 = i_shift + 5*neqns + 1'
c     `i_temp_1 = i_shift + 6*neqns + 1; i_temp_2 = i_shift + 7*neqns + 1'
c     where `i_shift = s_k_discrete + s_k_interp + s_defect'
c
c   Other work indexes.
      integer i_work, i_shift
c
c     `i_work' Index to free part of work array passed to various routines
c              for use as a work array within those routines.
c     `i_shift' Temporary variable used as a shift to simplify `work' index
c               calculations.
c
c***declaration of variables for the /IOCNTRL/ common block***
c   imports:
      integer    print_level_0, print_level_1, print_level_2
      parameter  (print_level_0 = 0, print_level_1 = 1,
     *                print_level_2 = 2)
      integer    profile
      common /IOCNTRL/ profile
c
c     `print_level_0' No output.
c     `print_level_1' Intermediate output.
c     `print_level_2' Full output.
c     `profile' Controls output, to standard output, of profiling infor-
c               mation such as Newton iteration counts, mesh selection,
c               relative defect estimate. This variable is identical to
c               the `output_control' parameter.
c
c
c     The following two interfaces are included only to allow convenient loading
c     of the integer work array upon successful exit. This work array is needed
c     when the user calls the `sol_eval' routine, subsequent to the return from
c     `mirkdc'.
c
c***declaration of variables from /RK_s/ common block***
c   imports:
      integer s
      common /RK_s/ s
c     `s' Number of discrete stages.
c
c***declaration of variables from /interp_s_star/ common block***
c   imports:
      integer s_star
      common /interp_s_star/ s_star
c     `s_star' Total number of stages required to form the
c              interpolant on each subinterval. It includes all the
c              stages of the discrete formula plus the additional
c              stages required for the interpolant.
c----------------------------------------------------------------------------
c     Called by : main
c     Calls to  : `rk_tableau',`interp_tableau',`init_de',`init_Y',`NewIter',
c                 `defect_estimate',`mesh_selector',`interp_eval',`half_mesh'
c----------------------------------------------------------------------------
c***initialization***
c
c     Transfer output_control to common block variable `profile'.
      profile = output_control
c
      if (profile .GT. print_level_0) then
c       Output required size of `work' and `iwork' based on given
c       values for `neqns' and `MxNsub'.
        write(6,*)
        write(6,*) 'The input value for the order of the Runge-Kutta'
        write(6,80) 'method is ',method,'.'
        write(6,90) 'The input value for the number of ODEs is '
     +                                              ,neqns,'.'
        write(6,*) 'The input value for the maximum number of '
        write(6,100) 'subintervals is ', MxNsub,'.'
        write(6,*)
        write(6,*) 'Based on these values the size of the double'
        write(6,110) 'precision work array should be at least',
     +   MxNsub*Mxs*neqns +  5*neqns*MxNsub + 2*neqns**2*MxNsub +
     +                    6*neqns + 4*neqns**2 + 2*neqns**2*Mxs,'.'
        write(6,*) 'and the size of the integer work array should'
        write(6,120) 'be at least', neqns*(MxNsub+1),'.'
        write(6,*)
  80    format(1X,a10,i2,a1)
  90    format(1x,a41,i3,a1)
 100    format(1x,a15,i5,a1)
 110    format(1x,a39,i7,a1)
 120    format(1x,a11,i6,a1)
      endif
c
c     Set up formula dependent coefficients.
      call rk_tableau(method)
      call interp_tableau(method)
c
c     Initialize remaining ODE dependent variables.
      call init_de(leftbc,a,b)
c
c     Setup initial mesh.
c     If `mesh_input'=1 or 2 then there is already a mesh stored in `mesh'.
c     If `mesh_input'=0 then setup uniform initial mesh of `Nsub' subintervals.
      if (mesh_input.eq.0) then
            mesh(0) = a
            mesh(Nsub) = b
            h = (b - a)/Nsub
            do 10 i = 1, Nsub-1
               mesh(i) = i*h + a
 10         continue
      endif
c
      if (profile .GT. print_level_0) then
          write(6,130) 'The size of the initial mesh is '
     +                 ,Nsub,' subintervals.'
          write(6,*)
 130      format(1x,a31,i5,a14)
      end if
      if (profile .GT. print_level_1) then
          write(6,140)'The initial mesh of',Nsub,' subintervals:'
          write(6,150)(mesh(i),i=0,Nsub)
          write(6,*)
  140     format(1x,a19,i5,a14)
  150     format(7F10.6)
      endif
c
c     Initialize approximate solution.
c     If `sol_input' = 2 then an approximate solution is already stored in `Y'.
c     If `sol_input' = 1 then call `init_Y' to get an initial solution.
      if (sol_input.eq.1) then
          call init_Y(Nsub,neqns,mesh,Y,ldY)
      endif
c
c     Transfer the 2-D array, `Y', to its 1-D equivalent for use within MIRKDC.
      do 15 i = 1,Nsub
           do 20 j = 1,neqns
             Y(i*neqns + j) = Y(i*ldY + j)
 20        continue
 15   continue
c
c     Define the Newton tolerance and maximum iteration count
      newton_tol = tol/100.0d0
      maxiter = 40
c
c***main outer loop to control discrete problem sequence***
c   REPEAT-UNTIL LOOP  (Implemented here as a GOTO-IF pair)
 1                                                   continue
c
c     Compute array sizes and indexes to provide access to the work
c     array during the call to `newiter'. `s_defect' is also computed
c     here since its value is required in more than one section of the code.
      s_k_discrete = Mxs * neqns * Nsub
      s_k_interp = Mxs * neqns * Nsub
      i_k_discrete = 1
      i_k_interp = s_k_discrete + 1
      s_PHI = neqns * (Nsub+1)
      s_delta = neqns * (Nsub+1)
      s_top = neqns**2
      s_bot = neqns**2
      s_blocks = 2 * neqns**2 * Nsub
      s_defect = Nsub*neqns
      i_PHI = s_k_discrete + s_k_interp + 1
      i_delta = s_k_discrete + s_k_interp + s_PHI + 1
      i_top = s_k_discrete + s_k_interp + s_PHI +
     +          s_delta + 1
      i_bot = s_k_discrete + s_k_interp + s_PHI +
     +          s_delta + s_top + 1
      i_blocks = s_k_discrete + s_k_interp + s_PHI +
     +             s_delta + s_top + s_bot + 1
      i_work = s_k_discrete + s_k_interp + s_PHI +
     +           s_delta + s_top + s_bot + s_blocks + 1
c
c     Call `NewIter' to compute a solution for the discrete problem
c     based on the current mesh.
c
      if (profile .GT. print_level_0) then
          write(6,*) 'Begin the Newton iteration:'
          write(6,*)
      end if
c
      call NewIter(neqns,leftbc,Nsub,mesh,Y,
     +       newton_tol,maxiter,info,work(i_PHI),work(i_delta),
     +       work(i_top),work(i_bot),work(i_blocks),iwork,
     +       work(i_k_discrete),work(i_work),fsub,gsub,dfsub,dgsub)
c
      if ((info .NE. 0) .AND. (profile .EQ. print_level_1)) then
              write(*,*)'the Newton iteration has failed.'
      endif
      if ((info .NE. 0 ) .AND. (profile .GT. print_level_1)) then
          if (info .EQ. 1) then
                write(6,*)'the Newton iteration failed because the',
     +          ' maximum number of allowed iterations was exceeded.'
          endif
          if (info .EQ. 2) then
                write(6,*)'the Newton iteration failed because a ',
     +          'singular matrix was encountered in the computation.'
          endif
          if (info .EQ. 3) then
                write(6,*)'the Newton iteration failed because a ',
     +          'sufficiently large damping factor could not be found.'
          endif
      endif
c
c
c     If the Newton iteration has converged, compute the defect estimate.
c
      if (info .EQ. 0) then
c
c         Set indexes to work array.
          i_defect = s_k_discrete + s_k_interp + 1
          i_shift = s_k_discrete + s_k_interp + s_defect
          i_z = i_shift + 1
          i_z_prime = i_shift + neqns + 1
          i_def_1 = i_shift + 2*neqns + 1
          i_def_2 = i_shift + 3*neqns + 1
          i_f_1 = i_shift + 4*neqns + 1
          i_f_2 = i_shift + 5*neqns + 1
          i_temp_1 = i_shift + 6*neqns + 1
          i_temp_2 = i_shift + 7*neqns + 1
          i_work = i_shift + 8*neqns + 1
c
          call defect_estimate(neqns,Nsub,mesh,Y,
     +         work(i_defect), defect_norm,info,
     +         work(i_z),work(i_z_prime),
     +         work(i_def_1),work(i_def_2), work(i_f_1),
     +         work(i_f_2),work(i_temp_1),work(i_temp_2),
     +         work(i_k_discrete),work(i_k_interp),
     +         work(i_work),fsub)
c
c         The `defect_estimate' routine returns the defect estimate
c         vector along with its max norm. It also checks the size of
c         the defect norm to make sure it is less than a certain
c         threshold value (see the `defect_estimate' routine for
c         further details.) If the defect norm is too large the
c         `defect_estimate' routine sets `info' to 4 to indicate that
c         the defect is too large and that the solution should not be trusted.
c
          if (profile .GT. print_level_0) then
                write(6,*)'the Newton iteration was successful. ',
     +          'Build a continuous approximation to the solution and',
     +          ' sample the defect.'
                write(6,*)'The norm of defect is ',defect_norm
                if (info .EQ. 4) then
                  write(6,*)'Since the defect is greater than ',
     +            '10% the solution is not acceptable.'
                endif
          endif
c
      endif
c     End of `if-endif (info .EQ. 0)'
c
c
c     If the Newton iteration failed or the defect was not acceptable,
c     `info' will not be 0 and we will proceed to the `else' clause of the
c     following `if-then-else' statement where mesh halving will be
c     attempted. Otherwise, we will proceed to the `then' clause below
c     and attempt mesh redistribution.
c
c     Set work array indexes for `mesh_selector' and `half mesh'.
      s_mesh_new = MxNsub + 1
      i_mesh_new = s_k_discrete+s_k_interp
     +                                          +s_defect+1
      i_work = s_k_discrete+s_k_interp+s_defect
     +                                          +s_mesh_new+1
      i_Y_new = s_k_discrete+s_k_interp
     +                           +s_defect+s_mesh_new+1
c
c
      if (info.EQ.0) then
c
c         If we have not yet satisfied the tolerance, we must try again
c         on a new mesh, if possible.
c
          if (defect_norm .GT. tol) then
c
            if (profile .GT. print_level_0) then
                  write(6,*) 'User defined tolerance',tol,' has not '
                  write(6,*) 'been satisfied.'
                  write(6,*) 'Construct a new mesh which ',
     +            'equidistributes the defect.'
            endif
c
c           Attempt to select a new mesh by calling `mesh_selector'.
            call mesh_selector(neqns,Nsub,mesh,work(i_defect),
     +         tol,Nsub_star,work(i_mesh_new),
     +           info,MxNsub,work(i_work))
c
c           If we were not successful in obtaining a new mesh, `info'
c           will not equal 0 and we should terminate the computation.
c           Otherwise we can compute a new solution estimate on the new
c           mesh, and then update `Nsub', `mesh' and `Y'.
c
            if (info .EQ. 0) then
              if (profile .GT. print_level_0) then
                  write(6,*)
                  write(6,*)
                  write(6,160)'New mesh will be of size ',
     +                         Nsub_star,' subintervals.'
 160              format(1x,a25,i5,a14)
              end if
              if (profile .GT. print_level_1) then
                   write(6,150)(work(i_mesh_new+i), i=0, Nsub_star)
                   write(6,*)
              end if
c
c             Use current computed solution to generate initial guess for next
c             problem. Compute solution approximations at new mesh points.
c
c             Set work array indexes.
              s_Y_new = neqns*(Nsub_star+1)
              i_Y_new = s_k_discrete+s_k_interp
     +                              +s_defect+s_mesh_new+1
              i_work = s_k_discrete+s_k_interp
     +                        +s_defect+s_mesh_new+s_Y_new+1
c
              do 30 i = 0, Nsub_star
                  call interp_eval(neqns,Nsub,mesh,Y,
     +              work(i_mesh_new+i),
     +              work(i_Y_new+i*neqns), work(i_work),
     +              work(i_k_discrete),work(i_k_interp))
 30           continue
c
c             Copy `mesh_new' to `mesh', `Y_new' to `Y'; update `Nsub'.
              call dcopy(Nsub_star+1,work(i_mesh_new),1,mesh,1)
              call dcopy((Nsub_star+1)*neqns,work(i_Y_new),1,Y,1)
              Nsub = Nsub_star
c
            endif
c           End of `if-endif (info .EQ. 0)'
c
          endif
c         End of `if-endif (defect_norm .GT. tol)'
c
          if ( (info. EQ. 0) .AND. (defect_norm .LT. tol) .AND.
     +         (profile .GT. print_level_0) ) then
                  write(6,*)'The user defined tolerance',tol,' has ',
     +            'been satisfied.'
                  write(6,*)'Successful computation. '
                  write(6,*)
          endif
c
      else
c     Else clause associated with `if-then-else (info .EQ. 0)'
c
          if (profile .GT. print_level_0) then
            write(6,*)'Cannot obtain a solution for the current mesh.'
          endif
c
c         Simply half the current mesh, if possible, and try again.
          if (2*Nsub .GT. MxNsub) then
c              New mesh would be too large. Terminate the computation.
               if (profile .GT. print_level_0) then
                  write(6,*)'New mesh would be too large.'
               end if
               info = -1
          else
c           Else clause associated with if_then_else `2*Nsub .GT. MxNsub'
c
            call half_mesh(Nsub,mesh,work(i_mesh_new))
            Nsub_star = 2*Nsub
c
c           Copy `mesh_new' to `mesh'
            call dcopy(Nsub_star+1,work(i_mesh_new),1,mesh,1)
            Nsub = Nsub_star
c
            if (profile .GT. print_level_0) then
                  write(6,*)'Double the mesh and try again.'
                  write(6,*)
                  write(6,*)
                  write(6,170)'New mesh will be of size', Nsub,'.'
  170             format(1x,a24,i5,a1)
            end if
            if (profile .GT. print_level_1) then
                   write(6,150)(mesh(i), i = 0, Nsub)
                   write(6,*)
            end if
c
c           Initialize the approximate solution - store in `Y_new'.
            if (sol_input .EQ. 1) then
                  call init_Y(Nsub,neqns,mesh,Y,ldY)
c                 Transferring the 2-D array, Y, to its 1-D equivalent for
c                 use with in MIRKDC.
                  do 35 i = 1,Nsub
                    do 40 j = 1,neqns
                      Y(i*neqns + j) = Y(i*ldY + j)
  40                continue
  35              continue
            endif
c           If sol_input = 2 then the initial guess on the initial mesh
c           was provided through `Y' and `mesh', respectively. Since this
c           initial guess is no longer usable, with the new mesh, we
c           continue by simply setting `Y_new' to zero.
            if (sol_input .EQ. 2) then
                  if (profile .GT. print_level_0) then
                    write(6,*)'Iteration failed and there is no ',
     +                    'Init_Y routine; Y will be set to zero.'
                    write(6,*)'It would be better to retry the ',
     +                 'continuation sequence with smaller steps'
                  endif
                  call dcopy((Nsub+1)*neqns,0.0d0,0,Y,1)
            endif
c
c           Reset `info' to 0 and `defect_norm' > `tol'
c           to force a restart. (See UNTIL conditions below.)
            info = 0
            defect_norm = 2* tol
c
          endif
c         End of `if-then-else (2*Nsub .GT. MxNsub)'
c
      endif
c     End of `if-then-else (info .EQ. 0)', `Acceptable solution obtained'
c
c   UNTIL( (info .NE. 0) .OR. (defect_norm .LE. tol))
                                       if((info .EQ. 0) .AND.
     +                                 (defect_norm .GT. tol))goto1
c
c***conclusion***
c
c     Set values of integer work array to allow user convenient
c     use of `sol_eval' routine, after a successful termination.
c     See the `sol_eval' routine for further details.
      if (info .EQ. 0) then
        iwork(1) = neqns
        iwork(2) = Nsub
        iwork(3) = s
        iwork(4) = s_star
        iwork(5) = method
c
c       Set values of double precision work array to allow user
c       convenient use of `sol_eval' routine after a successful termination.
        call dcopy(Nsub+1,mesh(0),1,work(2*Mxs*neqns*Nsub+1),1)
        call dcopy((Nsub+1)*neqns,Y,1,
     +                 work(2*Mxs*neqns*Nsub+(Nsub+1)+1),1)
c
c       Transfer the 1-D version of `Y' back to its 2-D equivalent.
        do 45 i = 1,(Nsub-1),-1
          do 50 j = 1,neqns
           Y(i*ldY + j) = Y(i*neqns + j)
 50       continue
 45     continue
      else
        if (profile.GT.print_level_0) then
          write(6,*)
          write(6,*) 'The computation has failed.'
          write(6,*)
          write(6,*) '1) If the profiling information begins with and ',
     +    'simply proceeds through a series of failed Newton ',
     +    'iterations each followed by a doubling of the mesh, then',
     +    'there may be an error in the coding of the problem.',
     +    'It is recommended that the routines fsub, gsub, dfsub and',
     +    ' dgsub be checked.'
          write(6,*)
          write(6,*) '2) If the profiling information shows a series',
     +    'of successful, converged Newton iterations with a halt',
     +    'simply because the value of MxNsub was exceeded, then it',
     +    'is recommended that the program be run again with a larger',
     +    ' value of MxNsub.'
          write(6,*)
          write(6,*) '3) More information may be obtained by choosing',
     +    'a larger value for "output_control".'
        endif
      endif
c
      return
      end

       SUBROUTINE COLROW(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,
     *             NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,B,X,IFLAG)
C
C***************************************************************
C
C  THIS PROGRAM SOLVES THE LINEAR SYSTEM  A*X = B  WHERE  A IS
C  AN ALMOST BLOCK DIAGONAL MATRIX OF THE FORM
C
C               TOPBLK
C               ARRAY(1)
C                     ARRAY(2)
C                          .
C                             .
C                                .
C                                   .
C                                    ARRAY(NBLOKS)
C                                           BOTBLK
C
C  WHERE
C           TOPBLK IS  NRWTOP  BY NOVRLP
C           ARRAY(K), K=1,NBLOKS, ARE NRWBLK BY NRWBLK+NOVRLP
C           BOTBLK IS NRWBOT BY NOVRLP,
C  AND
C           NOVRLP = NRWTOP + NRWBOT
C  WITH
C           NOVRLP.LE.NRWBLK .
C
C  THE LINEAR SYSTEM IS OF ORDER  N = NBLOKS*NRWBLK + NOVRLP.
C
C  THE METHOD IMPLEMENTED IS BASED ON GAUSS ELIMINATION WITH
C  ALTERNATE ROW AND COLUMN ELIMINATION WITH PARTIAL PIVOTING,
C  WHICH PRODUCES A STABLE DECOMPOSITION OF THE MATRIX  A
C  WITHOUT INTRODUCING FILL-IN.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  TO OBTAIN A SINGLE PRECISION VERSION OF THIS PACKAGE, REMOVE
C  ALL DOUBLE PRECISION STATEMENTS.  THERE IS ONE SUCH STATEMENT
C  IN C O L R O W, THREE IN C R D C M P, AND TWO IN C R S O L V.
C  IN ADDITION, REFERENCES TO BUILT-IN FUNCTIONS DABS AND DMAX1
C  MUST BE REPLACED BY ABS AND AMAX1, RESPECTIVELY.  DABS OCCURS
C  NINE TIMES, IN C R D C M P.  DMAX1 OCCURS FOUR TIMES, IN
C  C R D C M P.  FINALLY, ZERO IS INITIALISED TO 0.D0 IN A
C  DATA STATEMENT IN C R D C M P.  THIS MUST BE REPLACED BY:
C               DATA ZERO/0.0/
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C               N      - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM,
C                         GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         THE FIRST BLOCK OF THE ALMOST BLOCK
C                         DIAGONAL MATRIX A
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C                         BY NCLBLK BLOCK OF THE MATRIX A
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         THE LAST BLOCK OF THE MATRIX A
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C               TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C                        DESIRED DECOMPOSITION OF THE MATRIX A
C                        (IF IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C                    X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR (IF IFLAG = 0)
C
C               IFLAG  - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  AUXILIARY PROGRAMS  *****
C
C       CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
C    *     BOTBLK,NRWBOT,PIVOT,IFLAG)
C            - DECOMPOSES THE MATRIX  A  USING MODIFIED
C              ALTERNATE ROW AND COLUMN ELIMINATON WITH
C              PARTIAL PIVOTING, AND IS USED FOR THIS
C              PURPOSE IN C O L R O W.
C              THE ARGUMENTS ARE AS IN C O L R O W.
C
C       CRSLVE(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
C    *     BOTBLK,NRWBOT,PIVOT,B,X)
C            - SOLVES THE SYSTEM A*X = B ONCE A IS DECOMPOSED.
C              THE ARGUMENTS ARE ALLAS IN C O L R O W.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       THE SUBROUTINE  C O L R O W  AUTOMATICALLY SOLVES THE
C  INPUT SYSTEM WHEN IFLAG=0.  C O L R O W  IS CALLED ONLY ONCE
C  FOR A GIVEN SYSTEM. THE SOLUTION FOR A SEQUENCE OF P RIGHT
C  HAND SIDES CAN BE OBTAINED BY ONE CALL TO  C O L R O W  AND
C  P-1 CALLS TO CRSLVE ONLY. SINCE THE ARRAYS TOPBLK,ARRAY,
C  BOTBLK AND PIVOT CONTAIN THE DECOMPOSITION OF THE GIVEN
C  COEFFICIENT MATRIX AND PIVOTING INFORMATION ON RETURN FROM
C  C O L R O W , THEY MUST NOT BE ALTERED BETWEEN SUCCESSIVE
C  CALLS TO CRSLVE WITH THE SAME LEFT HAND SIDES. FOR THE
C  SAME REASON, IF THE USER WISHES TO SAVE THE COEFFICIENT
C  MATRIX, THE ARRAYS TOPBLK,ARRAY,BOTBLK MUST BE COPIED
C  BEFORE A CALL TO  C O L R O W .
C
C*************************************************************************
C
C             *****  SAMPLE CALLING PROGRAM  *****
C
C      THE FOLLOWING PROGRAM WILL EXERCISE COLROW, IN THE
C      CASE WHEN THE COEFFICIENT MATRIX IS NON-SINGULAR.
C
C       DOUBLE PRECISION TOP,AR,BOT,B,X
C       DOUBLE PRECISION ERROR,ERR
C       DIMENSION TOP(2,4),AR(4,8,2),BOT(2,4),B(12),X(12)
C       INTEGER PIVOT(12)
C       DATA N,NRWTOP,NOVRLP,NRWBLK,NCLBLK,NBLOKS,NRWBOT/12,2,4,4,8,2,2/
C       DATA TOP(1,1),TOP(1,2),TOP(1,3),TOP(1,4),
C    *       TOP(2,1),TOP(2,2),TOP(2,3),TOP(2,4)/
C    *0.0 D0,-0.98D0,-0.79D0,-0.15D0,
C    *-1.00D0, 0.25D0,-0.87D0, 0.35D0/
C       DATA AR(1,1,1),AR(1,2,1),AR(1,3,1),AR(1,4,1),
C    *       AR(1,5,1),AR(1,6,1),AR(1,7,1),AR(1,8,1)/
C    *0.78D0, 0.31D0,-0.85D0, 0.89D0,-0.69D0,-0.98D0,-0.76D0,-0.82D0/
C       DATA AR(2,1,1),AR(2,2,1),AR(2,3,1),AR(2,4,1),
C    *       AR(2,5,1),AR(2,6,1),AR(2,7,1),AR(2,8,1)/
C    *0.12D0,-0.01D0, 0.75D0, 0.32D0,-1.00D0,-0.53D0,-0.83D0,-0.98D0/
C       DATA AR(3,1,1),AR(3,2,1),AR(3,3,1),AR(3,4,1),
C    *       AR(3,5,1),AR(3,6,1),AR(3,7,1),AR(3,8,1)/
C    *-0.58D0, 0.04D0, 0.87D0, 0.38D0,-1.00D0,-0.21D0,-0.93D0,-0.84D0/
C       DATA AR(4,1,1),AR(4,2,1),AR(4,3,1),AR(4,4,1),
C    *       AR(4,5,1),AR(4,6,1),AR(4,7,1),AR(4,8,1)/
C    *-0.21D0,-0.91D0,-0.09D0,-0.62D0,-1.99D0,-1.12D0,-1.21D0, 0.07D0/
C       DATA AR(1,1,2),AR(1,2,2),AR(1,3,2),AR(1,4,2),
C    *       AR(1,5,2),AR(1,6,2),AR(1,7,2),AR(1,8,2)/
C    *0.78D0,-0.93D0,-0.76D0, 0.48D0,-0.87D0,-0.14D0,-1.00D0,-0.59D0/
C       DATA AR(2,1,2),AR(2,2,2),AR(2,3,2),AR(2,4,2),
C    *       AR(2,5,2),AR(2,6,2),AR(2,7,2),AR(2,8,2)/
C    *-0.99D0, 0.21D0,-0.73D0,-0.48D0,-0.93D0,-0.91D0, 0.10D0,-0.89D0/
C       DATA AR(3,1,2),AR(3,2,2),AR(3,3,2),AR(3,4,2),
C    *       AR(3,5,2),AR(3,6,2),AR(3,7,2),AR(3,8,2)/
C    *-0.68D0,-0.09D0,-0.58D0,-0.21D0, 0.85D0,-0.39D0, 0.79D0,-0.71D0/
C       DATA AR(4,1,2),AR(4,2,2),AR(4,3,2),AR(4,4,2),
C    *       AR(4,5,2),AR(4,6,2),AR(4,7,2),AR(4,8,2)/
C    *0.39D0,-0.99D0,-0.12D0,-0.75D0,-0.68D0,-0.99D0, 0.50D0,-0.88D0/
C       DATA BOT(1,1),BOT(1,2),BOT(1,3),BOT(1,4),
C    *       BOT(2,1),BOT(2,2),BOT(2,3),BOT(2,4)/
C    *0.71D0,-0.64D0, 0.0 D0, 0.48D0,
C    *0.08D0,100.0D0,50.00D0,15.00D0/
C       DATA B(1),B(2),B(3),B(4),B(5),B(6),B(7),B(8),B(9),B(10),B(11),
C    *       B(12)/
C    *-1.92D0,-1.27D0,-2.12D0,-2.16D0,-2.27D0,-6.08D0,-3.03D0,
C    *-4.62D0,-1.02D0,-3.52D0,.55D0,165.08D0/
C
C*************************************************************************
C
C   THE INPUT MATRIX IS GIVEN BY:
C
C  0.0  -0.98 -0.79 -0.15
C -1.00  0.25 -0.87  0.35
C  0.78  0.31 -0.85  0.89 -0.69 -0.98 -0.76 -0.82
C  0.12 -0.01  0.75  0.32 -1.00 -0.53 -0.83 -0.98
C -0.58  0.04  0.87  0.38 -1.00 -0.21 -0.93 -0.84
C -0.21 -0.91 -0.09 -0.62 -1.99 -1.12 -1.21  0.07
C                          0.78 -0.93 -0.76  0.48 -0.87 -0.14 -1.00 -0.59
C                         -0.99  0.21 -0.73 -0.48 -0.93 -0.91  0.10 -0.89
C                         -0.68 -0.09 -0.58 -0.21  0.85 -0.39  0.79 -0.71
C                          0.39 -0.99 -0.12 -0.75 -0.68 -0.99  0.50 -0.88
C                                                  0.71 -0.64  0.0   0.48
C                                                  0.08 100.0 50.00 15.00
C
C       THE RIGHT HAND SIDE IS GIVEN BY:
C
C         B = (-1.92,-1.27,-2.12,-2.16,-2.27,-6.08,-3.03,-4.62,
C              -1.02,-3.52,0.55,165.08)
C
C       THE SOLUTION OF THIS SYSTEM IS GIVEN BY;
C
C          X = (1,1,1,1,1,1,1,1,1,1,1,1)
C
C*************************************************************************
C
C       CALL COLROW(N,TOP,NRWTOP,NOVRLP,AR,NRWBLK,NCLBLK,NBLOKS,
C    *              BOT,NRWBOT,PIVOT,B,X,IFLAG)
C       IF(IFLAG.NE.0)GO TO 1000
C       ERROR = 0.D0
C       DO 10 I=1,N
C          ERR = 1.D0 - X(I)
C          ERROR = DMAX1(ERROR,DABS(ERR))
C          WRITE(6,100)X(I),ERR
C  10   CONTINUE
C       WRITE(6,200)ERROR
C 200   FORMAT(12H MAX ERROR = ,D15.7)
C 100   FORMAT(1H ,F15.7,D15.7)
C       RETURN
C1000   CONTINUE
C       WRITE(6,300)IFLAG
C 300   FORMAT(9H IFLAG =  ,I3)
C       RETURN
C       END
C
C***************************************************************
C
        DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK,B,X
        INTEGER PIVOT(1)
        DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),
     *          BOTBLK(NRWBOT,1),B(1),X(1)
        CALL CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
     *          BOTBLK,NRWBOT,PIVOT,IFLAG)
        IF(IFLAG.NE.0)RETURN
        CALL CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS,
     *          BOTBLK,NRWBOT,PIVOT,B,X)
        RETURN
        END
        SUBROUTINE CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,
     *             NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,IFLAG)
C
C***************************************************************
C
C  C R D C M P DECOMPOSES THE ALMOST BLOCK DIAGONAL MATRIX A
C  USING MODIFIED ALTERNATE ROW AND COLUMN ELIMINATION WITH
C  PARTIAL PIVOTING.  THE MATRIX  A  IS STORED IN THE ARRAYS
C  TOPBLK, ARRAY, AND BOTBLK.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C               N      - INTEGER
C                         THE ORDER OF THE LINEAR SYSTEM,
C                         GIVEN BY NBLOKS*NRWBLK + NOVRLP
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         THE FIRST BLOCK OF THE ALMOST BLOCK
C                         DIAGONAL MATRIX A TO BE DECOMPOSED
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         ARRAY(,,K) CONTAINS THE K-TH NRWBLK
C                         BY NCLBLK BLOCK OF THE MATRIX A
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         THE LAST BLOCK OF THE MATRIX A
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C               TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE
C                        DESIRED DECOMPOSITION OF THE MATRIX A
C                        (IF IFLAG = 0)
C
C                PIVOT - INTEGER(N)
C                         RECORDS THE PIVOTING INDICES DETER-
C                         MINED IN THE DECOMPOSITION
C
C               IFLAG  - INTEGER
C                         =  1, IF INPUT PARAMETERS ARE INVALID
C                         = -1, IF MATRIX IS SINGULAR
C                         =  0, OTHERWISE
C
C***************************************************************
C
        DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK
        DOUBLE PRECISION ROWMAX,ROWPIV,ROWMLT,COLMAX,COLPIV
        DOUBLE PRECISION SWAP,COLMLT,PIVMAX,ZERO,TEMPIV
        INTEGER PIVOT(1)
        DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),
     *          BOTBLK(NRWBOT,1)
        DATA ZERO/0.0D0/
C
C***************************************************************
C
C          ****  DEFINE THE CONSTANTS USED THROUGHOUT  ****
C
C***************************************************************
C
        IFLAG = 0
        PIVMAX = ZERO
        NRWTP1 = NRWTOP+1
        NROWEL = NRWBLK-NRWTOP
        NRWEL1 = NROWEL+1
        NVRLP0 = NOVRLP-1
C
C***************************************************************
C
C          ****  CHECK VALIDITY OF THE INPUT PARAMETERS....
C
C               IF PARAMETERS ARE INVALID THEN TERMINATE AT 10;
C                                         ELSE CONTINUE AT 100.
C
C***************************************************************
C
        IF(N.NE.NBLOKS*NRWBLK+NOVRLP)GO TO 10
        IF(NOVRLP.NE.NRWTOP+NRWBOT)GO TO 10
        IF(NCLBLK.NE.NOVRLP+NRWBLK)GO TO 10
        IF(NOVRLP.GT.NRWBLK)GO TO 10
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE ACCEPTABLE - CONTINUE AT 100.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        GO TO 100
10      CONTINUE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          PARAMETERS ARE INVALID.  SET IFLAG = 1, AND TERMINATE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IFLAG = 1
        RETURN
100     CONTINUE
C
C***************************************************************
C
C               ****  FIRST, IN TOPBLK....
C
C***************************************************************
C
C          ***  APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN
C                 PIVOTING ....
C
C***************************************************************
C
        DO 190 I = 1,NRWTOP
           IPLUS1 = I+1
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           IPVT = I
           COLMAX = DABS(TOPBLK(I,I))
           DO 110 J = IPLUS1,NOVRLP
              TEMPIV = DABS(TOPBLK(I,J))
              IF(TEMPIV.LE.COLMAX)GO TO 110
                 IPVT = J
                 COLMAX = TEMPIV
110        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           IF(PIVMAX+COLMAX.EQ.PIVMAX)GO TO 1000
           PIVMAX = DMAX1(COLMAX,PIVMAX)
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE COLUMNS
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           PIVOT(I) = IPVT
           IF(IPVT.EQ.I)GO TO 140
              DO 120 L = I,NRWTOP
                 SWAP = TOPBLK(L,IPVT)
                 TOPBLK(L,IPVT) = TOPBLK(L,I)
                 TOPBLK(L,I) = SWAP
120           CONTINUE
              DO 130 L = 1,NRWBLK
                 SWAP = ARRAY(L,IPVT,1)
                 ARRAY(L,IPVT,1) = ARRAY(L,I,1)
                 ARRAY(L,I,1) = SWAP
130           CONTINUE
140        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS AND PERFORM COLUMN
C                       ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           COLPIV = TOPBLK(I,I)
           DO 180 J = IPLUS1,NOVRLP
              COLMLT = TOPBLK(I,J)/COLPIV
              TOPBLK(I,J) = COLMLT
              IF(IPLUS1.GT.NRWTOP)GO TO 160
                 DO 150 L = IPLUS1,NRWTOP
                    TOPBLK(L,J) = TOPBLK(L,J)-COLMLT*TOPBLK(L,I)
150              CONTINUE
160           CONTINUE
              DO 170 L = 1,NRWBLK
                 ARRAY(L,J,1) = ARRAY(L,J,1)-COLMLT*ARRAY(L,I,1)
170           CONTINUE
180        CONTINUE
190     CONTINUE
C
C***************************************************************
C
C          ****  IN EACH BLOCK ARRAY(,,K)....
C
C***************************************************************
C
        INCR = 0
        DO 395 K = 1,NBLOKS
           KPLUS1 = K+1
C
C          *****************************************************
C
C          ***  FIRST APPLY NRWBLK-NRWTOP ROW ELIMINATIONS WITH
C                       ROW PIVOTING....
C
C          *****************************************************
C
           DO 270 J = NRWTP1,NRWBLK
              JPLUS1 = J+1
              JMINN = J-NRWTOP
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE ROW PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IPVT = JMINN
              ROWMAX = DABS(ARRAY(JMINN,J,K))
              LOOP = JMINN+1
              DO 210 I = LOOP,NRWBLK
                 TEMPIV = DABS(ARRAY(I,J,K))
                 IF(TEMPIV.LE.ROWMAX)GO TO 210
                 IPVT = I
                 ROWMAX = TEMPIV
210           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IF(PIVMAX+ROWMAX.EQ.PIVMAX)GO TO  1000
              PIVMAX = DMAX1(ROWMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE ROWS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              INCRJ = INCR+J
              PIVOT(INCRJ) = INCR+IPVT+NRWTOP
              IF(IPVT.EQ.JMINN)GO TO 230
                 DO 220 L = J,NCLBLK
                    SWAP = ARRAY(IPVT,L,K)
                    ARRAY(IPVT,L,K) = ARRAY(JMINN,L,K)
                    ARRAY(JMINN,L,K) = SWAP
220              CONTINUE
230           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLERS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              ROWPIV = ARRAY(JMINN,J,K)
              DO 240 I = LOOP,NRWBLK
                 ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV
240           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              DO 260 L = JPLUS1,NCLBLK
                 ROWMLT = ARRAY(JMINN,L,K)
                 DO 250 I = LOOP,NRWBLK
                    ARRAY(I,L,K) = ARRAY(I,L,K)
     *                                -ROWMLT*ARRAY(I,J,K)
250              CONTINUE
260           CONTINUE
270        CONTINUE
C
C          *****************************************************
C
C          ***  NOW APPLY NRWTOP COLUMN ELIMINATIONS WITH
C                      COLUMN PIVOTING....
C
C          *****************************************************
C
           DO 390 I = NRWEL1,NRWBLK
              IPLUSN = I+NRWTOP
              IPLUS1 = I+1
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE COLUMN PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IPVT = IPLUSN
              COLMAX = DABS(ARRAY(I,IPVT,K))
              LOOP = IPLUSN+1
              DO 310 J = LOOP,NCLBLK
                 TEMPIV = DABS(ARRAY(I,J,K))
                 IF(TEMPIV.LE.COLMAX)GO TO 310
                 IPVT = J
                 COLMAX = TEMPIV
310           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IF(PIVMAX+COLMAX.EQ.PIVMAX)GO TO 1000
              PIVMAX = DMAX1(COLMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE COLUMNS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              INCRN = INCR+IPLUSN
              PIVOT(INCRN) = INCR+IPVT
              IRWBLK = IPLUSN-NRWBLK
              IF(IPVT.EQ.IPLUSN)GO TO 340
                 DO 315 L = I,NRWBLK
                    SWAP = ARRAY(L,IPVT,K)
                    ARRAY(L,IPVT,K) = ARRAY(L,IPLUSN,K)
                    ARRAY(L,IPLUSN,K) = SWAP
315              CONTINUE
                 IPVBLK = IPVT-NRWBLK
                 IF(K.EQ.NBLOKS)GO TO 330
                    DO 320 L = 1,NRWBLK
                       SWAP = ARRAY(L,IPVBLK,KPLUS1)
                       ARRAY(L,IPVBLK,KPLUS1)
     *                                 = ARRAY(L,IRWBLK,KPLUS1)
                       ARRAY(L,IRWBLK,KPLUS1) = SWAP
320                 CONTINUE
                    GO TO 340
330              CONTINUE
                 DO 335 L = 1,NRWBOT
                    SWAP = BOTBLK(L,IPVBLK)
                    BOTBLK(L,IPVBLK) = BOTBLK(L,IRWBLK)
                    BOTBLK(L,IRWBLK) = SWAP
335              CONTINUE
340           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS AND PERFORM COLUMN
C                       ELIMINATION
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              COLPIV = ARRAY(I,IPLUSN,K)
              DO 380 J = LOOP,NCLBLK
                 COLMLT = ARRAY(I,J,K)/COLPIV
                 ARRAY(I,J,K) = COLMLT
                 IF(I.EQ.NRWBLK)GO TO 350
                    DO 345 L = IPLUS1,NRWBLK
                       ARRAY(L,J,K) = ARRAY(L,J,K)
     *                                -COLMLT*ARRAY(L,IPLUSN,K)
345                 CONTINUE
350              CONTINUE
                 JRWBLK = J-NRWBLK
                 IF(K.EQ.NBLOKS)GO TO 370
                    DO 360 L = 1,NRWBLK
                       ARRAY(L,JRWBLK,KPLUS1) =
     *                                  ARRAY(L,JRWBLK,KPLUS1)
     *                         -COLMLT*ARRAY(L,IRWBLK,KPLUS1)
360                 CONTINUE
                    GO TO 380
370              CONTINUE
                 DO 375 L = 1,NRWBOT
                    BOTBLK(L,JRWBLK) = BOTBLK(L,JRWBLK)
     *                              -COLMLT*BOTBLK(L,IRWBLK)
375              CONTINUE
380           CONTINUE
390        CONTINUE
           INCR = INCR + NRWBLK
395     CONTINUE
C
C***************************************************************
C
C          ****  FINALLY, IN BOTBLK....
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          ***  APPLY NRWBOT ROW ELIMINATIONS WITH ROW
C                  PIVOTING....
C
C               IF BOT HAS JUST ONE ROW GO TO 500
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF(NRWBOT.EQ.1)GO TO 500
           DO 470 J = NRWTP1,NVRLP0
              JPLUS1 = J+1
              JMINN = J-NRWTOP
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               DETERMINE ROW PIVOT AND PIVOT INDEX
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IPVT = JMINN
              ROWMAX = DABS(BOTBLK(JMINN,J))
              LOOP = JMINN+1
              DO 410 I = LOOP,NRWBOT
                 TEMPIV = DABS(BOTBLK(I,J))
                 IF(TEMPIV.LE.ROWMAX) GO TO 410
                 IPVT = I
                 ROWMAX = TEMPIV
410           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               TEST FOR SINGULARITY:
C
C                       IF SINGULAR THEN TERMINATE AT 1000;
C                                   ELSE CONTINUE.
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              IF(PIVMAX+ROWMAX.EQ.PIVMAX)GO TO 1000
              PIVMAX = DMAX1(ROWMAX,PIVMAX)
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               IF NECESSARY INTERCHANGE ROWS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              INCRJ = INCR+J
              PIVOT(INCRJ) = INCR+IPVT+NRWTOP
              IF(IPVT.EQ.JMINN)GO TO 430
                 DO 420 L = J,NOVRLP
                    SWAP = BOTBLK(IPVT,L)
                    BOTBLK(IPVT,L) = BOTBLK(JMINN,L)
                    BOTBLK(JMINN,L) = SWAP
420              CONTINUE
430           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               COMPUTE MULTIPLIERS
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              ROWPIV = BOTBLK(JMINN,J)
              DO 440 I = LOOP,NRWBOT
                 BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV
440           CONTINUE
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               PERFORM ROW ELIMINATION WITH COLUMN INDEXING
C
C             CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
              DO 460 L = JPLUS1,NOVRLP
                 ROWMLT = BOTBLK(JMINN,L)
                 DO 450 I = LOOP,NRWBOT
                    BOTBLK(I,L) = BOTBLK(I,L)-ROWMLT*BOTBLK(I,J)
450              CONTINUE
460           CONTINUE
470        CONTINUE
500     CONTINUE
C
C***************************************************************
C
C          DONE PROVIDED THE LAST ELEMENT IS NOT ZERO
C
C***************************************************************
C
        IF(PIVMAX+DABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN
C
C***************************************************************
C
C       ****  MATRIX IS SINGULAR - SET IFLAG = - 1.
C                                  TERMINATE AT 1000.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
1000    CONTINUE
        IFLAG = -1
        RETURN
        END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine criterion(neqns, leftbc, Nsub, mesh, Y, top,blocks,
     +           bot, pivot, PHI, delta, g,k_discrete,work,fsub,gsub)
c
c***overview***
c
c        This routine takes as input the current iterate, `Y' and computes
c     a related function known as the natural criterion function, `g' = g(Y) =
c     0.5*|inv(J)*PHI(Y)|**2, where PHI(Y) is the residual function evaluated
c     at `Y' and J is the Jacobian of the residual function evaluated at `Y',
c     or possibly at a previous value of the iterate if a Fixed Jacobian
c     iteration is employed. See [Ascher, Mattheij, and Russell, 1988,
c     Pg. 335]. The intermediate quantities that arise during the calculation
c     of `g', namely, `PHI = PHI(Y)' and `delta = inv(J)*PHI(Y)', are returned
c     because they can be useful during the next step of the iteration.
c     The mesh of `Nsub' subintervals is stored in `mesh'; `neqns' is the
c     number of differential equations and the factored Jacobian is provided
c     in the arrays `top', `blocks', `bot', and `pivot'. The number of boundary
c     conditions imposed at the left endpoint is `leftbc'. The stages values
c     generated by the `resid' routine are stored in the `k_discrete'
c     array are passed out of this routine for use by other parts of the code.
c     `work' is a work array passed to the `resid' routine and used here for
c     temporary storage. The array, `work', is a work array.
c     The subroutines, `fsub', and `gsub', are provided by the user to define
c     the boundary value ODE and the boundary conditions.
c------------------------------------------------------------------------------
c***declaration of constants***
      integer Mxs
      parameter(Mxs=10)
c
c     `Mxs' Maximum value for the total number of stages of the RK method.
c
c     We wish to be able to detect when the value of the criterion
c     function will overflow. Thus we set the overflow guard to
c     sqrt(8.5d37)=9.22d18.
c
      double precision overflow_guard
      parameter(overflow_guard=9.22d18)
c
c***declaration of parameters***
c   imports:
      integer                 neqns,leftbc,Nsub
      double precision        mesh(0:Nsub)
      double precision        Y(neqns*(Nsub+1))
c
      double precision        top(neqns*leftbc)
      double precision        blocks(2*neqns**2 * Nsub)
      double precision        bot(neqns*(neqns-leftbc))
      integer                 pivot(neqns*(Nsub+1))
c
c     `neqns' Number of differential equations
c     `leftbc' Number of boundary conditions at the left
c              end of the problem interval
c     `Nsub' Number of subintervals into which the problem
c            interval is partitioned
c     `mesh' Set of points which partition the problem interval
c     `Y'    Current discrete approximation to the solution,
c            evaluated at each meshpoint.
c     `top' Derivative information for left boundary conditions,
c           in factored form.
c     `blocks' Derivative information for residual function, in factored form.
c     `bot' Derivative information for right boundary conditions,
c           in factored form.
c     `pivot' Contains pivoting information associated with the factorization
c             of the Newton matrix.  Together `top', `blocks', `bot', and
c             `pivot' contain the almost block diagonal matrix representing
c             the coefficient matrix of the Newton system in factored form
c             subsequent to the call to `COLROW'.
c   exports:
      double precision      PHI(neqns*(Nsub+1)), delta(neqns*(Nsub+1))
      double precision      g, k_discrete(Mxs*neqns*Nsub)
c
c     `PHI' Value of the residual function evaluated at `Y' .
c     `delta' Solution to the linear system with the usual Newton matrix as its
c             coefficient matrix and the righthand side, `PHI'. That is,
c             `delta = inv(J)*PHI', where J is the Jacobian of the residual
c             function.
c     `g' Value of natural criterion function at `Y'.
c                        i.e. g = 0.5*|inv(J)*PHI(Y)|^2
c     `k_discrete' Vector containing the discrete stages for
c                  all subintervals. The ith set of `s*neqns'
c                  locations of `k_discrete' contains the `s'
c                  stages, each of length `neqns' corresponding
c                  to the ith subinterval. Passed from `resid'.
c
c***declaration of work space***
      double precision      work(neqns*(Nsub+1))
c
c     `work' Work array passed to `resid' and also used here for
c            a temporary copy of the contents of `PHI' provided
c            to `crslve' because `crslve' will overwrite the
c            contents of the `PHI' vector, which we need to export.
c
c***user-supplied subroutines***
      external fsub,gsub
c
c     `fsub' Defines f(t,y) for first order system of differential equations.
c     `gsub' Defines the boundary condition equations.
c
c***declaration of variables for /IOCNTRL/ common block***
      integer   print_level_0, print_level_1, print_level_2
      parameter (print_level_0 = 0, print_level_1 = 1,
     *              print_level_2 = 2)
      integer                 profile
      common /IOCNTRL/        profile
c
c     `print_level_0' No output.
c     `print_level_1' Intermediate output.
c     `print_level_2' Full output.
c     `profile' Controls output, to standard output, of profiling infor-
c               mation such as Newton iteration counts, mesh selection,
c               relative defect estimate. This variable is identical to
c               the `output_control' parameter.
c----------------------------------------------------------------------------
c      called by: `fixed_jacob', `damp_factor'
c      calls to: `resid', `dcopy', `crslve', `idamax'
       integer                  idamax
c-----------------------------------------------------------------------------
c***evaluate natural criterion function at Y ****
c     First we compute PHI(Y) by calling `resid'; the value will
c     be returned in the vector `PHI'.
c
      call resid(neqns,leftbc,Nsub,mesh,Y,PHI,
     +                                 k_discrete,work,fsub,gsub)
c
c     We next compute `delta = inv(J)*PHI(Y)' by solving
c     the corresponding linear system .  J, stored in `top',
c     `blocks', `bot', and `pivot' has already been
c     factored so we can simply call `crslve' to get the solution
c     corresponding to the right hand side, PHI(Y), stored in `PHI'.
c     Since `crslve' would overwrite `PHI', we actually provide
c     only a copy of `PHI' to `crslve'.
c
      call dcopy((Nsub+1)*neqns,PHI,1,work,1)
c
      call crslve(top, leftbc, neqns, blocks,
     +     neqns,2*neqns,Nsub,bot,neqns-leftbc,pivot,work,delta)
c
c     We complete this calculation by computing the value of the
c     natural criterion function at `Y'. In terms of the quantities we
c     have just computed, this is `g = 0.5*|delta|**2'.
c
      g =  dabs(delta(idamax((Nsub+1)*neqns,delta,1)))
c
c     Here we insert a test here to guard against overflow. If the
c     norm of the Newton correction is large enough that its square
c     overflows, the iteration is in difficulty anyway.
c     We signal difficulty by setting g to -1.0. The negative g
c     condition will be trapped by the calling routine, and the
c     iteration will be terminated.
c
      if (g .GT. overflow_guard) then
                  g = -1.0d0
                  if (PROFILE .GT. PRINT_LEVEL_1) then
                      write(6,*)'Natural criterion function overflow'
                  end if
      else
                  g = 0.5d0 * g**2
                  if (PROFILE .GT. PRINT_LEVEL_2) then
                        write(6,*)'Natural criterion function value',g
                  end if
      endif
c
      return
      end

        SUBROUTINE CRSLVE(TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,
     *             NCLBLK,NBLOKS,BOTBLK,NRWBOT,PIVOT,B,X)
C
C***************************************************************
C
C  C R S L V E  SOLVES THE LINEAR SYSTEM
C                       A*X = B
C  USING THE DECOMPOSITION ALREADY GENERATED IN  C R D C M P.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY  ...
C
C               TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP)
C                         OUTPUT FROM  C R D C M P
C
C               NOVRLP - INTEGER
C                         THE NUMBER OF COLUMNS IN WHICH SUCC-
C                         ESSIVE BLOCKS OVERLAP, WHERE
C                                NOVRLP = NRWTOP + NRWBOT
C
C               NRWTOP - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK TOPBLK
C
C               ARRAY  - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS)
C                         OUTPUT FROM  C R D C M P
C
C               NRWBLK - INTEGER
C                         NUMBER OF ROWS IN K-TH BLOCK
C
C               NCLBLK - INTEGER
C                         NUMBER OF COLUMNS IN K-TH BLOCK
C
C               NBLOKS - INTEGER
C                         NUMBER OF NRWBLK BY NCLBLK BLOCKS IN
C                         THE MATRIX A
C
C               BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP)
C                         OUTPUT FROM  C R D C M P
C
C               NRWBOT - INTEGER
C                         NUMBER OF ROWS IN THE BLOCK BOTBLK
C
C                PIVOT - INTEGER(N)
C                         THE PIVOT VECTOR FROM  C R D C M P
C
C                    B - DOUBLE PRECISION(N)
C                         THE RIGHT HAND SIDE VECTOR
C
C                    X - DOUBLE PRECISION(N)
C                         WORK SPACE
C
C       *** ON RETURN  ...
C
C
C	             X - DOUBLE PRECISION(N)
C                         THE SOLUTION VECTOR
C
C***************************************************************
C
        DOUBLE PRECISION TOPBLK,ARRAY,BOTBLK,X,B
        DOUBLE PRECISION DOTPRD,XJ,XINCRJ,BINCRJ,SWAP
        INTEGER PIVOT(1)
        DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),
     *          BOTBLK(NRWBOT,1),B(1),X(1)
C
C***************************************************************
C
C          ****  DEFINE THE CONSTANTS USED THROUGHOUT  ****
C
C***************************************************************
C
        NRWTP1 = NRWTOP+1
        NRWBK1 = NRWBLK+1
        NVRLP1 = NOVRLP+1
        NRWTP0 = NRWTOP-1
        NRWBT1 = NRWBOT+1
        NROWEL = NRWBLK-NRWTOP
        NRWEL1 = NROWEL+1
        NVRLP0 = NOVRLP-1
        NBLKS1 = NBLOKS+1
        NBKTOP = NRWBLK+NRWTOP
C
C***************************************************************
C
C               ****  FORWARD RECURSION  ****
C
C***************************************************************
C
C          ***  FIRST, IN TOPBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD SOLUTION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 130 J = 1,NRWTOP
           X(J) = B(J)/TOPBLK(J,J)
           IF(J.EQ.NRWTOP)GO TO 120
              XJ = -X(J)
              LOOP = J+1
              DO 110 I = LOOP,NRWTOP
                 B(I) = B(I)+TOPBLK(I,J)*XJ
110           CONTINUE
120        CONTINUE
130     CONTINUE
C
C       ********************************************************
C
C          ***  IN EACH BLOCK ARRAY(,,K)....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        INCR = 0
        DO 280 K = 1,NBLOKS
           INCRTP = INCR+NRWTOP
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD MODIFICATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 220 J = 1,NRWTOP
              INCRJ = INCR+J
              XINCRJ = -X(INCRJ)
              DO 210 I = 1,NRWBLK
                 INCRI = INCRTP+I
                 B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
210           CONTINUE
220        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 240 J = NRWTP1,NRWBLK
              INCRJ = INCR+J
              JPIVOT = PIVOT(INCRJ)
              IF(JPIVOT.EQ.INCRJ)GO TO 225
                 SWAP = B(INCRJ)
                 B(INCRJ) = B(JPIVOT)
                 B(JPIVOT) = SWAP
225           CONTINUE
              BINCRJ = -B(INCRJ)
              LOOP = J-NRWTP0
              DO 230 I = LOOP,NRWBLK
                 INCRI = INCRTP+I
                 B(INCRI) = B(INCRI)+ARRAY(I,J,K)*BINCRJ
230           CONTINUE
240        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD SOLUTION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 270 J = NRWBK1,NBKTOP
              INCRJ = INCR+J
              JRWTOP = J -NRWTOP
              X(INCRJ) = B(INCRJ)/ARRAY(JRWTOP,J,K)
              IF(J.EQ.NBKTOP)GO TO 260
                 XINCRJ = -X(INCRJ)
                 LOOP = J-NRWTP0
                 DO 250 I = LOOP,NRWBLK
                    INCRI = INCRTP+I
                    B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
250              CONTINUE
260           CONTINUE
270        CONTINUE
           INCR = INCR+NRWBLK
280     CONTINUE
C
C       ********************************************************
C
C          ***  FINALLY, IN BOTBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD MODIFICATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        INCRTP = INCR+NRWTOP
        DO 320 J = 1,NRWTOP
           INCRJ = INCR+J
           XINCRJ = -X(INCRJ)
           DO 310 I = 1,NRWBOT
              INCRI = INCRTP+I
              B(INCRI) = B(INCRI)+BOTBLK(I,J)*XINCRJ
310        CONTINUE
320     CONTINUE
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               FORWARD ELIMINATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        IF(NRWBOT.EQ.1)GO TO 350
           DO 340 J = NRWTP1,NVRLP0
              INCRJ = INCR+J
              JPIVOT = PIVOT(INCRJ)
              IF(JPIVOT.EQ.INCRJ)GO TO 325
                 SWAP = B(INCRJ)
                 B(INCRJ) = B(JPIVOT)
                 B(JPIVOT) = SWAP
325           CONTINUE
              BINCRJ = -B(INCRJ)
              LOOP = J-NRWTP0
              DO 330 I = LOOP,NRWBOT
                 INCRI = INCRTP+I
                 B(INCRI) = B(INCRI)+BOTBLK(I,J)*BINCRJ
330           CONTINUE
340        CONTINUE
350     CONTINUE
C
C***************************************************************
C
C               ****  BACKWARD RECURSION  ****
C
C***************************************************************
C
C          ***  FIRST IN BOTBLK....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD SOLUTION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 430 LL = 1,NRWBOT
           J = NVRLP1-LL
           INCRJ = INCR+J
           NRWBTL = NRWBT1-LL
           X(INCRJ) = B(INCRJ)/BOTBLK(NRWBTL,J)
           IF(LL.EQ.NRWBOT)GO TO 420
              XINCRJ = -X(INCRJ)
              LOOP = NRWBOT-LL
              DO 410 I = 1,LOOP
                 INCRI = INCRTP+I
                 B(INCRI) = B(INCRI)+BOTBLK(I,J)*XINCRJ
410           CONTINUE
420        CONTINUE
430     CONTINUE
C
C       ********************************************************
C
C          ***  THEN IN EACH BLOCK ARRAY(,,K)....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 490 L = 1,NBLOKS
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD ELIMINATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           K = NBLKS1-L
           INCR = INCR-NRWBLK
           DO 450 L1 = NRWEL1,NRWBLK
              I = NRWBLK+NRWEL1-L1
              IPLUSN = I+NRWTOP
              LOOP = IPLUSN+1
              INCRN = INCR+IPLUSN
              DOTPRD = X(INCRN)
              DO 440 J = LOOP,NCLBLK
                 INCRJ = INCR+J
                 DOTPRD = DOTPRD-ARRAY(I,J,K)*X(INCRJ)
440           CONTINUE
              X(INCRN) = DOTPRD
              IPVTN = PIVOT(INCRN)
              IF(INCRN.EQ.IPVTN)GO TO 445
                 SWAP = X(INCRN)
                 X(INCRN) = X(IPVTN)
                 X(IPVTN) = SWAP
445           CONTINUE
450        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD MODIFICATION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           INCRTP = INCR+NRWTOP
           DO 460 J = NRWBK1,NCLBLK
              INCRJ = INCR+J
              XINCRJ = -X(INCRJ)
              DO 455 I = 1,NROWEL
                 INCRI = INCRTP+I
                 B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
455           CONTINUE
460        CONTINUE
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD SOLUTION
C
C          CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
           DO 480 LL = 1,NROWEL
              J = NRWBK1-LL
              INCRJ = INCR+J
              NRWELL = NRWEL1-LL
              X(INCRJ) = B(INCRJ)/ARRAY(NRWELL,J,K)
              IF(LL.EQ.NROWEL)GO TO 470
                 XINCRJ = -X(INCRJ)
                 LOOP = NROWEL-LL
                 DO 465 I = 1,LOOP
                    INCRI = INCRTP+I
                    B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ
465              CONTINUE
470           CONTINUE
480        CONTINUE
490     CONTINUE
C
C       ********************************************************
C
C          ***  IN TOPBLK FINISH WITH....
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               BACKWARD ELIMINATION
C
C       CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
        DO 520 L = 1,NRWTOP
           I = NRWTP1-L
           LOOP = I+1
           DOTPRD = X(I)
           DO 510 J = LOOP,NOVRLP
              DOTPRD = DOTPRD-TOPBLK(I,J)*X(J)
510        CONTINUE
           X(I) = DOTPRD
           IPVTI = PIVOT(I)
           IF(I.EQ.IPVTI)GO TO 515
                 SWAP = X(I)
                 X(I) = X(IPVTI)
                 X(IPVTI) = SWAP
515        CONTINUE
520     CONTINUE
        RETURN
        END
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine damp_factor(neqns,leftbc,Nsub,mesh,Y,delta_0,
     +           g_0,top,bot,blocks,pivot,lambda_min,lambda,
     +           PHI,delta,g,Fixed_Jacobian,info,
     +           k_discrete,Y_0,work,fsub,gsub)
c
c***overview***
c
c           This routine takes a discrete solution approximation, `Y', provided
c     on a mesh of `Nsub' subintervals (contained in `mesh') and uses
c     an iteration to produce a suitable value for the damping
c     factor `lambda', for use in a damped Newton update of `Y'. The number
c     of differential equations is `neqns', while the number of boundary
c     conditions at the lefthand endpoint of the problem interval is `leftbc'.
c
c           With the iterate `Y' stored in `Y_0' and the corresponding Newton
c     correction input in `delta_0', this routine predicts a new value
c     for `lambda' and then computes a trial iterate value
c     `Y = Y_0 - lambda*delta_0'. Its acceptability is measured by
c     whether or not it produces a sufficient reduction in the natural
c     criterion function `g = 0.5*|delta|^2', where `delta = inv(J)*PHI',
c     where J is evaluated at `Y_0', the previous iterate, and PHI is the
c     value of the residual function at the trial iterate `Y'. The arrays
c     `top', `bot', `blocks', contain the derivative information associated
c     with the Jacobian matrix J and `pivot' contains the pivoting information
c     for the linear system solution. The new `g' value is compared with the
c     previous value stored in `g_0'. If the reduction is not sufficient, a new
c     damping factor is chosen, if possible, and then this evaluation
c     step is repeated.  If it is found that it is impossible to find a
c     suitable damping factor, this is interpreted as a signal that the
c     entire Newton iteration will not converge. This is indicated by
c     setting `info' = 3. `lambda_min' is a minimum value for the
c     damping factor. A successful return is indicated by `info'= 0. A large
c     natural criterion function value, `g', signifies difficulty in the
c     Newton iteration and is indicated by `info' = 3.
c
c           At the end of a successful search for a damping factor, it may be
c     appropriate to take one or more fixed_Jacobian steps. The criterion is
c     that the current damped Newton step was a full Newton step i.e.
c     `lambda' = 1. This is signaled by setting `Fixed_Jacobian' to TRUE.

c           A by-product of the call to `resid' is the calculation of
c     intermediate quantities associated with the underlying Runge-Kutta method
c     called stages. These are returned in the array called `k_discrete'.
c     This array is exported through the parameter list for later use
c     during the construction and evaluation of the interpolant associated
c     with the Runge-Kutta method. The array, `work', is a work array.
c     The subroutines, `fsub', and `gsub', are provided by the user to define
c     the boundary value ODE and the boundary conditions.
c------------------------------------------------------------------------------
c***declaration of constants***
      integer     Mxs
      parameter   (Mxs=10)
c
c     `Mxs' Maximum number of stages of the RK method.
c
c***declaration of parameters***
c   imports:
      integer            neqns, leftbc, Nsub
      double precision   mesh(0:Nsub), Y(neqns*(Nsub+1))
      double precision   delta_0(neqns*(Nsub+1)), g_0
      double precision   top(neqns*neqns), bot(neqns*neqns)
      double precision   blocks(2*neqns**2*Nsub)
      integer            pivot(neqns*(Nsub+1))
      double precision   lambda_min, lambda
c
c     `neqns' Number of differential equations.
c     `leftbc' Number of boundary conditions at the left
c              end of the problem interval.
c     `Nsub' Number of subintervals into which the problem
c            interval is partitioned.
c     `mesh' Set of points which partition the problem interval.
c     `Y' On input, the current discrete approximation to the solution.
c     `delta_0' Full Newton correction, equal to `inv(J)*PHI',
c               where J and PHI are evaluated at the input iterate, `Y'.
c     `g_0' Value of the natural criterion function at the input iterate.
c     `top' Storage of derivative info for left boundary conditions.
c     `bot' Storage of derivative information for right boundary conditions.
c     `blocks' Storage of derivative information for residual function.
c     `pivot' Pivoting information during linear system solution.
c             Together `top', `blocks', and `bot' contain
c             the almost block diagonal matrix representing the
c             coefficient matrix of the Newton system, in factored form.
c     `lambda_min' Minimum value for the damping factor. If the
c                  algorithm prescribes a value for lambda below this
c                  minimum, it is an indication that the Jacobian is
c                  "effectively singular" and that the Newton iteration
c                  will not converge.
c     `lambda' On input, initial guess for the value of the damping
c              factor used to perform the damped Newton update.
c
c   exports:
c     double precision   Y(neqns*(Nsub+1)), lambda
      double precision   PHI(neqns*(Nsub+1)), delta(neqns*(Nsub+1)), g
      logical            Fixed_Jacobian
      integer            info
      double precision   k_discrete(Mxs*neqns*Nsub)
c
c     `Y' On output, the updated discrete approximation to the
c         solution, evaluated at each meshpoint.
c     `lambda' On output, the updated damping factor used to perform
c              the damped Newton update.
c     `PHI' Value of the residual function associated with the updated iterate.
c     `delta' Newton correction vector, equal to inv(J)*PHI,
c             where PHI is evaluated at the new iterate and J is
c             evaluated at the previous iterate.
c     `g' Value of natural criterion function at the updated iterate.
c     `Fixed_Jacobian' .TRUE. if the next step to be taken should be
c                      a fixed Jacobian step. .FALSE. if the next step
c                      to be taken should be a damped Newton step.
c     `info' Communication flag:
c                        0 - successful iteration converged to within user
c                            specified tolerance.
c                        3 - unsuccessful termination - it was impossible to
c                            obtain a suitable damping factor for the Newton
c                            update (indicative of an "effectively singular"
c                            Jacobian) or evaluation of natural criterion
c                            function overflowed
c     `k_discrete' Discrete stages for all subintervals. The i-th set of
c                  `s*neqns' locations of `k_discrete' contains the `s' stages,
c                  each of length `neqns' corresponding to the i-th subinterval.
c                  Computed by the `resid' routine within the call to
c                  `criterion'.
c   work space:
      double precision   Y_0(neqns*(Nsub+1))
      double precision     work(neqns*(Nsub+1))

c     `Y_0' Current iterate value is saved here during the search
c           for a suitable damping factor and new iterate value.
c     `work' Passed to the `criterion' routine.
c
c***user-supplied subroutines***
      external fsub,gsub
c
c     `fsub' Defines f(t,y) for first order system of differential equations.
c     `gsub' Defines the boundary condition equations.
c
c***declaration of local variables***
      logical            accept
      double precision   sigma, tau
c
c     `accept' Used to monitor the search for a new damping factor.
c     `sigma, tau'  Parameters used to the control selection
c                   of the damping factor - see descriptions below.
c
c***declaration of variables for /IOCNTRL/ common block***
c   imports:
      integer   print_level_0,print_level_1,print_level_2
      parameter (print_level_0 = 0,print_level_1 = 1,
     +              print_level_2 = 2)
      integer   profile
      common /IOCNTRL/ profile
c
c     `print_level_0' No output.
c     `print_level_1' Intermediate output.
c     `print_level_2' Full output.
c     `profile' Controls output, to standard output, of profiling infor-
c               mation such as Newton iteration counts, mesh selection,
c               relative defect estimate. This variable is identical to
c               the `output_control' parameter.
c------------------------------------------------------------------------------
c     Called by: `damped_newt'
c     Calls to : `criterion', `dcopy', `daxpy'
c------------------------------------------------------------------------------
c***initialization***
c     Two parameters which control the search for a suitable damping
c     factor are defined here:
c
c     (a)  `sigma' ensures that the reduction in the size of the natural
c           criterion function `g', evaluated at the new iterate, will be
c           nonnegligible. This keeps the iteration from stalling.
c
c     (b)  `tau' controls how much the values for `lambda' are allowed to
c           change from one step to the next in the iteration performed by
c           this routine.
      sigma = 0.01d0
      tau = 0.1d0
c
c     During the calculation of trial iterates for the determination
c     of a suitable damping factor, `Y' will be overwritten, so
c     save the current iterate, stored in `Y', in the vector, `Y_0'.
      call dcopy((Nsub+1)*neqns,Y,1,Y_0,1)
c
c***iterative selection of new damping factor***
c
c     Use a repeat-until loop to control the search for the new
c     damping factor. The loop will be controlled by the flags,
c     `accept' and `info': `accept' = .TRUE. or info .NE. 0
c     will terminate the loop. (Implement the loop in FORTRAN
c     using an if-goto pair.)
c
      accept = .FALSE.
c     REPEAT
 2                                                         continue
c
c     ***compute trial iterate using current damping factor***
c
c            Compute an updated iterate value, stored in `Y', using the
c            current value for the damping factor, stored in `lambda'
c            `Y = Y_0 - lambda*delta_0'.
c
             call dcopy((Nsub+1)*neqns,Y_0,1,Y,1)
             call daxpy((Nsub+1)*neqns,-lambda,delta_0,1,Y,1)
c
c     ***evaluate natural criterion function at this new iterate***
c
c            Measure the suitability of the current value of the
c            damping factor by seeing how much the trial iterate
c            reduces the value of the natural criterion function. Compute
c            the value of the criterion function at `Y' by calling the
c            `criterion' routine. The value is returned in `g'; the
c            residual function, evaluated at `Y' is returned in `PHI'.
c            The value `inv(J)*PHI', (with `J' evaluated at the previous
c            iterate), is returned in `delta'.
c
             call criterion(neqns,leftbc,Nsub,mesh,Y,top,
     +         blocks, bot, pivot,PHI,delta,g,k_discrete,work,fsub,gsub)
c
c            If the natural criterion function value, `g', has overflowed
c            the `criterion' routine will trap that condition and set
c            `g = -1.0' ( a value which cannot otherwise arise).
c            Since the natural criterion function is related to the
c            size of the Newton correction, a large `g' value signifies
c            difficulty in the Newton iteration. Set `info'=3 and exit.
c
             if (g .LT. 0.0d0) then
                 info = 3
             else
c               (g >= 0.0d0)
c                Natural criterion function has not overflowed.
c
c            ***accept current damping factor or compute a new one,
c                                    depending on criterion function value***
c
c                Check value of `g' to ensure that the reduction in size is
c                sufficient. If so `lambda' and thus `Y' are acceptable.
c                If not, select a new `lambda', if possible.
c
                 if (g .LE. (1.0d0 - 2.0d0*lambda*sigma)*g_0) then
                   accept = .TRUE.
c
c                  This damped Newton step is successful. The new iterate
c                  value has already been stored in `Y' so no update is
c                  needed at this point.
c
c                  Decide whether or not to take a fixed Jacobian step
c                  next time.
c
                   if (lambda .EQ. 1.0d0) Fixed_Jacobian = .TRUE.
                 else
c                  `(g > (1.0d0-2.0d0*lambda*sigma)*g_0')
c                  Compute a new damping factor.
                   lambda = max(tau*lambda,
     +                     (lambda**2 * g_0)/((2*lambda-1.0d0)*g_0+g))
c
                   if (profile .GT. print_level_2) then
                     write(6,*)'Damping factor lambda=', lambda
                   endif
c
c                  Check this new `lambda' value. If it is too small, then
c                  the Newton iteration will not converge, set `info' = 3
c                  to signal this.
c
                   if (lambda .LT. lambda_min) info = 3
                 endif
c                End of `if (g.LE.(1.0d0-2.0d0*lambda*sigma)*g_0)'
             end if
c            End of `if (g.LT.0.0d0)'
c
c     UNTIL `(accept .OR. (info .GT. 0))'
c
                                  if ((.NOT. accept) .AND.
     +                              (info .LE. 0)) go to 2
c
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine damped_newt(neqns,leftbc,Nsub,mesh,Y,newtol,
     +    lambda,PHI,top,bot,blocks,pivot,Fixed_Jacobian,
     +    convrg,delta,delta_0_norm,g,info,k_discrete,delta_0,
     +    work,fsub,gsub,dfsub,dgsub)
c
c***overview***
c
c            This routine takes a mesh of `Nsub' subintervals, partitioning
c     the problem interval and a discrete initial approximation provided
c     at the meshpoints, (contained in the vectors `Y' and `mesh',
c     respectively) and  uses a full or damped Newton step to refine a
c     solution approximation for the associated discrete nonlinear system.
c     Convergence occurs if the Newton correction is less than the
c     tolerance `newtol' and is indicated by `convrg' = TRUE. The
c     number of differential equations is `neqns', while the number of
c     boundary conditions at the lefthand endpoint of the problem interval
c     is `leftbc'. A successful return is indicated by `info'= 0. If this
c     routine encounters a singular matrix, it signals this be setting
c     `info' = 2. A large natural criterion function value, `g', signifies
c     difficulty and is indicated by `info' = 3. On the other hand, if the
c     iteration appears to be converging well, this routine will signal that
c     the iteration can switch over to a fixed Jacobian step (see below).
c
c            At the beginning of each damped Newton step, it is assumed that
c     the current iterate value, `Y', and the value of the residual function,
c     `PHI', evaluated at `Y' are available. If the previous step was also
c     a damped Newton step, this routine also expects `delta = inv(J)*PHI',
c     where J, the Jacobian, is evaluated at the previous iterate and
c     `delta_0_norm', the norm of the previous value of `delta_0 = inv(J)*PHI',
c     where both J and `PHI' are associated with the previous iterate. The
c     step begins with the evaluation of then Jacobian of the residual function
c     at the current iterate, `Y'. The Newton system, consisting of the
c     Jacobian as the coefficient matrix and `PHI', as the right hand side, is
c     then solved to obtain the full Newton correction, which we store in the
c     vector `delta_0'. The factored Jacobian is available in the vectors
c     `top', `bot', `blocks'. Pivoting information for the solution of the
c     linear system is stored in `pivot'. During the call to `damp_factor' the
c     current iterate is copied into `Y_0', a damping factor, `lambda' is
c     computed, and then a damped Newton step is taken to produce a new
c     iterate, `Y = Y_0 - lambda*delta_0'.
c
c           At the end of a successful damped Newton step, it may be
c     appropriate to take one or more fixed Jacobian steps . The criterion
c     is that the current damped Newton step was a full Newton step
c     (`lambda'=1). This is signaled by setting `Fixed_Jacobian' to TRUE.
c     The fixed Jacobian step, for reasons associated with efficient data
c     management, expects to receive a number of quantities : (i) the new
c     iterate, `Y', (ii) the residual function evaluated at this new iterate,
c     `PHI', (iii) `delta = inv(J)*PHI', where J, the Jacobian, has been
c     evaluated rather at the previous iterate  and (iv) the value of the
c     natural criterion function, `g = 0.5*|delta|**2'. All these values are
c     to be naturally computed during the current damped Newton step and
c     thus will be already set up for the fixed Jacobian step if appropriate.
c
c           A by-product of the evaluation of the right hand side of the
c     Newton system is the calculation of intermediate quantities associated
c     with the underlying Runge-Kutta method called stages. These are
c     returned to the `damped_newt' routine when it calls the `damp_factor'
c     routine and are stored in the array called `k_discrete'.
c     This array is exported through the parameter list for later use
c     during the construction and evaluation of the interpolant associated
c     with the Runge-Kutta method. The array, `work', is a work array provided
c     to the `damp_factor' and `newmat' routines. The subroutines, `fsub',
c     `gsub', `dfsub', and `dgsub', are provided by the user to define the
c     boundary value ODE and the boundary conditions and their derivatives.
c------------------------------------------------------------------------------
c***declaration of constants***
      integer    Mxs
      parameter  (Mxs=10)
c
c     `Mxs'    is the maximum number of stages of the RK method.
c
c     We wish to be able to detect when the criterion function value
c     will overflow. Thus the overflow_guard to 9.22d18.
c
      double precision       overflow_guard
      parameter   (overflow_guard = 9.22d18)
c
c***declaration of parameters***
c   imports:
      integer            neqns, leftbc, Nsub
      double precision   mesh(0:Nsub), Y(neqns*(Nsub+1))
      double precision   newtol, lambda, PHI(neqns*(Nsub+1))
      double precision   delta(neqns*(Nsub+1)), delta_0_norm
c
c     `neqns' Number of differential equations.
c     `leftbc' Number of boundary conditions at the left
c              end of the problem interval.
c     `Nsub' Number of subintervals into which the problem
c            interval is partitioned.
c     `mesh' Set of points which partition the problem interval.
c     `Y' On input, the current discrete approximation to the solution.
c     `newtol' Value to be applied to the Newton iteration,
c              i.e. the iteration stops if the Newton correction is
c              less than or equal to the user specified tolerance.
c              A blended relative/absolute tolerance is applied.
c              If `delta_j' is the Newton correction for, `Y_j', the jth
c              component of `Y', then we require
c                    `(|delta_j| / |Y_j|) <= newtol'.
c     `lambda' On input, the damping factor used in the previous damped
c              Newton step. If its value is set to 0 on input, the
c              previous step was not a damped Newton step.
c     `PHI' On input, the value of the residual function on each
c           subinterval evaluated at the current iterate `Y'.
c     `delta' On input, it is `inv(J)*PHI' where J is evaluated at the
c             previous iterate and `PHI' is evaluated at the input
c             iterate, if the previous step was a damped Newton step.
c     `delta_0_norm' On input, is |inv(J)*PHI| where J and `PHI' were
c                    evaluated at the previous iterate if the
c                    previous step was a damped Newton step.
c                    It is undefined if the previous step was not a
c                    damped Newton step.
c   exports:
c     double precision   Y(neqns*(Nsub+1)), lambda
c     double precision   PHI(neqns*(Nsub+1))
      double precision   top(neqns*neqns)
      double precision   bot(neqns*neqns)
      double precision   blocks(2*neqns**2*Nsub)
      integer            pivot(neqns*(Nsub+1))
      logical            Fixed_Jacobian, convrg
c     double precision   delta(neqns*(Nsub+1))
c     double precision   delta_0_norm
      double precision   g
      integer            info
      double precision   k_discrete(Mxs*neqns*Nsub)
c
c     `Y' On output, the updated discrete approximation to the solution.
c     `lambda' On output, the value for the new damping factor.
c     `PHI' On output, the value of the residual function on each
c           subinterval evaluated at the updated iterate `Y'.
c     `top' Storage of derivative information for left boundary
c           conditions.
c     `bot' Storage of derivative information for right boundary
c           conditions.
c     `blocks' Storage of derivative information for residual function.
c     `pivot' Pivoting information during linear system solution.
c             Together `top', `blocks', and `bot' contain
c             the almost block diagonal matrix representing the
c             coefficient matrix of the Newton system, in factored form.
c     `Fixed_Jacobian' .TRUE. if the next step to be taken should be
c                      a fixed Jacobian step. .FALSE. if the next step
c                      to be taken should be a damped Newton step.
c     `convrg' Used to check for convergence of the Newton iteration.
c     `delta' Newton correction vector, equal to `inv(J)*PHI(Y)',
c             where J is evaluated at the input iterate and `PHI'
c             and `PHI' is evaluated at the new iterate `Y'.
c     `delta_0_norm' on output, |inv(J)*PHI| where J and `PHI' are
c                    evaluated at the input iterate.
c     `g' Value of natural criterion function at the updated iterate.
c     `info' Internal communication flag:
c            0 - successful iteration converged to within user
c                specified tolerance.
c            2 - unsuccessful termination - a singular coefficient
c                matrix was encountered during the attempted solution
c                of the Newton system.
c            3 - unsuccessful termination - it was impossible to obtain
c                a suitable damping factor for the Newton update
c                (indicative of an "effectively singular" Jacobian) or
c                evaluation of natural criterion function overflowed.
c     `k_discrete' Vector containing the discrete stages for all
c                  subintervals. The i-th set of `s*neqns' locations of
c                  `k_discrete' contains the `s' stages, each of length
c                  `neqns' corresponding to the i-th subinterval.
c   work space:
      double precision   delta_0(neqns*(Nsub+1))
      double precision
     +                  work(neqns*(2*Nsub+3+2*(Mxs+1)*neqns))
      integer            i_Y_0, i_work

c     `delta_0' Full Newton correction, equal to `inv(J)*PHI',
c               where J and `PHI' are evaluated at the input iterate.
c     `work' Work array provided to the `damp_factor' and
c            `newmat' routines.
c     `i_Y_0' Index to the start of the part of `work' to be passed
c             to `damp_factor' for use as the work array `Y_0'.
c     `i_work' Index to the start of the remainder of `work' to be passed
c              to `damp_factor' for use as the work array `work'.
c
c***user-supplied subroutines***
      external fsub,gsub,dfsub,dgsub
c
c     `fsub' Defines f(t,y) for the first order
c            system of differential equations, y' = f(t,y).
c     `gsub' Defines the boundary condition equations.
c     `dfsub' Defines the Jacobian of the system
c             of differential equations.
c     `dgsub' Defines the Jacobian of the boundary conditions.
c
c***declaration of local variables***
      integer            factor
      double precision   g_0, lambda_min,mu
      integer            j
c
c     `factor' Used by the `colrow' routine to indicate trouble; a
c              return with `factor' equal to -1 indicates that a singular
c              Jacobian was encountered.
c     `g_0' Value of the natural criterion function at the input iterate.
c     `lambda_min' Parameter used to the control selection of the
c                  damping factor - see descriptions below.
c     `mu' Predicted value for the new damping factor.
c     `j' Loop index, 1,...,`(Nsub+1)*neqns'.
c
c***declaration of variables for /IOCNTRL/ common block***
c   imports:
      integer   print_level_0,print_level_1,print_level_2
      parameter (print_level_0 = 0,print_level_1 = 1,
     +               print_level_2 = 2)
      integer   profile
      common /IOCNTRL/ profile
c
c     `print_level_0' No output.
c     `print_level_1' Intermediate output.
c     `print_level_2' Full output.
c     `profile' Controls output, to standard output, of profiling infor-
c               mation such as Newton iteration counts, mesh selection,
c               relative defect estimate. This variable is identical to
c               the `output_control' parameter.
c------------------------------------------------------------------------------
c     Called by : `NewIter'
c     Calls to  : `NewMat', `colrow', `damp_factor', `daxpy', `idamax'
               integer   idamax
c------------------------------------------------------------------------------
c***initialization of work array indexes***
      i_Y_0 = 1
      i_work = neqns*(Nsub+1)+1

c***perform a damped newton step
c     The flag `info' is used to monitor for trouble; It is initialized to 0.
      info = 0
c
c     `lambda_min' is a parameter which controls the search for a
c     suitable damping factor; it gives the minimum value
c     for the damping factor. It is used here in the initial guess for the
c     damping factor. It is also used with the `damp_factor' routine.
      lambda_min = 0.01d0
c
c     First construct the Newton system at the input iterate.
c     It is assumed that the value of  the residual function
c     at the current iterate was obtained at the end of the previous
c     iteration step and is available in the vector `PHI'.  The
c     coefficient  matrix  is  computed  by calling the `NewMat'
c     routine and is returned in the  vectors  `top', `blocks', `bot'.
c
      call NewMat(neqns,leftbc,Nsub,mesh,Y,top,blocks,bot,
     +                               k_discrete,work,dfsub,dgsub)
c
c     The coefficient matrix constructed by `NewMat' has an
c     "almost block diagonal" structure. `colrow', a package of
c     subroutines designed to solve such systems is now employed.
c     This routine will return with `factor' = -1 if  the coefficient
c     matrix of the Newton system is singular and `factor'=0 otherwise.
c     The solution of the linear system, which gives the Newton correction,
c     is returned in the vector `delta_0'.
c
      call colrow((Nsub+1)*neqns,top,leftbc,neqns,blocks,
     +              neqns,2*neqns,Nsub,bot,neqns-leftbc,
     +              pivot,PHI,delta_0,factor)
c
      if ( factor .EQ. -1 ) then
            info = 2
c           This is a signal that a singular Jacobian matrix was detected.
      end if
c
c
      if (info .EQ. 0) then
c
c          `colrow' was successful in solving the linear system to obtain
c          the Newton correction. Attempt to compute the corresponding
c          new iterate using a damped Newton step. At this time, the value
c          of the natural criterion function for the current iterate is
c          also computed and stored in `g_0'. (It can be shown
c          that, for this case, `g_0 = 0.5 * |delta_0|**2'.)
c
c          Check here to see if the natural criterion function
c          value will overflow. If so, trap the condition. Since the
c          natural criterion function is related to the size of the
c          Newton correction, a large `g' value signifies difficulty in
c          the Newton iteration. Set `info' = 3 and exit.
c
           g_0 = dabs(delta_0(idamax((Nsub+1)*neqns, delta_0,1)))
c
           if (g_0 .GT. overflow_guard) then
                 if (profile .GT. print_level_1) then
                     write(6,*)'Natural criterion function overflow'
                     info = 3
                 end if
           else
                 if (profile .GT. print_level_1) then
                     write(6,*)'Norm of Newton correction',g_0
                 endif
                 g_0 = 0.5d0 * g_0**2
           endif
c
c          Test for convergence now. If the norm of the scaled Newton
c          correction, `delta_0', is sufficiently small, then the iteration
c          has converged.
c
           convrg = .TRUE.
           do 20 j = 1, (Nsub+1)*neqns
             if( dabs(delta_0(j)) .GT. newtol*(dabs(Y(j))+1.0d0))
     +                convrg = .FALSE.
  20       continue
c
           if (convrg) then
c              Update iterate with full Newton correction.
c              The current iterate is still available in `Y' and
c              the Newton correction is in `delta_0'.
               call daxpy((Nsub+1)*neqns,-1.0d0,delta_0,1,Y,1)
c
c              The stages of the RK method are updated in the next call to
c              the `interp_setup' routine with this new iterate value.
c
           else
c              Attempt to update the Newton iterate with a damped Newton
c              correction. Compute a prediction for the new damping factor.
c              The value for `lambda' must be between `lambda_min' and 1.0d0.
c              If the input `lambda' value is 0.0d0, this is a signal that the
c              previous iteration was not a damped Newton iteration and
c              appropriate information for the prediction of a new `lambda'
c              value is not available, in which case set the predicted
c              value to 1.0d0.
c
               if (lambda .EQ. 0.0d0 ) then
                   lambda = 1.0d0
               else
c                  For temporary storage overwrite `delta = delta - delta_0'.
                   call daxpy((Nsub+1)*neqns,-1.0d0,delta_0,1,delta,1)
c
c                  Use the norm of this difference in the prediction for
c                  `lambda'.
                   mu = (lambda*delta_0_norm)/
     +                  dabs(delta(idamax((Nsub+1)*neqns,delta,1)))
                   lambda = max(lambda_min, min(1.0d0, mu))
               endif
c
c              Print out current value for `lambda', if appropriate.
               if (profile .GT. print_level_1) then
                 write(6,*)'Initial damping factor lambda =', lambda
               endif
c
c              The previous value of the norm of the Newton correction,
c              stored in `delta_0_norm', was just used in the prediction
c              of the new `lambda' value, and is no longer needed.
c              Update the value of `delta_0_norm' using the current value
c              of `delta_0' computed by `colrow' above (the norm of `delta_0'
c              was already computed during the calculation of `g_0').
c
               delta_0_norm = dsqrt(2.0d0*g_0)
c
c              Compute a new iterate by updating the current iterate with
c              the (possibly damped) Newton correction.
c
               call damp_factor(neqns,leftbc,Nsub,mesh,Y,delta_0,
     +              g_0,top,bot,blocks,pivot,lambda_min,lambda,
     +              PHI,delta,g,Fixed_Jacobian,info,k_discrete,
     +              work(i_Y_0),work(i_work),fsub,gsub)

               if ((info .EQ. 0) .AND.
     +               (profile .GT. print_level_2)) then
                     write(6,*)'Norm of damped Newton correction',
     +                                         lambda*delta_0_norm
               endif
c
c              If successful, this routine will compute an appropriate
c              damping factor, returned in `lambda', the updated iterate,
c              returned in `Y', the corresponding residual function value,
c              returned in `PHI', `delta = inv(J)*PHI', and `g', the natural
c              criterion function, evaluated at the new iterate.
c
           endif
c          End of `if (convrg)'
c
      endif
c     End of `if (info.EQ.0)'
c
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine defect_estimate(neqns,Nsub,mesh,Y,defect,
     +    defect_norm,info,z,z_prime,def_1,def_2,f_1,f_2,
     +    temp_1,temp_2,k_discrete,k_interp,work,fsub)
c
c***overview***
c
c        This routine uses the discrete solution approximation, `Y',
c     obtained on the mesh stored in `mesh', plus the stages of the underlying
c     RK method, stored in `k_discrete', plus some new stages, stored in
c     `k_interp', to construct an interpolant for `Y'. `neqns' is the number of
c     differential equations; `Nsub' is the the number of subintervals into
c     which the mesh divides the problem interval. This interpolant is then used
c     to compute an estimate of the defect on each subinterval, for each
c     solution component. The i-th group of `neqns' locations of `defect' is
c     used to store the relative defect information associated with the i-th
c     subinterval. The max norm of the relative defect is returned in
c     `defect_norm'. A successful return is indicated by `info = 0'. If the
c     `defect_norm' is bigger than the constant parameter `defect_threshold'
c     then `info' is set to 4.
c
c         The primary information required to define the interpolant on the i-th
c     subinterval is the current solution approximation at the (i-1)-st and
c     i-th mesh points (which can be obtained from `Y'), and some
c     stages associated with the i-th interval.  The first `s' stages come
c     from the discrete formula, and have already been computed. They are
c     available in the `s' vector locations making up the i-th set of
c     `s*neqns' locations of the vector,'k_discrete'. The remaining
c     `s_star - s' stages are required only for the interpolant and must
c     now be computed . For each subinterval, this is done through a call
c     to the `interp_setup' routine, which returns the new stage evaluations
c     for that subinterval. The resultant stages are then stored in the
c     appropriate locations of the `k_interp' vector (which has the same
c     format as the `k_discrete' vector).
c
c           On the i-th subinterval the evaluation of the interpolant at
c     `t = mesh(i-1)+tau*hi' is given by `z(t) = yim1 + hi sum br(tau)*kr',
c     where `hi' is the subinterval size, 'yim1' is the solution approximation
c     at `mesh(i-1)', `kr' is the r-th stage, `br(tau)' is the corresponding
c     weight polynomial evaluated at `tau' and the sum runs from 1 to `s_star'.
c
c            Once the interpolant  is available, we can sample
c     the defect at the point, `mesh(i-1)+tau*hi', by substituting the
c     interpolant into the differential equation and subtracting the right
c     side from the left side, i.e. the defect on the i-th subinterval is
c     `delta = z'(mesh(i-1)+tau*hi)-f(mesh(i-1)+tau*hi,z(mesh(i-1)+tau*hi))'.
c     A blended relative/absolute defect measure, `defect', equal to
c     `|delta|/(|f|+1)' is computed, where `f' is the corresponding right hand
c     side of the ODE.
c
c            The sample point is available from the /sample_point/ common
c     block and the subsequent evaluation of the weights is done by a call
c     to the `interp_weights' routine. To improve the reliability of the
c     defect estimate, we also sample at `1-tau' and choose the larger of
c     the two defect values for the estimate, `delta'.
c
c            `z','z_prime,'def_1','def_2','f_1','f_2','temp_1','temp_2'
c     are work arrays used for temporary storage here while `work'
c     is a work array passed to `interp_setup'. The subroutine, `fsub',
c     is provided by the user to define the boundary value ODEs.
c------------------------------------------------------------------------------
c***declaration of constants***
      integer Mxs
      parameter(Mxs=10)
c
c     `Mxs'    Maximum value for the total number of stages of the interpolant.
c
      double precision defect_threshold
      parameter (defect_threshold = 1.0d-1)
c
c     `defect_threshold' Used to monitor the size of the defect.
c
c***declaration of parameters***
c   imports:
      integer neqns, Nsub
      double precision mesh(0:Nsub), Y(neqns*(Nsub+1))
      double precision k_discrete(Mxs*neqns*Nsub)
c
c     `neqns' Number of differential equations.
c     `Nsub' Number of subintervals.
c     `mesh' Array of points defining the subintervals.
c     `Y' Contains the discrete solution at the mesh points
c         the i-th set of `neqns' locations of `Y' contains
c         the solution approximation at the i-th mesh point.
c     `k_discrete' Contains the `s' discrete stage values for each subinterval.
c
c   exports:
      double precision         defect(Nsub*neqns), defect_norm
      integer                  info
      double precision k_interp(Mxs*neqns*Nsub)
c
c     `defect' Contains the corresponding blended relative/absolute
c              defect measure.
c     `defect_norm' Max norm of the defect measure.
c     `info' Communication flag :
c                        0 - Successful return.
c                        4 - Size of the defect norm was too large and the
c                            solution from the current discrete system can not
c                            be trusted.
c     `k_interp' Contains the `s_star-s' extra stage values for each
c                subinterval need to form the interpolant.
c
c   work space:
      double precision z(neqns),z_prime(neqns)
      double precision def_1(neqns),def_2(neqns)
      double precision f_1(neqns),f_2(neqns)
      double precision temp_1(neqns),temp_2(neqns)
      double precision  work(neqns)

c     `z' Used to accumulate argument for fsub evaluation
c     `z_prime' Used to save derivative of interpolant
c     `def_1,def_2' Contain the values of the defect at the sample
c                   points `tau' and `1-tau' for each subinterval.
c     `f_1,f_2' Contain the corresponding fsub value at the
c               the two sample points.
c     `temp_1', `temp_2' Temporary location for storage of scaled defect
c                        at `tau' and `1-tau'.
c     `work' Work array passed to `interp_setup'.
c
c***user-supplied subroutines***
      external fsub
c
c     `fsub' Defines f(t,y) for the first order system of differential
c            equations, y' = f(t,y).
c
c***declaration of local variables***
      integer i,j
      double precision hi, weights_1(Mxs), weights_2(Mxs)
      double precision weights_1_prime(Mxs), weights_2_prime(Mxs)
      double precision estimate_1, estimate_2
c
c     `i' Loop index from 1 to Nsub.
c     `j' Loop index from 1 to `neqns'.
c     `hi' Size of the i-th subinterval.
c     `f' Used to store `fsub' evaluation.
c     `weights_1', `weights_2' Contains the weight polynomials of the
c                              interpolant evaluated at the sample point `tau'
c                              and `1-tau'.
c     `weights_1_prime', `weights_2_prime' Contains the first derivatives of
c                                          the weight polynomials of the
c                                          interpolant evaluated at the sample
c                                          points `tau' and `1-tau'.
c     `estimate_1', `estimate_2' Maximum scaled defect at `tau' and `1-tau'
c
c***declaration of variables from /RK_s/ common block***
      integer s
      common /RK_s/ s
c     `s' Number of discrete stages.
c
c***declaration of variables from /interp_s_star/ common block***
      integer s_star
      common /interp_s_star/ s_star
c     `s_star' Total number of stages required to form the
c              interpolant on each subinterval. It includes all the
c              stages of the discrete formula plus the additional
c              stages required for the interpolant.
c
c***declaration of variables from /sample_point/ common block***
      double precision tau_star
      common /sample_point/ tau_star
c
c     `tau_star' Relative position of one of the sample points within each
c                subinterval. The other sample point is `1-tau_star'.
c-----------------------------------------------------------------------------
c     Called by: `mirkdc'
c     Calls to : `interp_weights', `interp_setup', `sum_stages', `fsub'
c                `daxpy', `dcopy'
      integer    idamax
c------------------------------------------------------------------------------
c     Setup the weights, evaluated at the first sample point.
c
      call interp_weights(s_star,tau_star,weights_1,weights_1_prime)
c
c     Setup the weights, evaluated at the second sample point.
c
      call interp_weights(s_star,1.0d0-tau_star,weights_2,
     +                                        weights_2_prime)
c
      do 30 i = 1, Nsub
            hi = mesh(i) - mesh(i-1)
c
c           We must first setup the new stages for this subinterval.
c
            call interp_setup(neqns,mesh(i-1),hi,Y((i-1)*neqns+1),
     +            Y(i*neqns+1),s,k_discrete((i-1)*s*neqns+1),
     +            s_star,k_interp((i-1)*(s_star-s)*neqns+1),work,fsub)
c
c           The new stages for the i-th subinterval have been computed
c           and stored in the ith set of `(s_star-s)*neqns' locations of
c           `k_interp'.
c
c     ***Sample Point 1***
c
c            With all the information needed to form the interpolant on the
c            i-th subinterval already available, we can compute the value of
c            the interpolant, and its first derivative at the sample point
c            by calling the `sum_stages' routine.
c
             call sum_stages(neqns,hi,Y((i-1)*neqns+1),
     +                        s,k_discrete((i-1)*s*neqns+1),s_star,
     +                        k_interp((i-1)*(s_star-s)*neqns+1),
     +                        weights_1,weights_1_prime,z,z_prime)
c
c            This yields the value of the interpolant, in `z', and its first
c            derivative, in `z_prime'.   We must next evaluate the righthand
c            side of the system of ode's at `z', i.e. compute f(z) and store
c            in `f_1'.
c
             call fsub(neqns,mesh(i-1)+tau_star*hi,z,f_1)
c
c            Since the calculation of `z_prime' was already complete after
c            the return from `sum_stages', we can now compute the defect which
c            we store in `def_1'.
c
             call daxpy(neqns,-1.0d0,f_1,1,z_prime,1)
             call dcopy(neqns,z_prime,1,def_1,1)
c
c            Compute the norm of the scaled defect for Sample Point 1.
c
             do 10 j =1,neqns
                     temp_1(j) =
     +                        def_1(j)/(dabs(f_1(j))+1.0d0)
  10         continue
c
             estimate_1=dabs(temp_1(idamax(neqns,temp_1,1)))
c
c     ***Sample Point 2***
c
c            With all the information needed to form the interpolant on the
c            i-th subinterval already available, we can compute the value of
c            the interpolant, and its first derivative at the sample point
c            by calling the `sum_stages' routine.
c
             call sum_stages(neqns,hi,Y((i-1)*neqns+1),
     +                       s,k_discrete((i-1)*s*neqns+1),s_star,
     +                       k_interp((i-1)*(s_star-s)*neqns+1),
     +                       weights_2,weights_2_prime,z,z_prime)
c
c            This yields the value of the interpolant, in `z', and its first
c            derivative, in `z_prime'.  We must evaluate the righthand side
c            of the system of ode's at `z', i.e. compute f(z).
c
             call fsub(neqns,mesh(i-1)+(1.0d0-tau_star)*hi,z,f_2)
c
c            Since the calculation of `z_prime' was already complete after
c            the return from `sum_stages', we can now compute the defect,
c            which we store in `def_2'.
c
             call daxpy(neqns,-1.0d0,f_2,1,z_prime,1)
             call dcopy(neqns,z_prime,1,def_2,1)
c
             do 20 j =1,neqns
                  temp_2(j) =
     +                     def_2(j)/(dabs(f_2(j))+1.0d0)
 20          continue
c
             estimate_2=dabs(temp_2(idamax(neqns,temp_2,1)))
c
c     ***Compare defect estimates from the two sample points***
c
             if (estimate_1 .GT. estimate_2) then
                call dcopy(neqns,temp_1,1,defect((i-1)*neqns+1),1)
             else
                call dcopy(neqns,temp_2,1,defect((i-1)*neqns+1),1)
             endif
c
c            The scaled defect for the i-th subinterval has now been copied
c            into the ith group of `neqns' locations of the vector `defect'.
c
 30   continue
c
c     Compute the max norm of the defect estimate.
      defect_norm = dabs(defect(idamax(Nsub*neqns,defect,1)))
c
c     If this defect_norm is bigger than `defect_threshold' then
c     this is an indication that, the defect is no longer a
c     small relative perturbation of the original system of ODEs
c     and the solution from the current discrete system should
c     not be trusted. It is probably better to start over as if the
c     Newton iteration had failed. Here, we check the size of `defect_norm',
c     and set `info' = 4 (to cause a mesh doubling), if it is too large.
      if (defect_norm .GT. defect_threshold) info = 4
c
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine fixed_jacob(neqns,leftbc,Nsub,mesh,Y,newtol,
     +       delta,g,PHI,top,bot,blocks,pivot,Fixed_Jacobian,lambda,
     +       convrg,info,k_discrete,Y_hat,PHI_hat,work,fsub,gsub)
c
c***overview***
c
c           This routine takes a mesh of `Nsub' subintervals, partitioning
c     the problem interval and a discrete initial approximation provided
c     at the meshpoints, (contained in the vectors `Y' and `mesh',
c     respectively) and uses a Fixed Jacobian step, as part of a Newton
c     iteration, to produce an updated solution approximation. Convergence
c     occurs if the Newton correction is less than the tolerance `newtol'
c     and is indicated by `convrg' = TRUE. The number of differential equations
c     is `neqns', while the number of boundary conditions at the lefthand
c     endpoint of the problem interval is `leftbc'. A successful return is
c     indicated by `info'= 0. A large natural criterion function value,'g',
c     signifies difficulty in the Newton iteration and is indicated by
c     `info' = 3.
c
c          When a fixed Jacobian step is to be taken, the following
c     quantities are assumed to be available from the previous iteration
c     step  (i) the iterate, `Y', (ii) the residual function evaluated at this
c     iterate, `PHI', (iii) `delta = inv(J)*PHI', where J, the Jacobian ,
c     available in factored form in the arrays `top', `bot', `blocks', and
c     `pivot', has not been evaluated at `Y', but rather at some previous
c     iterate and (iv) the value of the natural criterion function, `g', at
c     the iterate, `Y'.
c
c          The fixed Jacobian step consists of first computing  a (potential)
c     new iterate, `Y_hat = Y - delta'. The acceptability of `Y_hat' is
c     determined  by evaluating the natural criterion function at `Y_hat',
c     giving `g_hat'. During the calculation of `g_hat, it is necessary to
c     evaluate the residual function at `Y_hat' and store it in the vector,
c     `PHI_hat', and a new value for `delta = inv(J)*PHI_hat', both of which
c     will then be available for the next step of the iteration. The
c     acceptability of the new iterate value is determined by comparing `g_hat'
c     with the value of the natural criterion function for the previous iterate,
c     stored in `g'. If `g_hat < rho*g' then the new iterate value is accepted
c     and another fixed Jacobian step is attempted. If `g_hat > rho*g', the
c     code returns to taking damped Newton steps (the value of `rho' is set in
c     this routine). This is signaled by setting `Fixed_Jacobian' to FALSE.
c     The value of `lambda', the damping factor, is set to 0 as well. This is
c     to indicate appropriate information for the prediction of a new damping
c     factor is not available for the upcoming damped Newton step.
c
c           A by-product of the evaluation of the right hand side of the
c     Newton system is the calculation of intermediate quantities associated
c     with the underlying Runge-Kutta method called stages. These are
c     returned in the array called `k_discrete', which is exported through
c     the parameter list for later use during the construction and evaluation
c     of the interpolant associated with the Runge-Kutta method. The array,
c     `work', is a work array. The subroutines, `fsub', and `gsub', are provided
c     by the user to define the boundary value ODE and the boundary conditions.
c------------------------------------------------------------------------------
c***declaration of constants***
      integer  Mxs
      parameter (Mxs=10)
c
c     `Mxs'    is the maximum number of stages of the RK method.
c
c***declaration of parameters***
c   imports:
      integer            neqns, leftbc, Nsub
      double precision   mesh(0:Nsub), Y(neqns*(Nsub+1))
      double precision   newtol, delta(neqns*(Nsub+1)), g
      double precision   PHI(neqns*(Nsub+1)), top(neqns*neqns)
      double precision   bot(neqns*neqns), blocks(2*neqns**2*Nsub)
      integer            pivot(neqns*(Nsub+1))
c
c     `neqns'   the number of differential equations
c     `leftbc'  the number of boundary conditions at the left
c               end of the problem interval
c     `Nsub'    the number of subintervals into which the problem
c               interval is partitioned
c     `mesh'    the set of points which partition the problem interval.
c     `Y' on input, the current discrete approximation to the solution.
c     `newtol' the tolerance value to be applied to the Newton iteration,
c              i.e. the iteration stops if the Newton correction
c              is less than or equal to the user specified
c              tolerance. A blended relative/absolute tolerance
c              is applied. If `delta_j' is the Newton correction
c              for the jth component of Y, Y_j, then we require
c                          |delta_j| / (|Y_j|+1) <= newtol.
c     `delta'   on input, the Newton correction vector, equal to
c               inv(J)*PHI(Y), where J is evaluated at a previous iterate.
c     `g'       on input, the value of natural criterion function at the
c               current iterate, `Y'.
c     `PHI'    on input, is the value of the residual function on each
c              subinterval, evaluated at `Y'.
c     `top'     storage of derivative information for left boundary conditions.
c     `bot'     storage of derivative information for right boundary conditions.
c     `blocks'  storage of derivative information for residual function.
c     `pivot'   pivoting information associated with the factorization
c               of J. Together `top', `blocks', and `bot' contain
c               the almost block diagonal matrix representing the
c               coefficient matrix of the Newton system, in factored form.
c   exports:
c     double precision   Y(neqns*(Nsub+1)),delta(neqns*(Nsub+1))
c     double precision   g, PHI(neqns*(Nsub+1))
      logical            Fixed_Jacobian
      double precision   lambda
      logical            convrg
      integer            info
      double precision   k_discrete(Mxs*neqns*Nsub)
c
c     `Y' On output, the updated discrete approximation to the solution.
c     `delta' On output, equal to inv(J)*PHI where the
c             Jacobian J is evaluated at the previous iterate
c             and PHI is evaluated at the updated iterate `Y_hat'
c     `g' On output, value of the criterion function evaluated
c         at the updated iterate.
c     `PHI' On output, is value of the residual function
c           on each subinterval evaluated at the updated iterate.
c     `Fixed_Jacobian' .TRUE. if the next step to be taken should be a
c                      fixed Jacobian step. .FALSE. if the next step
c                      to be taken should be a damped Newton step.
c     `lambda' Set to 0 to indicate appropriate information for the
c              prediction of a new lambda is not available during the
c              upcoming damped Newton step
c     `convrg' Used to check for convergence of the Newton iteration.
c     `info' Communication flag:
c                 0 - successful iteration converged to
c                     within user specified tolerance.
c                 3 - unsuccessful termination - A large natural
c                     criterion function value, g, signifies
c                     difficulty in the Newton iteration.
c     `k_discrete' Discrete stages for all subintervals. The ith set of
c                  `s*neqns' locations of `k_discrete' contains the `s'
c                  stages, each of length `neqns', corresponding
c                  to the ith subinterval. Computed by `resid'
c                  from within `criterion`.
c   work space:
      double precision    Y_hat(neqns*(Nsub+1))
      double precision    PHI_hat(neqns*(Nsub+1))
      double precision    work(neqns*(Nsub+1))

c     `Y_hat'  Trial iterate value.
c     `PHI_hat' Residual function value associated with the trial iterate.
c     `work' Work space passed to `criterion'.
c
c***user - supplied subroutines***
      external fsub,gsub
c     `fsub' Defines f(t,y) for the first order
c            system of differential equations, y' = f(t,y).
c     `gsub' Defines the boundary condition equations.
c
c***declaration of local variables***
      double precision   g_hat, rho
      integer            j
c
c     `g_hat' Value of the natural criterion function at the trial iterate.
c     `rho' Used in measuring performance of fixed Jacobian iteration.
c     `j' Loop index, 1,...,(Nsub+1)*neqns.
c
c***declaration of variables for /IOCNTRL/ common block***
c   imports:
      integer   print_level_0,print_level_1,print_level_2
      parameter (print_level_0 = 0,print_level_1 = 1,
     +              print_level_2 = 2)
      integer   profile
      common /IOCNTRL/ profile
c
c     `print_level_0' No output.
c     `print_level_1' Intermediate output.
c     `print_level_2' Full output.
c     `profile' Controls output, to standard output, of profiling infor-
c               mation such as Newton iteration counts, mesh selection,
c               relative defect estimate. This variable is identical to
c               the `output_control' parameter.
c------------------------------------------------------------------------------
c     Called by : `Newiter'
c     Calls to  : `criterion', `dcopy', `daxpy', `idamax'
      integer       idamax
c------------------------------------------------------------------------------
c     The flag `info' is used to monitor for trouble; it is initialized to 0.
      info = 0
c
c     When a fixed Jacobian step has been taken, the decision about
c     whether or not to follow with another fixed Jacobian step or not
c     will be determined by the size of the decrease in the value of the
c     natural criterion function. This decrease must be at least a factor
c     of `rho'. Here `rho' is initialized.
      rho = 0.5d0
c
c     In addition to the new iterate `Y' being available, we also
c     have the corresponding evaluation of the residual function,
c     `PHI', the Newton correction, `delta', for the current fixed
c     Jacobian step, and the evaluation of the criterion function
c     for the this iterate, `g = 0.5*|inv(J)*PHI|**2 = 0.5*|delta|**2'.
c
c     If the Newton correction, `delta', is less than the tolerance, the
c     iteration has converged. A blended relative/absolute tolerance is
c     employed.
c
      if (profile .GT. print_level_1) then
           write(6,*)'Norm of Newton correction',
     +       dabs(delta(idamax((Nsub+1)*neqns,delta,1)))
      endif
c
      convrg = .TRUE.
      do 10 j = 1, (Nsub+1)*neqns
        if(dabs(delta(j)).GT.newtol*(dabs(Y(j))+1.0d0)) convrg=.FALSE.
  10  continue
c
      if (convrg) then
c          Update current iterate before exiting.
           call daxpy((Nsub+1)*neqns,-1.0d0,delta,1,Y,1)
c
c          The stages of the RK method are updated in the next
c          call to the `interp_setup' routine with this new
c          iterate value.
c
      else
c          If the iteration has not converged, some care must be taken
c          in using the fixed Newton step `delta' to  update the
c          current iterate. The correction will not be useful if it
c          does not lead to an iterate, `Y', which provides a sufficient
c          reduction in the size of the natural criterion function.
c          The acceptability of the new iterate is determined by
c          examining the resultant reduction in the size of the natural
c          criterion function. The following calls compute the trial iterate
c          `Y_hat = Y - delta'.
c
           call dcopy((Nsub+1)*neqns,Y,1,Y_hat,1)
           call daxpy((Nsub+1)*neqns,-1.0d0,delta,1,Y_hat,1)
c
c          The natural criterion function is now evaluated to provide a
c          measure of how good this trial iterate is. By-products
c          of this calculation, are the quantities `PHI_hat', the evaluation
c          of the residual function at `Y_hat', and `delta = 0.5*
c          |inv(J)*PHI_hat|' for the NEXT iterate. These MAY be useful for the
c          next iteration step. If `Y_hat' turns out to be not acceptable `PHI',
c          the residual function at `Y', will be needed, so here it is not
c          overwritten. Also, `g', the value of the natural criterion function
c          at `Y', will be needed for comparison with `g_hat', the value of
c          the natural criterion function at `Y_hat', so `g' is not overwritten
c          either. The previous value of `delta = inv(J)*PHI' is no longer
c          needed in either case, so it is overwritten with the new value
c          `inv(J)*PHI_hat'.
c
           call criterion(neqns,leftbc,Nsub,mesh,Y_hat,top,blocks,bot,
     +          pivot,PHI_hat,delta,g_hat,k_discrete,work,fsub,gsub)
c
c          If the natural criterion function value is about to overflow
c          the `criterion' routine will trap that condition and set
c          g_hat = -1.0 ( a value which cannot otherwise arise).
c          Since the natural criterion function is related to the
c          size of the Newton correction, a large value signifies
c          difficulty in the Newton iteration. Set `info'=3 and exit
c          the loop.
c
           if (g_hat .LT. 0.0d0) then
               info = 3
           else
c             The Natural criterion function value has not overflowed.
c             If the new iterate results in a sufficient reduction in
c             the size of the natural criterion function, accept it.
c
              if ( g_hat .LE. rho * g ) then
c                 The new iterate is acceptable.
c                 Update the appropriate values `Y','PHI','g' in
c                 preparation for another fixed Jacobian step.
c                 (The flag `Fixed_Jacobian' remains TRUE.)
c
                  call dcopy((Nsub+1)*neqns,Y_hat,1,Y,1)
                  call dcopy((Nsub+1)*neqns,PHI_hat,1,PHI,1)
                  g = g_hat
              else
c                 The new iterate is not acceptable. Set the flag so that
c                 next iteration step is a damped Newton step. This also
c                 requires setting `lambda=0'.
                  Fixed_Jacobian = .FALSE.
                  lambda = 0.0d0
c
c                 If `g_hat' is greater than `rho' * `g' but less than `g'
c                 itself, it may be reasonable to use the new iterate
c                 and corresponding residual function value to begin the
c                 subsequent full Newton step. These new values, `Y_hat'
c                 and `PHI_hat', are copied to `Y' and `PHI'.
c
                  if (g_hat .LT. g) then
                      call dcopy((Nsub+1)*neqns,Y_hat,1,Y,1)
                      call dcopy((Nsub+1)*neqns,PHI_hat,1,PHI,1)
                  else
c                     Begin the full Newton step using the values in
c                     `Y' and `PHI' from the beginning of this step.
c                     A call to `resid' is made here to obtain the `k_discrete'
c                     values, evaluated at `Y', as these must be up-to-date
c                     with the `Y' and `PHI' values before the damped Newton
c                     step begins.
c
                      call resid(neqns,leftbc,Nsub,mesh,Y,PHI,
     +                                 k_discrete,work,fsub,gsub)
                  endif
c                 `if (g_hat .LT. g)'
              endif
c             `if (g_hat .LE. rho*g)'
           endif
c          `if (g_hat .LT. 0.0d0)'
      endif
c     `if (convrg)'
c
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine half_mesh(Nsub,mesh_current,mesh_new)
c
c***overview***
c
c     This routine takes as input a mesh of `Nsub' subintervals which
c     partitions the problem interval, stored in the vector `mesh_current',
c     The routine returns a new mesh, stored in `mesh_new', which is exactly
c     twice as fine as the given one, i.e. has `2*Nsub' subintervals, obtained
c     by splitting each subinterval of the given mesh in half.
c------------------------------------------------------------------------------
c***declaration of parameters***
c   imports:
      integer                Nsub
      double precision       mesh_current(0:Nsub)
c
c     `Nsub' Number of subintervals in the input mesh.
c     `mesh_current' Input meshpoints.
c
c   exports:
      double precision       mesh_new(0:2*Nsub)
c
c     `mesh_new' Meshpoints making up the new mesh.
c
c***declaration of local variables***
      integer       i
c
c     `i' loop index from 1 to `Nsub'.
c
c------------------------------------------------------------------------------
c     Called by: `mirkdc', `mesh_selector'.
c------------------------------------------------------------------------------
c
      mesh_new(0) = mesh_current(0)
c
      do 10 i = 1, Nsub
c
           mesh_new(2*i) = mesh_current(i)
c
c          Generate a new mesh point at location 2*i-1 using the average of
c          the adjacent meshpoint values.
c
           mesh_new(2*i-1) = (mesh_current(i) + mesh_current(i-1))/2.0d0
c
 10   continue
c
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine interp_eval(neqns,Nsub,mesh,Y,t,z,z_prime,
     +                                     k_discrete,k_interp)
c***overview***
c
c          This routine uses the current discrete solution approximation, `Y',
c     obtained on a mesh of `Nsub' subintervals, stored in `mesh', plus
c     additional stage information associated with the underlying discrete RK
c     method, to compute an evaluation of an interpolant to the discrete
c     solution.  This interpolant is evaluated at `t', with the value returned
c     in `z', and the value of the first derivative returned in `z_prime'.
c
c          The interpolant is computed as follows: Given `t', assumed to
c     be somewhere in the problem interval, the subinterval of `mesh' that
c     contains `t' is determined first, i.e. we find i such that
c                `mesh(i-1) =< t < mesh(i) for t < mesh(Nsub)' or
c                `i = Nsub for t = mesh(Nsub)'.
c     The fraction of the way `t' is through the i-th subinterval, a quantity
c     which we call `tau', is also determined during this step. Then the
c     interpolant on this i-th subinterval is evaluated - the interpolant has
c     the form `z(mesh(i-1) + tau*hi) = yim1 + hi sum br(tau)*kr'.
c     The weights, `br(tau)', can be computed with a call to `interp_weights'.
c
c          The remaining information required to define the interpolant
c     on the i-th subinterval is the current solution approximation at the
c     (i-1)st mesh point, `yim1', which can be obtained from `Y', and the
c     stages associated with the i-th subinterval. The first `s' stages,
c     contained in `k_discrete', come from the discrete formula and
c     remaining `s_star - s' stages, contained in `k_interp', are necessary
c     for the interpolant. Each stage is a vector of length `neqns' where
c     `neqns' is the number of differential equations.
c
c          `k_discrete' is the vector containing the discrete stages for
c     all subintervals. The i-th set of `s*neqns' locations of `k_discrete'
c     contains the `s' stages, each of length `neqns' corresponding to the
c     i-th subinterval. `k_interp' plays a similar role for the `s_star-s'
c     stages associated with the interpolant.
c------------------------------------------------------------------------------
c***declaration of constants***
      integer Mxs
      parameter(Mxs=10)
c
c     `Mxs' Maximum value for the total number of stages of the interpolant.
c
c***declaration of parameters***
c   imports:
      integer neqns, Nsub
      double precision mesh(0:Nsub), Y(neqns*(Nsub+1)), t
      double precision k_discrete(Mxs*neqns*Nsub)
      double precision k_interp(Mxs*neqns*Nsub)
c
c     `neqns' Number of differential equations.
c     `Nsub' Number of subintervals.
c     `mesh' Array of points defining the subintervals.
c     `Y' Contains the discrete solution at the mesh points;
c         the i-th set of `neqns' locations of `Y' contains
c         the solution approximation at the i-th mesh point.
c     `t' Point at which the interpolant is to be evaluated.
c     `k_discrete' Contains the `s' discrete stage values for each subinterval.
c     `k_interp' Contains the `s_star-s' extra stage values for each
c                subinterval need to form the interpolant.
c   exports:
      double precision z(neqns), z_prime(neqns)
c
c     `z' Value of the interpolant at `t'.
c     `z_prime' Value of the first derivative of the interpolant at `t'.
c
c***declaration of local variables***
      integer i
      double precision hi
      double precision tau, weights(Mxs), weights_prime(Mxs)
c
c     `i' Index for subinterval containing `t', i ranges from 1 to `Nsub'.
c     `hi' Size of the i-th subinterval.
c     `tau' Equal to (t-mesh(i-1))/hi.
c     `weights' Contains the weight polynomials of the interpolant
c               evaluated at `tau'.
c     `weights_prime' Contains the first derivatives of the weight
c                     polynomials of the interpolant evaluated at `tau'.
c
c***declaration of variables from /RK_s/ common block***
      integer s
      common /RK_s/ s
c     `s' Number of discrete stages.
c
c***declaration of variables from /interp_s_star/ common block***
      integer s_star
      common /interp_s_star/ s_star
c     `s_star' Total number of stages required to form the
c              interpolant on each subinterval. It includes all the
c              stages of the discrete formula plus the additional
c              stages required for the interpolant.
c-----------------------------------------------------------------------------
c     Called by: `mirkdc'.
c     Calls to: `interval', `interp_weights', `sum_stages'.
c------------------------------------------------------------------------------
c     Call `interval' to determine the subinterval containing t.
c     The index of the subinterval is returned in `i'.
c
      call interval(Nsub,mesh,t,i)
c
c     Compute `tau' for use in the evaluation of the weight polynomials.
c
      hi = mesh(i) - mesh(i-1)
      tau = (t - mesh(i-1))/hi
c
c     Setup the weights, evaluated at `tau'.
c
      call interp_weights(s_star, tau, weights,weights_prime)
c
c     We need to evaluate the interpolant and its derivative at `tau'.
c     This is done by taking weighted sums of the stages, using the weights
c     we have just computed, and the stages corresponding to the i-th
c     subinterval, and calling the `sum_stages' routine.
c
      call sum_stages(neqns,hi,Y((i-1)*neqns+1),s,
     +                    k_discrete((i-1)*s*neqns+1),s_star,
     +                    k_interp((i-1)*(s_star-s)*neqns+1),
     +                    weights, weights_prime, z, z_prime)
c
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine interp_setup(neqns,tim1,hi,yim1,yi,s,ki_discrete,
     +                                     s_star,ki_interp,y,fsub)
c***overview***
c
c         The purpose of this routine is to perform some preliminary
c     calculations to set up stages for later use in the evaluation of the
c     interpolant. The interpolant is defined on the i-th subinterval by,
c          `Z_i(tim1 + hi*tau) = yim1 + hi sum b_r(tau)*k_r',
c     where the stages, `k_r', depend on `yim1' and `yi', the solution
c     approximations at the left and right endpoints, `hi', the size of the
c     subinterval, 'tim1' the left endpoint, and the coefficients of the
c     underlying Runge-Kutta method. The weight polynomials, `b_r(tau)', can be
c     evaluated through a call to the `interp_weights' routine.
c     The first `s' stages of each subinterval were already obtained during the
c     computation of the discrete solution, so only the last `s_star - s' stages
c     need to be computed now.
c
c         This routine will construct the `s_star-s' extra stages needed by
c     the interpolant for a single subinterval. This routine takes as input
c     the vector, `ki_discrete', containing the `s' stages associated with
c     the discrete formula on the i-th subinterval, and then uses the
c     Runge-Kutta coefficients for the extra stages, available from the
c     `/interp_coeffs/' common block, to construct the `s_star' new stages.
c     These are returned in the vector `ki_interp'. Each stage is a vector of
c     length `neqns', where `neqns' is the number of differential equations.
c     `y' is a work array. The subroutine, `fsub', is provided by the user to
c     define the boundary value ODE.
c-------------------------------------------------------------------------------
c***declaration of constants***
      integer  Mxs
      parameter (Mxs=10)
c
c     `Mxs' Maximum value for the total number of stages of the interpolant.
c
c***declaration of parameters***
c   imports:
      integer            neqns
      double precision   tim1, hi, yim1(neqns), yi(neqns)
      integer            s, s_star
      double precision   ki_discrete(s*neqns)
c
c     `neqns' Number of ordinary differential equations.
c     `tim1' Left endpoint of the ith subinterval.
c     `hi' Size of the current subinterval.
c     `yim1' Solution approximation at the left endpoint of the current
c                                                             subinterval.
c     `yi' Solution approximation at the right endpoint of the current
c                                                            subinterval.
c     `s' Number of discrete stages.
c     `s_star' Total number of stages required to form the
c              interpolant on each subinterval. It includes all the
c              stages of the discrete formula plus the additional
c              stages required for the interpolant.
c     `ki_discrete' Vector of `s*neqns' components. The j-th
c                   component is of length `neqns' and corresponds
c                   to the j-th of the `s' stage values associated
c                   with the discrete formula on the i-th subinterval.
c   exports:
      double precision ki_interp((s_star-s)*neqns)
c
c     `ki_interp' Vector with same structure as `ki_discrete'.
c                 The j-th vector component corresponds to the j-th
c                 of the new stages needed to form the interpolant.
c
c   work space:
      double precision y(neqns)
c
c     `y' Used to hold the argument for the function evaluation
c         for each new stage as it is being computed.
c
c***user-supplied subroutines***
      external fsub
c
c     `fsub' Defines f(t,y) for the first order system
c            of differential equations, y' = f(t,y).
c
c***declaration of local variables***
      integer          r, j

c     `r,j' Loop indexes from 1 to `s_star-s'.
c
c***declaration of variables from common block /interp_coeff/***
c   imports:
      double precision c_star(Mxs), v_star(Mxs), x_star(Mxs*Mxs)
      common /interp_coeff/ c_star,v_star,x_star
c
c     `c_star' First component contains the abscissa for the first new stage.
c              The second component contains the abscissa for the second new
c              stage, and so on.
c
c     `v_star' First component contains the blending coefficient for the first
c              new stage. The second component contains the blending coefficient
c              for the second new stage, and so on.
c
c     `x_star' Contains the coupling coefficients for each of the new stages.
c              The Runge-Kutta coefficient matrix, X_star, is related to the
c              data structure, `x_star' as follows (the matrix is stored
c              columnwise in the vector):
c
c                     X_star(i,j) = x_star( (j-1)*s + i )
c
c------------------------------------------------------------------------------
c    Called by: `defect_estimate'.
c    Calls to : `fsub', `dcopy', `daxpy'.
c------------------------------------------------------------------------------
c***construction of the new stages***
c
      do 20 r = 1 , s_star-s
c
c          initialize `y' to zero
           call dcopy(neqns,0.0d0,0,y,1)
c
c          Accumulate the argument for the rth new stage in `y'.
c
c          The contributions from the discrete stages.
           do 5 j = 1 , s
c
c                The contribution from the j-th discrete stage.
                 call daxpy(neqns,x_star((j-1)*(s_star-s)+r),
     +                          ki_discrete((j-1)*neqns+1),1,y,1)
 5         continue
c
c          The contributions from the previous new stages.
c
            do 10 j = 1 , r-1
c
c                 The contribution from the j-th new stage.
                  call daxpy(neqns,x_star((j+s-1)*(s_star-s)+r),
     +                           ki_interp((j-1)*neqns+1),1,y,1)
 10         continue
c
c           Multiply the sum of the stages by `hi'.
            call dscal(neqns,hi,y,1)
c
c           Add on the contribution from `yim1'.
            call daxpy(neqns,(1-v_star(r)),yim1,1,y,1)
c
c           Add on the contribution from `yi'.
            call daxpy(neqns,v_star(r),yi,1,y,1)
c
c           All contributions to the argument for r-th new stage have
c           been made. Simply evaluate `fsub' at this argument to get the
c           new stage. The result is copied into the appropriate locations
c           of `ki_interp'.
c
            call fsub(neqns,tim1+c_star(r)*hi,y,
     +                      ki_interp((r-1)*neqns+1))
 20   continue
c
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine interp_tableau (method)
c
c***overview***
c
c        This routine defines all the extra coefficients needed to augment
c      the discrete RK formula with an interpolant. `method' defines the MIRK
c      method to be used. The possible values for `method' and the corresponding
c      MIRK scheme associated with each of these values are listed below.
c
c             Value of `method' |   MIRK formula
c            -------------------------------------
c                  121          |   MIRK121 scheme
c                  221          |   MIRK221 scheme
c                  232          |   MIRK232 scheme (abscissa = 0,2/3)
c                  2321         |   MIRK232 scheme (abscissa = 1,1/3)
c                  343          |   MIRK343 scheme
c                  453          |   MIRK453 scheme
c                  563          |   MIRK563 scheme
c                  5631         |   An improved MIRK563 scheme
c
c      The interpolant is defined on a subinterval by subinterval basis, as
c      follows: When t can be expressed as t_i + hi*tau for some i and tau,
c      with hi = t_i+1 - t_i then the interpolant, Z(t) = Zi(t_i+ hi*tau) =
c      yi + hi sum br(tau)*kr, where the kr depend on information from the ith
c      subinterval as well as certain blending and coupling coefficients that
c      are specific to the MIRK formula and interpolant. br(tau) is the rth
c      weight polynomial of the interpolant evaluated at tau.
c
c        The main purpose of this routine is to define the parameter which
c      gives,`s_star', the total number of stages for the interpolant and the
c      coefficients for the `s_star-s' extra stages needed by the interpolant.
c      Once defined by this routine, the coefficients are then
c      available from the /interp_coeff/ common block, while `s_star' is
c      available from the /interp_s_star/ common block. The `order' of the
c      interpolant is available through the /interp_order/ common block.
c
c         This routine also defines the sample point, relative to the size
c      of a subinterval, at which the defect should be sampled. This is the
c      point at which the interpolant must be evaluated and is thus specific
c      to the particular interpolant being implemented.  The defect sample
c      point, which we call `tau_star', is made available through the
c      /sample_point/ common block.
c
c            The final portion of formula dependent information involves
c      the expressions for the weight polynomials, arising as part of the
c      definition of the interpolant for the discrete formula. This information
c      is provided through the `interp_weights' routine which evaluates all
c      the weight polynomials and their derivatives at a given point. This
c      routine can be called by other parts of the code whenever the weight
c      polynomials and derivatives are required - i.e. whenever an evaluation
c      of the interpolant is required.
c-----------------------------------------------------------------------------
c***declaration of constants***
      integer    Mxs
      parameter  (Mxs=10)
c
c     `Mxs'     is the maximum number of stages of the RK method.
c
c***declaration of parameters***
c   imports:
      integer       method
c
c     `method'   defines the method to be used
c
c***declaration of variables from common block /interp_s_star/***
      integer s_star
      common /interp_s_star/ s_star
c
c     `s_star'  is the total number of stages required to form the
c               interpolant on each subinterval. It includes all the
c               stages of the discrete formula plus the additional
c               stages required for the interpolant.
c
c***declaration of variables from common block /interp_coeff/***
c     exports:
        double precision c_star(Mxs), v_star(Mxs), x_star(Mxs*Mxs)
        common /interp_coeff/ c_star,v_star,x_star
c
c       The first component of `c_star' contains the abscissa
c       for the first new stage. The second component of `c_star' contains
c       the abscissa for the second new stage, and so on.
c
c       The first component of `v_star' contains the blending coefficient
c       for the first new stage. The second component of `v_star' contains
c       the blending coefficient for the second new stage, and so on.
c
c       The matrix, X_star, is related to the data structure, `x_star'
c       as follows (the matrix is stored columnwise in the vector):
c
c                      X_star(i,j) = x_star( (j-1)*s + i )
c
c       `x_star' contains the coupling coefficients for each of the new
c       stages.
c
c***declaration of variables for the /sample_point/ common block***
c     exports:
        double precision tau_star
        common /sample_point/ tau_star
c
c       `tau_star'  is the point within each subinterval at which the
c                       defect is to be sampled.
c
c***declaration of variables for the /rk_method/ common block***
c     exports:
        integer   amethod
        common   /rk_method/ amethod
c
c       `amethod'   the same as `method', defines the method to be
c       used. It is defined twice here so that it may be
c       passed into a common block.
c
c***declaration of variables for the /interp_order/ common block***
c     exports:
        integer p
        common /interp_order/ p
c
c       `p' is the order of the interpolant => local error is O(h**(p+1)).
c------------------------------------------------------------------------------
c      called by: `mirkdc'
c      calls to: `dcopy'
c------------------------------------------------------------------------------
c      Translate short form of method values.
       if (method .EQ. 2) method = 221
       if (method .EQ. 4) method = 343
       if (method .EQ. 6) method = 563

c      Assign the value of `method' to `amethod' so that it may be
c      passed into a common block.
       amethod = method

c      Beginning method definition based on value of `method'.
c
       if (method. EQ. 121) then
c
c        SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c        Here we give coefficients defining the interpolant for a 1-stage,
c        2nd order MIRK, stage order 1, with a 2nd order interpolant.
c
         s_star = 3
c
         c_star(1) = 0.0d0
         c_star(2) = 1.0d0
c
         v_star(1) = 0.0d0
         v_star(2) = 1.0d0
c
c        Column 1 of x_star
c
         x_star(1) = 0.0d0
         x_star(2) = 0.0d0
c
c        Column 2 of x_star
c
         x_star(3) = 0.0d0
         x_star(4) = 0.0d0
c
c        Column 3 of x_star
c
         x_star(5) = 0.0d0
         x_star(6) = 0.0d0
c
c        Define the sample point for the interpolant of the 1-stage, 2nd
c        order MIRK discrete formula.
c
         tau_star = 0.25d0
c
c        Define the value of the order of the interpolant, `p'.
c
         p = 2
c
      elseif (method. EQ. 221) then
c
c        SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c        Here we give coefficients defining the interpolant for a 2-stage,
c        2nd order MIRK, stage order 1, with a 2nd order interpolant.
c
         s_star = 2
c
c        Define the sample point for the interpolant of the 2-stage, 2nd
c        order MIRK discrete formula.
c
         tau_star = 0.25d0
c
c        Define the value of the order of the interpolant, `p'.
c
         p = 2
c
      elseif (method. EQ. 232) then
c
c        SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c        Here we give coefficients defining the interpolant for a 2-stage,
c        3rd order MIRK, stage order 2, with a 3rd order interpolant.
c
         s_star = 3
c
         c_star(1) = 1.0d0
c
         v_star(1) = 1.0d0
c
c        Column 1 of x_star
c
         x_star(1) = 0.0d0
c
c        Column 2 of x_star
c
         x_star(2) = 0.0d0
c
c        Column 3 of x_star
c
         x_star(3) = 0.0d0
c
c        Define the sample point for the interpolant of the 2-stage, 3rd
c        order MIRK discrete formula.
c
         tau_star  = 0.25d0
c
c        Define the value of the order of the interpolant, `p'.
c
         p = 3
c
      elseif (method. EQ. 2321) then
c
c        SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c        Here we give coefficients defining the interpolant for the
c        reflection of the above 2-stage, 3rd order MIRK, stage order 2,
c        with a 3rd order interpolant.
c
         s_star = 3
c
         c_star(1) = 0.0d0
c
         v_star(1) = 0.0d0
c
c        Column 1 of x_star
c
         x_star(1) = 0.0d0
c
c        Column 2 of x_star
c
         x_star(2) = 0.0d0
c
c        Column 3 of x_star
c
         x_star(3) = 0.0d0
c
c        Define the sample point for the interpolant of the 2-stage, 3rd
c        order MIRK discrete formula.
c
         tau_star = 0.25d0
c
c        Define the value of the order of the interpolant, `p'.
c
         p = 3
c
      elseif (method. EQ. 343) then
c
c        SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c        Here we give coefficients defining the interpolant for a 3-stage,
c        4th order MIRK, stage order 3, with a 4th order interpolant.
c
c        There is one extra stage. The corresponding row of the RK tableau
c        is    3/4 | 27/32 |  3/64  -9/64  0   0
c
         s_star = 4
c
         c_star(1) = 3.0d0/4.0d0
c
         v_star(1) = 27.0d0/32.0d0
c
c        Column 1 of x_star
c
         x_star(1) = 3.0d0/64.0d0
c
c        Column 2 of x_star
c
         x_star(2) = -9.0d0/64.0d0
c
c        Column 3 of x_star
c
         x_star(3) = 0.0d0
c
c        Column 4 of x_star
c
         x_star(4) = 0.0d0
c
c        Define the sample point for the 3-stage, 4th order MIRK
c        discrete formula.
c
         tau_star = 0.226d0
c
c        Define the value of the order of the interpolant, `p'.
c
         p = 4
c
      elseif (method .EQ. 453) then
c
c        SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c        Here we give coefficients defining the interpolant for a 4-stage,
c        5th order MIRK, stage order 3, with a 5th order interpolant.
c
         s_star = 6
c
         c_star(1) = 4.0d0/5.0d0
         c_star(2) = 13.0d0/23.0d0
c
         v_star(1) = 4.0d0/5.0d0
         v_star(2) = 13.0d0/23.0d0
c
c        Column 1 of x_star
c
         x_star(1) = 14.0d0/1125.0d0
         x_star(2) = 1.0d0/2.0d0
c
c        Column 2 of x_star
c
         x_star(3) = -74.0d0/875.0d0
         x_star(4) = 4508233.0d0/1958887.0d0
c
c        Column 3 of x_star
c
         x_star(5) = -128.0d0/3375.0d0
         x_star(6) = 48720832.0d0/2518569.0d0
c
c        Column 4 of x_star
c
         x_star(7) = 104.0d0/945.0d0
         x_star(8) = -27646420.0d0/17629983.0d0
c
c        Column 5 of x_star
c
         x_star(9)  = 0.0d0
         x_star(10) = -11517095.0d0/559682.0d0
c
c        Column 6 of x_star
c
         x_star(11) = 0.0d0
         x_star(12) = 0.0d0
c
c        Define the sample point for the interpolant of the 4-stage, 5th
c        order MIRK discrete formula.
c
         tau_star = 0.30d0
c
c        Define the value of the order of the interpolant, `p'.
c
         p = 5
c
      elseif (method .EQ. 4531) then
c
c        SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c        Here we give coefficients defining the interpolant for a 4-stage,
c        5th order MIRK, stage order 3, with a 5th order interpolant.
c
         s_star = 6
         c_star(1) = 7.0d0/13.0d0
         c_star(2) = 13.0d0/23.0d0
c
         v_star(1) = 4.0d0/5.0d0
         v_star(2) = 23.0d0/33.0d0
c
c        Column 1 of x_star
c
         x_star(1) = 36445.0d0/1542294.0d0
         x_star(2) = 1.0d0/2.0d0
c
c        Column 2 of x_star
c
         x_star(3) = -125759.0d0/1999270.0d0
         x_star(4) = 1499810527.0d0/2327157756.0d0
c
c        Column 3 of x_star
c
         x_star(5) = -3224608.0d0/11567205.0d0
         x_star(6) = -9339357584.0d0/2742721641.0d0
c
c        Column 4 of x_star
c
         x_star(7) = 915050.0d0/16194087.0d0
         x_star(8) = -11845932250.0d0/4918765257.0d0
c
c        Column 5 of x_star
c
         x_star(9)  = 0.0d0
         x_star(10) = 514365992257.0d0/113365827828.0d0
c
c        Column 6 of x_star
c
         x_star(11) = 0.0d0
         x_star(12) = 0.0d0
c
c        Define the sample point for the interpolant of the 4-stage, 5th
c        order MIRK discrete formula.
c
         tau_star = 0.30d0
c
c        Define the value of the order of the interpolant, `p'.
c
         p = 5
c
      elseif (method .EQ. 563) then
c
c        SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c        Here we give coefficients defining the interpolant for a 5-stage,
c        6th order MIRK, stage order 3, with a 6th order interpolant.
c
c        There are four extra stages. The corresponding rows of the RK tableau:
c        7/16 | 7/16 |  1547/32768  -1225/32768  749/4096   -287/2048 -861/16384
c        3/8  | 3/8  |  83/1536      -13/384     283/1536   -167/1536  -49/512
c        9/16 | 9/16 |  1225/32768  -1547/32768  287/2048 -749/4096    861/16384
c        1/8  | 1/8  |  233/3456     -19/1152       0         0           0
c                                            ...   -5/72      7/72      -17/216
c
         s_star = 9
c
         c_star(1) = 7.0d0/16.0d0
         c_star(2) = 3.0d0/8.0d0
         c_star(3) = 9.0d0/16.0d0
         c_star(4) = 1.0d0/8.0d0
c
         v_star(1) = 7.0d0/16.0d0
         v_star(2) = 3.0d0/8.0d0
         v_star(3) = 9.0d0/16.0d0
         v_star(4) = 1.0d0/8.0d0
c
c        Column 1 of x_star
c
         x_star(1) = 1547.0d0/32768.0d0
         x_star(2) = 83.0d0/1536.0d0
         x_star(3) = 1225.0d0/32768.0d0
         x_star(4) = 233.0d0/3456.0d0
c
c        Column 2 of x_star
c
         x_star(5) = -1225.0d0/32768.0d0
         x_star(6) = -13.0d0/384.0d0
         x_star(7) = -1547.0d0/32768.0d0
         x_star(8) = -19.0d0/1152.0d0
c
c        Column 3 of x_star
c
         x_star(9)  = 749.0d0/4096.0d0
         x_star(10) = 283.0d0/1536.0d0
         x_star(11) = 287.0d0/2048.0d0
         x_star(12) = 0.0d0
c
c        Column 4 of x_star
c
         x_star(13) = -287.0d0/2048.0d0
         x_star(14) = -167.0d0/1536.0d0
         x_star(15) = -749.0d0/4096.0d0
         x_star(16) = 0.0d0
c
c        Column 5 of x_star
c
         x_star(17) = -861.0d0/16384.0d0
         x_star(18) = -49.0d0/512.0d0
         x_star(19) = 861.0d0/16384.0d0
         x_star(20) = 0.0d0
c
c        Column 6 of x_star
c
         x_star(21) = 0.0d0
         x_star(22) = 0.0d0
         x_star(23) = 0.0d0
         x_star(24) = -5.0d0/72.0d0
c
c        Column 7 of x_star
c
         x_star(25) = 0.0d0
         x_star(26) = 0.0d0
         x_star(27) = 0.0d0
         x_star(28) = 7.0d0/72.0d0
c
c        Column 8 of x_star
c
         x_star(29) = 0.0d0
         x_star(30) = 0.0d0
         x_star(31) = 0.0d0
         x_star(32) = -17.0d0/216.0d0
c
c        Column 9 of x_star
c
         x_star(33) = 0.0d0
         x_star(34) = 0.0d0
         x_star(35) = 0.0d0
         x_star(36) = 0.0d0
c
c        Define the sample point for the interpolant of the 5-stage, 6th
c        order MIRK discrete formula.
c
         tau_star = 0.7156d0
c
c        Define the value of the order of the interpolant, `p'.
c
         p = 6
c
      elseif (method .EQ. 5631) then
c
c        SETUP OF THE COEFFICIENTS FOR THE NEW STAGES
c
c        Here we give coefficients defining the interpolant for a 5-stage,
c        6th order MIRK, stage order 3, with a 6th order interpolant.
c
c        There are four extra stages. The corresponding rows of the RK tableau:
c        7/16 | 7/16 |  1547/32768  -1225/32768  749/4096   -287/2048 -861/16384
c        3/8  | 3/8  |  83/1536      -13/384     283/1536   -167/1536  -49/512
c        9/16 | 9/16 |  1225/32768  -1547/32768  287/2048 -749/4096    861/16384
c        1/8  | 1/8  |  233/3456     -19/1152       0         0           0
c                                            ...   -5/72      7/72      -17/216
c
         s_star = 8
c
         c_star(1) = 1.0d0/2.0d0
         c_star(2) = 1.0d0/2.0d0-sqrt(7.0d0)/14.0d0
         c_star(3) = 87.0d0/100.0d0
c
         v_star(1) = 1.0d0/2.0d0
         v_star(2) = 1.0d0/2.0d0-sqrt(7.0d0)/14.0d0
         v_star(3) = 87.0d0/100.0d0
c
c        (The elements of the matrix containing the stage coefficients are
c        stored columnwise in the vector `x_star'.
c        Column 1 of x_star
c
      x_star(1) = 1.D0/64.D0
      x_star(2) = 3.D0/112.D0+9.D0/1960.D0*dsqrt(7.D0)
      x_star(3) = 2707592511.D0/1000000000000.D0-1006699707.D0/
     #1000000000000.D0*dsqrt(7.D0)
c
c        Column 2 of x_star
c
      x_star(4) = -1.D0/64.D0
      x_star(5) = -3.D0/112.D0+9.D0/1960.D0*dsqrt(7.D0)
      x_star(6) = -51527976591.D0/1000000000000.D0-1006699707.D0/
     #1000000000000.D0*dsqrt(7.D0)
c
c        Column 3 of x_star
c
      x_star(7) = 7.D0/192.D0*dsqrt(21.D0)
      x_star(8) = 11.D0/840.D0*dsqrt(7.D0)+3.D0/112.D0*dsqrt(7.D0)
     #*dsqrt(3.D0)
      x_star(9) = -610366393.D0/75000000000.D0+7046897949.D0/
     #1000000000000.D0*dsqrt(7.D0)+14508670449.D0/1000000000000.D0
     #*dsqrt(7.D0)*dsqrt(3.D0)
c
c        Column 4 of x_star
c
      x_star(10) = -7.D0/192.D0*dsqrt(21.D0)
      x_star(11) = 11.D0/840.D0*dsqrt(7.D0)-3.D0/112.D0*dsqrt(7.D0)
     #*dsqrt(3.D0)
      x_star(12) = -610366393.D0/75000000000.D0+7046897949.D0/
     #1000000000000.D0*dsqrt(7.D0)-14508670449.D0/1000000000000.D0
     #*dsqrt(7.D0)*dsqrt(3.D0)
c
c        Column 5 of x_star
c
      x_star(13) = 0
      x_star(14) = 88.D0/5145.D0*dsqrt(7.D0)
      x_star(15) = -12456457.D0/1171875000.D0+1006699707.D0/
     #109375000000.D0*dsqrt(7.D0)
c
c        Column 6 of x_star
c
      x_star(16) = 0
      x_star(17) = -18.D0/343.D0*dsqrt(7.D0)
      x_star(18) = 3020099121.D0/437500000000.D0*dsqrt(7.D0)+
     #47328957.D0/625000000.D0
c
c        Column 7 of x_star
c
      x_star(19) = 0
      x_star(20) = 0
      x_star(21) = -7046897949.D0/250000000000.D0*dsqrt(7.D0)
c
c        Column 8 of x_star
c
      x_star(22) = 0
      x_star(23) = 0
      x_star(24) = 0
c
c        Define the sample point for the interpolant of the 5-stage, 6th
c        order MIRK discrete formula.
c
         tau_star = 0.4d0
c
c        Define the value of the order of the interpolant, `p'.
c
         p = 6
c
      endif
c
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
          subroutine interp_weights(s_star,tau,w,wp)
c
c***overview***
c
c          This routine evaluates the weight polynomials (and their first
c      derivatives) corresponding to the `s_star' stages used to form the
c      interpolant for the discrete formula defined by the `RK_tableau'
c      routine. The evaluation takes place at `tau'. The results are returned
c      in the `s_star' locations of the vectors `w' and `wp'. The
c      interp_weights routine will be dependent on the parameter, `method'.
c      `method' defines the method to be used and will be accessed via the
c      /rk_method/ common block. The possible values for `method' and the
c      corresponding MIRK method associated with each of these values are
c      listed below.
c
c               Value of `method'   |   MIRK formula
c               -------------------------------------
c                      121          |   MIRK121 scheme
c                      221          |   MIRK221 scheme
c                      232          |   MIRK232 scheme (abscissa = 0,2/3)
c                      2321         |   MIRK232 scheme (abscissa = 1,1/3)
c                      343          |   MIRK343 scheme
c                      453          |   MIRK453 scheme
c                      563          |   MIRK563 scheme
c                      5631         |   An improved MIRK563 scheme
c
c          The presence of the interp_weights routine demonstrates the use
c      of an "information hiding" principle. The more common approach would be
c      to simply provide the coefficients of the polynomial weights in
c      an array made available to the rest of the software package through
c      the /interp_coeffs/ common block. In the approach employed here we
c      "hide" the polynomial expression for each weight inside the
c      `interp_weights' routine whose function is to provide an evaluation
c      of the weights at a given point. An advantage of this approach is
c      that changes can be made in the method by which the weights are
c      determined with no subsequent change needed in the rest of the
c      software.
c----------------------------------------------------------------------------
c***declaration of parameters***
c      imports:
         integer s_star
         double precision tau
c
c        `s_star' is the number of weight polynomials.
c        `tau'    is the point at which the weight polynomials
c                 are to be evaluated.
c      exports:
         double precision w(s_star), wp(s_star)
c
c        `w'    contains the values of the weight polynomials
c               of the interpolant, evaluated at `tau'.
c
c        `wp'   contains the values of the first derivative
c               of the weight polynomials of the interpolant,
c               evaluated  at `tau'.
c
c***declaration of local variables***
         double precision t2, tm1, t4m3, t2m1
c
c        Factors appearing in the expressions for the weight polynomials
c        and derivatives, defined below.
c
c***declaration of variables for the /rk_method/ common block***
c      exports:
         integer   method
         common   /rk_method/ method
c
c        `method'  defines the method to be used
c------------------------------------------------------------------------------
c     called by: `defect_estimate', `interp_eval'
c------------------------------------------------------------------------------
c      The polynomials in `tau' that define the interpolant are very
c      dependent on the underlying discrete formula and the  interpolant
c      derived. The following part of this routine will be different for
c      each formula implemented.
c
       if (method. EQ. 121) then
c
c         WEIGHT POLYNOMIAL IMPLEMENTATION for a 1-stage, 2nd order MIRK formula
c
          w(1) = 0.0d0
c
          w(2) = tau*(1.0d0-tau/2.0d0)
c
          w(3) = tau**2/2.0d0
c
c         Derivative polynomials.
c
          wp(1) = 0.0d0
c
          wp(2) = 1.0d0 - tau
c
          wp(3) = tau
c
       elseif (method. EQ. 221) then
c
c         WEIGHT POLYNOMIAL IMPLEMENTATION for a 2-stage, 2nd order MIRK formula
c
          w(1) = tau*(1.0d0-tau/2.0d0)
c
          w(2) = tau**2/2.0d0
c
c         Derivative polynomials.
c
          wp(1) = 1.0d0 - tau
c
          wp(2) = tau
c
       elseif (method .EQ. 232) then
c
c         WEIGHT POLYNOMIAL IMPLEMENTATION for a 2-stage, 3rd order MIRK formula
c
          w(1) = tau*(2.0d0*tau**2 - 5.0d0*tau + 4.0d0)/4.0d0
c
          w(2) = -3.0d0*tau**2*(2.0d0*tau - 3.0d0)/4.0d0
c
          w(3) = tau**2*(tau - 1.0d0)
c
c         Derivative polynomials.
c
          wp(1) = (tau - 2.0d0/3.0d0)*(tau - 1.0d0)/(2.0d0/3.0d0)
c
          wp(2) = tau*(tau - 1.0d0)/(-2.0d0/9.0d0)
c
          wp(3) = tau*(tau - 2.0d0/3.0d0)/(1.0d0/3.0d0)
c
       elseif (method .EQ. 2321) then
c
c         WEIGHT POLYNOMIAL IMPLEMENTATION for the reflection of the above
c         2-stage, 3rd order MIRK formula
c
          w(1) = 1.0d0/2.0d0*tau**3 - 1.0d0/4.0d0*tau**2
c
          w(2) = -3.0d0/2.0d0*tau**3 + 9.0d0/4.0d0*tau**2
c
          w(3) = tau**3 - 2.0d0*tau**2 + tau
c
c         Derivative polynomials.
c
          wp(1) = 3.0d0*tau*(tau-1.0d0/3.0d0)/2.0d0
c
          wp(2) = -9.0d0*tau*(tau-1.0d0)/2.0d0
c
          wp(3) = 3.0d0*(tau-1.0d0)*(tau-1.0d0/3.0d0)
c
       elseif (method .EQ. 343) then
c
c         WEIGHT POLYNOMIAL IMPLEMENTATION for a 3-stage, 4th order MIRK formula
c
          t2 = tau*tau
          tm1 = tau - 1.0d0
          t4m3 = tau*4.0d0 - 3.0d0
          t2m1 = tau*2.0d0 - 1.0d0
c
          w(1) = -tau*(2.0d0*tau-3.0d0)*(2.0d0*t2-3.0d0*tau+2.0d0)/6.0d0
c
          w(2) = t2*(12.0d0*t2-20.0d0*tau + 9.0d0)/6.0d0
c
          w(3) = 2.0d0*t2*(6.0d0*t2-14.0d0*tau+9.0d0)/3.0d0
c
          w(4) = -16.0d0*t2*tm1*tm1/3.0d0
c
c         Derivative polynomials.
c
          wp(1) = -tm1*t4m3*t2m1/3.0d0
c
          wp(2) = tau*t2m1*t4m3
c
          wp(3) = 4.0d0*tau*t4m3*tm1
c
          wp(4) = -32.0d0*tau*t2m1*tm1/3.0d0
c
      elseif (method .EQ. 453) then
c
      w(1) = tau*(22464.0d0-83910.0d0*tau+143041.0d0*tau**2
     + -113808.0d0*tau**3+33256.0d0*tau**4)/22464.0d0
c
      w(2) = tau**2*(-2418.0d0+12303.0d0*tau-19512.0d0*tau**2
     ++10904.0d0*tau**3)/3360.0d0
c
      w(3) = -8.0d0/81.0d0*tau**2*(-78.0d0+209.0d0*tau-204.0d0
     +*tau**2+8.0d0*tau**3)
c
      w(4) = -25.0d0/1134.0d0*tau**2*(-390.0d0+1045.0d0*tau
     +-1020.0d0*tau**2+328.0d0*tau**3)
c
      w(5) = -25.0d0/5184.0d0*tau**2*(390.0d0+255.0d0*tau-
     +1680.0d0*tau**2+2072.0d0*tau**3)
c
      w(6) = 279841.0d0/168480.0d0*tau**2*(-6.0d0+21.0d0*tau
     +-24.0d0*tau**2+8.0d0*tau**3)
c
c     Derivative polynomials.
c
      wp(1) = 1.0d0-13985.0d0/1872.0d0*tau+143041.0d0/7488.0d0
     +*tau**2-2371.0d0/117.0d0*tau**3+20785.0d0/2808.0d0*tau**4
c
      wp(2) = -403.0d0/280.0d0*tau+12303.0d0/1120.0d0*tau**2
     +-813.0d0/35.0d0*tau**3+1363.0d0/84.0d0*tau**4
c
      wp(3) = 416.0d0/27.0d0*tau-1672.0d0/27.0d0*tau**2+2176.0d0
     +/27.0d0*tau**3-320.0d0/81.0d0*tau**4
c
      wp(4) = 3250.0d0/189.0d0*tau-26125.0d0/378.0d0*tau**2+
     +17000.0d0/189.0d0*tau**3-20500.0d0/567.0d0*tau**4
c
      wp(5) = -1625.0d0/432.0d0*tau-2125.0d0/576.0d0*tau**2
     ++875.0d0/27.0d0*tau**3-32375.0d0/648.0d0*tau**4
c
      wp(6) = -279841.0d0/14040.0d0*tau+1958887.0d0/18720.0d0*
     +tau**2-279841.0d0/1755.0d0*tau**3+279841.0d0/4212.0d0*tau**4
c
      elseif (method .EQ. 4531) then
c
      w(1) = tau*(3857490.0d0-14075970.0d0*tau+24780520.0d0
     +*tau**2-20977215.0d0*tau**3+7479584.0d0*tau**4)/3857490.0d0
c
      w(2) = tau**2*(-959595.0d0+3459675.0d0*tau-4676415.0d0
     +*tau**2+1893368.0d0*tau**3)/989100.0d0
c
      w(3) = -32.0d0/63585.0d0*tau**2*(-13650.0d0+42100.0d0*tau
     +-47175.0d0*tau**2+10376.0d0*tau**3)
c
      w(4) = -250.0d0/89019.0d0*tau**2*(-2730.0d0+8420.0d0*tau
     +-9435.0d0*tau**2+4336.0d0*tau**3)
c
      w(5) = 28561.0d0/791280.0d0*tau**2*(390.0d0+255.0d0*tau
     +-1680.0d0*tau**2+2072.0d0*tau**3)
c
      w(6) = -279841.0d0/2449200.0d0*tau**2*(210.0d0-225.0d0*
     +tau-180.0d0*tau**2+536.0d0*tau**3)
c
      wp(1) = 1.0d0-938398.0d0/128583.0d0*tau+2478052.0d0/128583.0d0
     +*tau**2-399566.0d0/18369.0d0*tau**3+534256.0d0/55107.0d0*tau**4
c
      wp(2) = -9139.0d0/4710.0d0*tau+46129.0d0/4396.0d0*tau**2-
     +311761.0d0/16485.0d0*tau**3+473342.0d0/49455.0d0*tau**4
c
      wp(3) = 58240.0d0/4239.0d0*tau-269440.0d0/4239.0d0*tau**2+
     +402560.0d0/4239.0d0*tau**3-332032.0d0/12717.0d0*tau**4
c
      wp(4) = 65000.0d0/4239.0d0*tau-2105000.0d0/29673.0d0*tau**2
     ++3145000.0d0/29673.0d0*tau**3-5420000.0d0/89019.0d0*tau**4
c
      wp(5) = 371293.0d0/13188.0d0*tau+485537.0d0/17584.0d0*tau**2
     +-114244.0d0/471.0d0*tau**3+1056757.0d0/2826.0d0*tau**4
c
      wp(6) = -1958887.0d0/40820.0d0*tau+2518569.0d0/32656.0d0*
     +tau**2+839523.0d0/10205.0d0*tau**3-18749347.0d0/61230.0d0*tau**4
c
      elseif (method .EQ. 563) then
c
c     WEIGHT POLYNOMIAL IMPLEMENTATION for a 5-stage, 6th order MIRK formula
c
      w(1) =
     +tau-28607.D0/7434.D0*tau**2-166210.D0/33453.D0*tau**3+334780.D0/1
     +1151.D0*tau**4-1911296.D0/55755.D0*tau**5+406528.D0/33453.D0*tau*
     +*6
c
      w(2) =
     +777.D0/590.D0*tau**2-2534158.D0/234171.D0*tau**3+0.208858D7/78057
     +.D0*tau**4-10479104.D0/390285.D0*tau**5+0.11328512D8/0.1170855D7*
     +tau**6
c
      w(3) =
     +-1008.D0/59.D0*tau**2+222176.D0/1593.D0*tau**3-180032.D0/531.D0
     +*tau**4+876544.D0/2655.D0*tau**5-180224.D0/1593.D0*tau**6
c
      w(4) =
     +-1008.D0/59.D0*tau**2+222176.D0/1593.D0*tau**3-180032.D0
     +/531.D0*tau**4+876544.D0/2655.D0*tau**5-180224.D0/1593.D0*tau**6
c
      w(5) =
     +-378.D0/59.D0*tau**2+27772.D0/531.D0*tau**3-22504.D0/177.D0*tau*
     +*4+109568.D0/885.D0*tau**5-22528.D0/531.D0*tau**6
c
      w(6) =
     +-95232.D0/413.D0*tau**2+0.62384128D8/33453.D0*tau**3-49429504.D0
     +/11151.D0*tau**4+0.46759936D8/11151.D0*tau**5-46661632.D0/33453.
     +D0*tau**6
c
      w(7) =
     +896.D0/5.D0*tau**2-4352.D0/3.D0*tau**3+3456*tau**4-16384.D0/5.D0
     +*tau**5+16384.D0/15.D0*tau**6
c
      w(8) =
     +50176.D0/531.D0*tau**2-179554304.D0/234171.D0*tau**3+0.143363072D9
     +/78057.D0*tau**4-136675328.D0/78057.D0*tau**5+0.137363456D9
     +/234171.D0*tau**6
c
      w(9) =
     +16384.D0/441.D0*tau**3-16384.D0/147.D0*tau**4+16384.D0/147.D0
     +*tau**5-16384.D0/441.D0*tau**6
c
c     Derivative polynomials.
c
      wp(1) =
     +1-28607.D0/3717.D0*tau-166210.D0/11151.D0*tau**2+0.133912D7/
     +11151.D0*tau**3-1911296.D0/11151.D0*tau**4+813056.D0/11151.
     +D0*tau**5
c
      wp(2) =
     +777.D0/295.D0*tau-2534158.D0/78057.D0*tau**2+0.835432D7
     +/78057.D0*tau**3-10479104.D0/78057.D0*tau**4+0.22657024D8
     +/390285.D0*tau**5
c
      wp(3) =
     +-2016.D0/59.D0*tau+222176.D0/531.D0*tau**2-720128.D0/531.D0
     +*tau**3+876544.D0/531.D0*tau**4-360448.D0/531.D0*tau**5
c
      wp(4) =
     +-2016.D0/59.D0*tau+222176.D0/531.D0*tau**2-720128.D0/531.D0
     +*tau**3+876544.D0/531.D0*tau**4-360448.D0/531.D0*tau**5
c
      wp(5) =
     +-756.D0/59.D0*tau+27772.D0/177.D0*tau**2-90016.D0/177.D0*
     +tau**3+109568.D0/177.D0*tau**4-45056.D0/177.D0*tau**5
c
      wp(6) =
     +-190464.D0/413.D0*tau+0.62384128D8/11151.D0*tau**2-1977
     +18016.D0/11151.D0*tau**3+0.23379968D9/11151.D0*tau**4-
     +93323264.D0/11151.D0*tau**5
c
      wp(7) =
     +1792.D0/5.D0*tau-4352*tau**2+13824*tau**3-16384*tau**4
     ++32768.D0/5.D0*tau**5
c
      wp(8) =
     +100352.D0/531.D0*tau-179554304.D0/78057.D0*tau**2+0.57345
     +2288D9/78057.D0*tau**3-683376640.D0/78057.D0*tau**4+0.274
     +726912D9/78057.D0*tau**5
c
      wp(9) =
     +16384.D0/147.D0*tau**2-65536.D0/147.D0*tau**3+81920.D0/147
     +.D0*tau**4-32768.D0/147.D0*tau**5
c
      elseif (method .EQ. 5631) then
c
c     WEIGHT POLYNOMIAL IMPLEMENTATION for an improved 5-stage, 6th order
c     MIRK formula
c
      w(1) =
     #-(12233+1450*dsqrt(7.D0))*(800086000*tau**5+63579600*dsqrt
     #(7.D0)*tau**4-2936650584.D0*tau**4+4235152620.D0*tau**3
     #-201404565*dsqrt(7.D0)*tau**3+232506630*dsqrt(7.D0)*tau**2
     #-3033109390.D0*tau**2+1116511695*tau-116253315*dsqrt(7.D0)
     #*tau+22707000*dsqrt(7.D0)-191568780)*tau/2112984835740.D0
c
      w(2) =
     #-(-10799+650*dsqrt(7.D0))*(24962000*tau**4
     #+473200*dsqrt(7.
     #D0)*tau**3-67024328*tau**3-751855*dsqrt(7.D0)*tau**2
     #+66629600*tau**2-29507
     #250*tau+236210*dsqrt(7.D0)*tau+5080365+50895*dsqrt(7.D0))
     #*tau**2/29551834260.D0
c
      w(3) =
     #7.D0/1274940.D0*(259+50*dsqrt(7.D0))
     #*(14000*tau**4-48216*tau
     #**3+1200*dsqrt(7.D0)*tau**3-3555*dsqrt(7.D0)
     #*tau**2+62790*tau**2+3610*ds
     #qrt(7.D0)*tau-37450*tau+9135-1305*dsqrt(7.D0))*tau**2
c
      w(4) =
     #7.D0/1274940.D0*(259+50*dsqrt(7.D0))*(14000*tau**4
     #-48216*tau
     #**3+1200*dsqrt(7.D0)*tau**3-3555*dsqrt(7.D0)
     #*tau**2+62790*tau**2+3610*ds
     #qrt(7.D0)*tau-37450*tau+9135-1305*dsqrt(7.D0))*tau**2
c
      w(5) =
     #16.D0/2231145.D0*(259+50*dsqrt(7.D0))
     #*(14000*tau**4-48216*
     #tau**3+1200*dsqrt(7.D0)*tau**3-3555*dsqrt(7.D0)*tau**2
     #+62790*tau**2+3610*d
     #sqrt(7.D0)*tau-37450*tau+9135
     #-1305*dsqrt(7.D0))*tau**2
c
      w(6) =
     #4.D0/1227278493.D0*(740*dsqrt(7.D0)-6083)*(1561000*tau**2-
     #2461284*tau-109520*dsqrt(7.D0)*tau+979272+86913*
     #dsqrt(7.D0))*(tau-1)**2*tau**2
c
      w(7) =
     #-49.D0/63747.D0*dsqrt(7.D0)*(20000*tau**2
     #-20000*tau+3393)*(tau-1)**2*tau**2
c
      w(8) =
     #-1250000000.D0/889206903.D0*(28*tau**2-28*tau+9)
     #*(tau-1)**2*tau**2
c
c     Derivative polynomials.
c
      wp(1) =
     #(1450*dsqrt(7.D0)+12233)*(14*tau-7
     #+dsqrt(7.D0))*(tau-1)*(-4
     #00043*tau+75481+2083*dsqrt(7.D0))*(100*tau-87)
     #*(2*tau-1)/493029795006.D0
c
      wp(2) =
     # -(650*dsqrt(7.D0)-10799)*(14*tau-7+dsqrt(7.D0))*(37443*tau-
     #13762-2083*dsqrt(7.D0))*(100*tau-87)
     #*(2*tau-1)*tau/20686283982.D0
c
      wp(3) =
     # 7.D0/42498.D0*(259+50*dsqrt(7.D0))*(14*tau-7+dsqrt(7.D0))
     #*(tau-1)*(100*tau-87)*(2*tau-1)*tau
c
      wp(4) =
     # 7.D0/42498.D0*(259+50*dsqrt(7.D0))*(14*tau-7+dsqrt(7.D0))
     #*(tau-1)*(100*tau-87)*(2*tau-1)*tau
c
      wp(5) =
     # 32.D0/148743.D0*(259+50*dsqrt(7.D0))*(14*tau-7+dsqrt(7.D0
     #))*(tau-1)*(100*tau-87)*(2*tau-1)*tau
c
      wp(6) =
     # 4.D0/1227278493.D0*(740*dsqrt(7.D0)-6083)*(14*tau-7+dsqrt
     #(7.D0))*(tau-1)*(100*tau-87)*(6690*tau-4085
     #-869*dsqrt(7.D0))*tau
c
      wp(7) =
     # -98.D0/21249.D0*dsqrt(7.D0)*(tau-1)*(100*tau-13)
     #*(100*tau-87)*(2*tau-1)*tau
c
      wp(8) =
     # -1250000000.D0/2074816107.D0*(14*tau-7+dsqrt(7.D0))
     #*(tau-1)*(14*tau-7-dsqrt(7.D0))*(2*tau-1)*tau
c
      endif
c
      return
      end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine interval(Nsub, mesh, t, i)
c
c***overview***
c
c     This routine finds `i' such that `mesh(i-1) =< t < mesh(i)' for
c     `t < mesh(Nsub)'. The value of `Nsub' is returned when `t = mesh(Nsub)'.
c     `Nsub' is the number of subintervals and `mesh' is the list of points
c     defining the subintervals.
c----------------------------------------------------------------------------
c***declaration of parameters***
c   imports:
      integer Nsub
      double precision mesh(0:Nsub), t
c
c     `Nsub' Number of subintervals.
c     `mesh' Array of points defining the subintervals.
c     `t' Point to be located within `mesh'.
c
c   exports:
      integer i
c
c     `i' Mesh index such that `t' satisfies `mesh(i-1) =< t < mesh(i)'
c         for `t < mesh(Nsub)' or `i = Nsub for t = mesh(Nsub)'.
c
c***declaration of local variables***
      double precision     meshhalfnew
      integer              inew, iold, iright, ileft
c
c     `meshhalfnew' New mesh value resulting from halving the current interval.
c     `inew' Current subinterval containing `meshhalfnew'.
c     `iold' Previous subinterval containing `meshhalfnew'.
c     `ileft' Left hand mesh point bounding `meshhalfnew'.
c     `iright' Right hand mesh point bounding `meshhalfnew'.
c---------------------------------------------------------------------------
c      Called by: `interp_eval' or `sol_eval'.
c---------------------------------------------------------------------------
c     Handle special cases for the location of `t' first
c
      if ((t .LT. mesh(0)) .OR. (t .GT. mesh(Nsub))) then
         write(6,*)'Error-t is outside of problem interval'
      endif
c
      if (t .EQ. mesh(Nsub)) then
                  i = Nsub
                  return
      end if
c
      if (t .EQ. mesh(0)) then
                  i = 1
                  return
      end if
c
c     At this point we know that `mesh(0) < t < mesh(Nsub)'.
c     Binary search to determine subinterval containing `t'.
c
      inew = (Nsub/2)
      meshhalfnew = mesh(inew)
c
c     Determine if `t' is equal to `meshhalfnew'.
c
      if (t .EQ. meshhalfnew) then
         i = inew + 1
         return
      endif
c
c     `t' is not equal to the meshpoint.
c
      if (t .LT. meshhalfnew) then
c
c          `t' is on the left hand side of the interval.
c
           dowhile (t .LT. meshhalfnew)
              iold = inew
              inew = iold/2
              meshhalfnew = mesh(inew)
           enddo
c
           iright = iold
           ileft = inew
c
      else
c
c          `t' is on the right hand side of the interval.
c
           dowhile (t .GT. meshhalfnew)
              iold = inew
              inew = nint(real((iold + (Nsub - iold)/2)) + 0.5)
              meshhalfnew = mesh(inew)
           enddo
c
           ileft  = iold
           iright = inew
c
      endif
c
      i = iright
c
c     Determine if `t' is equal to `meshhalfnew'.
c
      if (t .EQ. meshhalfnew) then
         i = inew + 1
         return
      endif
c
c     `t' is not equal to the meshpoint.
c
c     `t' is bounded by `ileft' and `iright'.
c
      dowhile (abs(iright - ileft) .NE. 1)
           inew = (ileft + iright)/2
           if(mesh(inew) .LT. t) then
                 ileft = inew
           elseif (mesh(inew) .GT. t) then
                 iright = inew
           else
                 i = inew + 1
                 return
           endif
      enddo
c
      i = iright
c
c     `t' is in the `i'-th subinterval.
c
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine jacblk(neqns,hi,ti,yi,yip1,k,LLi,RRi,
     +                            yr,dfdy,ur,dkdyi,dkdyip1,dfsub)
c
c***overview***
c
c       ****Evaluation of the Jacobian for one subinterval****
c
c           This routine computes the left and right blocks (of the Jacobian of
c       the residual function) associated with the i-th subinterval. The left
c       and right block values are returned in the parameters `LLi' and `RRi'.
c       Each block is of size `neqns**2', where `neqns' is the number of
c       differential equations. We have,
c
c                         s                               s
c                         --                              --
c       LLi =  -I - hi * >  br * dkr/dyi,  RRi = I - hi * >  br * dkr/dyi+1
c                         --                              --
c                         r=1                             r=1
c
c      Each partial derivative, dkr/dyi, is computed as follows,
c
c                                          r-1
c                                          --
c      dkr/dyi = df/dy(tr,yr)*((1-vr)*I+hi*>  xrj * dkj/dyi)
c                                          --
c                                          j=1,
c
c      with a similar expression for each dkr/dyi+1.
c
c         The solution value, yr, at which df/dy is evaluated depends on
c      the input parameters `yi' and `yip1', the solution approximations
c      at each end of the ith subinterval, and the previously computed
c      stages provided by the input parameter `k'. The size of the ith
c      subinterval is `hi' and the endpoint is `ti'. The computation of
c      `yr' also depends on the coefficients of the RK method provided
c      through the /RK_s/,/RK_coeff/ common blocks.
c
c            (All matrices employed in this routine are actually stored
c      columnwise within vectors.)
c--------------------------------------------------------------------------------
c***declaration of constants***
            integer Mxs
            parameter(Mxs=10)
c
c            `Mxs' is the maximum value for the total number of stages
c             of the interpolant.
c
c***declaration of parameters***
c         imports:
             integer                   neqns
             double precision          hi, ti, yi(neqns),yip1(neqns)
             double precision          k(Mxs*neqns)
c
c            `neqns' is the number of differential equations.
c
c            `hi'    is the size of the current subinterval.
c
c            `ti'    is the value of the lefthand endpoint of the
c                       ith subinterval.
c
c            `yi'    is the current solution approximation at `ti'.
c
c            `yip1'  is the current solution approximation at `ti+h', i.e.
c                    at the right hand endpoint of the subinterval.
c
c            `k'     contains the `s' stage evaluations, each of size
c                       `neqns'associated with the `ith' subinterval.
c         exports:
              double precision          LLi(neqns**2),RRi(neqns**2)
c
c            `LLi' is the partial derivative of the ith component of
c                  the residual function with respect to yi
c
c            `RRi' is the partial derivative of the ith component of the
c                  residual function with respect to yi+1
c
c***declaration of work arrays***
            double precision        yr(neqns)
            double precision        dfdy(neqns**2)
            double precision        ur(neqns**2)
            double precision        dkdyi(neqns**2*Mxs)
            double precision        dkdyip1(neqns**2*Mxs)

c            `yr'      the argument for the evaluation of the `rth' stage.
c
c            `dfdy'    during the `rth' time through the main loop below, it
c                      contains the derivative of the right hand side of the
c                      system of ODE's, evaluated at the arguments of the
c                      'rth' stage of the RK scheme.
c
c            `ur'      used to store the intermediate values arising in the
c                      computation of each partial derivative.
c
c            `dkdyi'   The `rth' set of `neqns**2' locations of this vector
c                      will contain the value of the rth partial derivative
c                      with respect to `yi'.
c
c            `dkdyip1' The `rth' set of `neqns**2' locations of this vector
c                      will contain the value of the rth partial derivative
c                      with respect to `yip1'.
c
c***user - supplied subroutines***
            external dfsub
c
c           dfsub - This routine defines the Jacobian of the system of
c                   differential equations with respect to y, i.e.
c                   df/dy.
c
c***declaration of local variables***
            integer                 l,i,m,r, j,i1,i2
            double precision        sum
            integer                 neqnssq
            integer                 alpha,beta,sigma
            double precision        tr
            logical                 dkdyi_zero(Mxs)
            logical                 dkdyip1_zero(Mxs)
            logical                 matrix_mult_dkdyi
            logical                 matrix_mult_dkdyip1
c
c            `l,i,m,r,j'     are loop index variables
c
c            `sum'     is an accumulator
c
c            `i1',`i2' are indexes used in calculating ur and yr
c
c            `neqnssq' is the value of neqns squared
c
c            `alpha',`beta',`sigma' are indexes used in the
c                     multiplication of the two matrices `ur' and `dfdy'
c
c            `tr'     the abscissa for the evaluation of the `rth' stage.
c
c            `dkdyi_zero' is a logical array. If the computation of the rth
c                         stage has been completed then a value of `.TRUE.' in
c                         `dkdyi_zero(r)' indicates that every component of
c                         the rth partial derivative with respect to `yi'
c                         is zero. If the computation of the rth stage has
c                         not been completed then a value of `.TRUE.'
c                         indicates that the `ur' components of the rth
c                         partial derivative are currently zero. A value of
c                         `.FALSE.' in `dkdyi_zero(r)' implies the rth
c                         partial derivative is non-zero in at least one of
c                         its components.
c
c            `dkdyip1_zero' is a logical array. If the computation of the rth
c                           stage has been completed then a value of `.TRUE.'
c                           in `dkdyip1_zero(r)' indicates that every component
c                           of the rth partial derivative with respect to
c                           `yip1' is zero. If the computation of the rth
c                           stage has not been completed then a value of
c                           `.TRUE.' indicates that the `ur' components of
c                           the rth partial derivative are currently zero.
c                           A value of `.FALSE.' in `dkdyip1_zero(r)'
c                           implies the rth partial derivative is non-zero
c                           in at least one of its components. `dkdyi_zero'
c                           and `dkdyip1_zero' are both used to avoid
c                           multiplications by zero.
c
c            `matrix_mult_dkdyi' is a logical variable. For the compuation
c                                of a given stage, a value of `.FALSE.'
c                                indicates that  the partial derivative
c                                with respect to `yi'
c                                has the form `(1-v(r))*dfdy' and that
c                                multiplying the matrix `ur' by the matrix
c                                `dfdy' is unnecessary. A value of `.TRUE.'
c                                indicates that the partial derivative
c                                does not have the form `(1-v(r))*dfdy' and
c                                that multiplying the matrix `ur' by the
c                                matrix `dfdy' may be necessary.
c
c            `matrix_mult_dkdyip1' is a logical variable. For the compuation
c                                  of a given stage, a value of `.FALSE.'
c                                  indicates that  the partial derivative
c                                  with respect to `yip1'
c                                  has the form `v(r)*dfdy' and that
c                                  multiplying the  matrix `ur' by the matrix
c                                  `dfdy' is unnecessary. A value of `.TRUE.'
c                                  indicates that the partial derivative
c                                  does not have the form `v(r)*dfdy' and
c                                  that multiplying the matrix `ur' by the
c                                  matrix `dfdy' may be necessary.
c
c***declaration of variables in common block /RK_s/***
            integer                 s
            common /RK_s/ s
c
c           `s' is the number of discrete stages of the Runge-Kutta formula.
c
c***declaration of variables in common block /RK_coeff/***
            double precision        v(Mxs), b(Mxs)
            double precision        c(Mxs), x(Mxs**2)
            common /RK_coeff/ v, b, c, x
c
c           `v', `b', `c',and `x' are coefficients that define
c           the discrete Runge-Kutta formula.
c
c           The first component of `c' contains the abscissa for the first
c           new stage. The second component of `c' contains the abscissa for
c           the second new stage, and so on.
c
c           The first component of `v' contains the blending coefficient
c           for the first new stage. The second component of `v' contains
c           the blending coefficient for the second new stage, and so on.
c
c           The matrix, X, is related to  the data structure, `x'
c           as follows (the matrix is stored columnwise in the vector):
c
c                 X(i,j) = x( (j-1)*s + i )
c
c           `x' contains the coupling coefficients for each of the new
c           stages
c
c------------------------------------------------------------------------------
c            called by: `newmat'
c            calls to: `dfsub'
c------------------------------------------------------------------------------
c***COMPUTATION OF THE PARTIAL DERIVATIVES***
c
c     We compute the partial derivatives in pairs: dkr/dyi and dkr/dyi+1
c     during the `rth' pass of the following loop. The results are stored
c     in the appropriate locations of the vectors `dkdyi' and `dkdyip1',
c     respectively.
c
c     Initialize the components `dkdyi_zero' and `dkdyi_zero' to `.TRUE.'
c     indicating that all of the partial derivatives are currently zero.
c
      do 1 j = 1,s
        dkdyi_zero(j)   = .TRUE.
        dkdyip1_zero(j) = .TRUE.
 1    continue
c
      neqnssq = neqns**2
c
      do 40 r = 1, s
c
c         Initialize `matrix_mult_dkdyi' to `.TRUE.' indicating that the
c         partial derivative is not of the form `(1-v(r))*dfdy' and that
c         multiplying the matrix `ur' by `dfdy' may still be necessary.
c
          matrix_mult_dkdyi = .TRUE.
c
c         Construct the arguments for the `rth' stage, needed here for the
c         evaluation of `df/dy'
c
c         Compute the abscissa
c
          tr = ti + c(r)*hi
c
c         Begin the construction of `yr' by setting it to zero.
c
          do 2 j=1, neqns
            yr(j) = 0.0d0
 2        continue
c
c          Add on weighted contributions from first `r-1' stages.
c          The `rth' stage of the current subinterval is stored in
c          `k((r-1)*neqns+1)'.
c
           do 5 j = 1, r-1
c
            i1 = (j-1)*s + r
            i2 = (j-1)*neqns
c
            if ( x(i1) .NE. 0.0d0 ) then
c
c             If `x(i1) = 0' then there is no need to perform the following
c             multiplication.

              do 3 l = 1, neqns
                yr(l) = yr(l) + x(i1)*k(i2+l)
 3            continue
c
             endif
c
 5        continue
c
c          Multiply yr by hi.
c
           do 7 l =1, neqns
              yr(l) = hi*yr(l)
 7         continue
c
c          Add on weighted contributions from `yi' and `yip1'.
c
           do 10 l=1,neqns
              yr(l) = (1.0d0 - v(r))*yi(l) + v(r)*yip1(l) + yr(l)
 10        continue
c
c          The arguments for the `rth' stage, `tr' and `yr' have been
c          constructed. Evaluate df/dy, the derivative of the system of ODE's
c          by first setting it equal to zero then calling the user defined
c          subroutine `dfsub' at these arguments to fill in the non-zero
c          entries.
c
            do 12 i = 1,neqnssq
               dfdy(i) = 0.0d0
 12         continue
c
            call dfsub(neqns,tr,yr,dfdy)
c
c          The values of the `rth' partial derivatives with respect to `yi'
c          and `yip1' are obtained by multiplying `dfdy' by a weighted linear
c          combination of the partial derivatives of the previous `r-1'
c          partial derivatives. We will store these linear combinations, as we
c          construct them, in the matrix, `ur'.
c
c          Construct the `ur' matrix for dkr/dyi
c
           do 14 j = 1, neqnssq
                 ur(j) = 0.0d0
 14        continue
c
c          Add on weighted multiples of the previous partial derivatives.
c          The `jth' previous partial derivative is stored in
c          `dkdyi((j-1)*neqns**2+1)'. This can be done through a call to
c          `daxpy' since we are storing all matricies in vector form.
c
           do 15 j = 1, r-1
c
            i1 = (j-1)*s + r
            i2 = (j-1)*neqnssq
c
            if ( x(i1) .NE. 0  .AND. ( .NOT. dkdyi_zero(j) ) ) then
c
c             If `x(i1) = 0' or the jth partial derivative equals zero then
c             there is no need to perform the following multiplication.
c
              dkdyi_zero(r) = .FALSE.

c             Indicates the matrix `ur' is non zero in at least one of the
c             components and we will require a matrix multiply later.
c             It also indicates that the rth partial derivative of
c             `dkdyi' will be non zero in one of the components.
c
              do 16 l= 1, neqnssq
                 ur(l) = ur(l) + x(i1)*dkdyi(i2+l)
 16           continue
c
            endif
c
 15        continue
c
c          Multiply ur by hi
c
           if ( .NOT. dkdyi_zero(r) ) then
c
c            If the `ur' matrix for the rth partial derivative equals zero
c            then there is no need to perform the following multiplication.

             do 18 j = 1, neqnssq
               ur(j) = hi*ur(j)
 18          continue
c
           endif
c
c          Add on weighted contribution (1.0-v(r))*I.
c
           alpha = (r-1)*neqns**2

           if ( v(r) .NE. 1.0d0 ) then
c
c            If `1-v(r) = 0' then there is no need to perform the following
c            multiplication.
c
             if ( dkdyi_zero(r) ) then

               matrix_mult_dkdyi = .FALSE.

c              Signifies that multiplying the matrix `ur' by the
c              matrix `dfdy' is not necessary later on
c              because the matrix dkdyi can be obtained simply
c              by multiplying the scalar (1-v(r)) times the matrix dfdy.
c              We proceed to perform the scalar multiplication here.
c
               do j = 1,neqnssq
                 dkdyi(alpha + j) = ( 1.0d0-v(r) )*dfdy(j)
               enddo
c
             else

               do 20 j =1, neqnssq, neqns+1
                 ur(j) = 1.0d0-v(r) + ur(j)
 20            continue
c
             endif
c
             dkdyi_zero(r) = .FALSE.
c
c            Indicates the matrix `ur' is non zero in at least one of the
c            components and we will require a matrix multiply later.
c            It also indicates that the rth partial derivative of
c            `dkdyi' will be non zero in one of the components.
c
           endif
c
c           `ur' now contains the proper linear combination of previously
c           computed partial derivatives with respect to `yi'.
c           The value of dkr/dyi is obtained simply by multiplying `ur' and
c           `dfdy' together. The resultant matrix value is stored in
c           the `neqns*neqns' locations of `dkdyi((r-1)*neqns**2+1)'.
c
c           Multiply `ur' and `dfdy'
c
           if ( matrix_mult_dkdyi ) then
c
            if ( .NOT. dkdyi_zero(r) ) then
c
c            If the `ur' matrix for the rth partial derivative equals zero
c            or rth partial derivative could be obtained by multiplying
c            the scalar `(1-v(r)' by `dfdy' then there is no need to perform
c            the following matrix multiplication.
c
             do 26 i =1, neqns
                 beta = alpha + i
c
                 do 27 j = 1, neqns
                   sigma = (j-1)*neqns
                   sum = 0.0d0
c
                   do 28 m = 1,neqns
                      sum = sum + dfdy((m-1)*neqns+i)*ur(sigma+m)
  28               continue
c
                   dkdyi(sigma+beta) = sum
  27             continue
c
  26         continue
c
            else

c             If the `ur' matrix is equal to zero then the corresponding rth
c             partial derivative is set to zero. As well, the value of
c             `dkdyi_zero(r)' will left as `.TRUE.' to avoid
c             multiplying by zero in subsequent stage computations.
c
              do 25 j = 1,neqnssq
               dkdyi(alpha + j) = 0.0d0
 25           continue
c
            endif
c                 `(.NOT. dkdyi_zero(r) )'
c
           endif
c               `( matrix_mult_dkdyi) '
c
c          Initialize `matrix_mult_dkdyip1' to `.TRUE.' indicating that the
c          partial derivative is currently not of the form `v(r)*dfdy' and
c          that multiplying the matrix `ur' by `dfdy' may still be necessary.
c
            matrix_mult_dkdyip1 = .TRUE.
c
c           Construct the `ur' matrix for dkr/dyip1
c
            do 24 j = 1, neqnssq
                  ur(j) = 0.0d0
 24         continue
c
c           Add on weighted multiples of the previous partial derivatives.
c           The `jth' previous partial derivative is stored in
c           `dkdyip1((j-1)*neqns**2+1)'.
c
            do 30 j = 1, r-1
c
            i1 = (j-1)*s+ r
            i2 = (j-1)*neqnssq
c
            if ( x(i1) .NE. 0  .AND.( .NOT. dkdyip1_zero(j)) ) then

              dkdyip1_zero(r) = .FALSE.

c             Indicates the matrix `ur' is non zero in at least one of the
c             components and we will require a matrix multiply later.
c             It also indicates that the rth partial derivative of
c             `dkdyip1' will be non zero in one of the components.
c
              do 31 l = 1, neqnssq
                ur(l) = ur(l) + x(i1)*dkdyip1(i2+l)
 31           continue
c
             endif
c
 30         continue
c
c           Multiply ur by hi
c
            if ( .NOT. dkdyip1_zero(r) ) then
c
c            If the `ur' matrix for the rth partial derivative equals zero
c            then there is no need to perform the following multiplication.

             do 35 j = 1, neqnssq
               ur(j) = hi*ur(j)
 35          continue
c
            endif
c
c           Add on weighted contribution v(r)*I
c
            if ( v(r) .NE. 0.0d0 ) then
c
c            If `v(r) = 0' then there is no need to perform the following
c            multiplication.
c
               if ( dkdyip1_zero(r) ) then
c
                matrix_mult_dkdyip1 = .FALSE.

c               Signifies that multiplying the matrix `ur' by the
c               matrix `dfdy' is not necessary later on
c               because the matrix dkdyi can be obtained simply
c               by multiplying the scalar v(r) times the matrix dfdy.
c               We proceed to perform the scalar multiplication here.
c
                do j = 1,neqnssq
                  dkdyip1(alpha + j) = v(r) * dfdy(j)
                enddo
c
               else

                do 36 j = 1, neqnssq, neqns+1
                  ur(j) = v(r) + ur(j)
 36             continue
c
               endif
c
               dkdyip1_zero(r) = .FALSE.
c
c              Indicates the matrix `ur' is non zero in at least one of the
c              components and we will require a matrix multiply later.
c              It also indicates that the rth partial derivative of
c              `dkdyip1' will be non zero in one of the components.
c
            endif
c
c           `ur' now contains the proper linear combination of previously
c           computed partial derivatives with respect to `yip1'.
c           The value of dkr/dyi+1 is obtained simply by multiplying `ur'
c           and `dfdy' together. The resultant matrix value is stored in
c           the `neqns*neqns' locations of `dkdyip1((r-1)*neqns**2+1)'.
c
c          Multiply `ur' by `dfdy'
c
          if ( matrix_mult_dkdyip1 ) then

           if ( .NOT. dkdyip1_zero(r) ) then

c            If the `ur' matrix for the rth partial derivative equals zero
c            or rth partial derivative could be obtained by multiplying
c            the scalar `v(r)' by `dfdy' then there is no need to perform
c            the following matrix multiplication.
c
             do 39 i = 1, neqns
                   beta = alpha + i
c
                   do 38 j = 1 ,neqns
                    sum=0.0d0
                    sigma = (j-1)*neqns
c
                    do 37 m = 1, neqns
                     sum = sum + dfdy((m-1)*neqns+i)*ur(sigma+m)
 37                 continue
c
                    dkdyip1(sigma+beta) = sum
 38                continue
c
 39          continue
c
            else
c
c            If the `ur' matrix is equal to zero then the corresponding rth
c            partial derivative is set to zero. As well, the value of
c            `dkdyip1_zero(r)' will left as `.TRUE.' to avoid
c            multiplying by zero in subsequent stage computations.
c
             do 41 j = 1,neqnssq
              dkdyip1(alpha + j) = 0.0d0
 41          continue
c
            endif
c               `( .NOT. dkdyip1_zero(r) )'
c
          endif
c               `( matrix_mult_dkdyip1 )'
c
 40     continue
c
c      All partial derivatives have now been computed.
c
c      COMPUTATION OF THE BLOCKS, LLi AND RRi.
c
c      The left and right blocks are now constructed from linear combinations
c      of the partial derivatives.
c
c      Initialize LLi and RRi to 0
c
            do 44 j = 1, neqnssq
                  LLi(j) = 0.0d0
                  RRi(j) = 0.0d0
 44         continue
c
c      Add appropriately weighted partial derivative values to LLi and RRi.
c
       do 50 r = 1, s
c
        i1 = (r-1)*neqnssq
c
        do 51 l = 1, neqnssq
         LLi(l) = LLi(l) - b(r)*dkdyi(i1+l)
         RRi(l) = RRi(l) - b(r)*dkdyip1(i1+l)
 51     continue
c
 50    continue
c
c      Multiply LLi and RRi by hi
c
       do 53 j = 1, neqnssq
         LLi(j) = hi*LLi(j)
         RRi(j) = hi*RRi(j)
 53    continue
c
c      Add on -I and I, respectively
c
       do 55 j= 1, neqnssq, neqns+1
         LLi(j) = -1.0d0 + LLi(j)
         RRi(j) =  1.0d0 + RRi(j)
 55    continue
c
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mesh_selector(neqns,Nsub,mesh_current,defect,tol,
     *                Nsub_star,mesh_new,info,MxNsub,s_hat)
c
c***overview***
c
c     The purpose of this routine is to examine the relative defect
c     associated with a continuous solution approximation in order to decide
c     upon a good mesh. The `Nsub+1' points describing the current mesh are
c     passed in through the vector `mesh_current'. An estimate of the blended
c     relative/absolute defect, is passed in through the vector `defect',
c     which has `Nsub' components, each of size `neqns', where `neqns' is
c     the number of differential equations. The user defined tolerance, which
c     will be used to assess the relative defect, is passed in through the
c     parameter, `tol'.
c
c       After examining the blended relative/absolute defect, this routine
c     will call either the `half_mesh' routine, to generate a new mesh by
c     halving each subinterval of the current mesh, or the routine
c     `redistribute'. `redistribute' generates a new mesh with a different
c     number of points than the current mesh and with a different distribution
c     of points, depending how the relative defect is distributed over the
c     current mesh.
c
c       The redistribution process involves monitoring the mesh function
c     values which are based on the defect. A mesh redistribution is performed
c     if the maximum mesh function value, is sufficiently greater than
c     the average value. If so, a new mesh is determined that approximately
c     equidistributes the mesh function over the subintervals of the new mesh.
c     The number of points in the new mesh is determined from accuracy
c     considerations.
c
c          The algorithm implemented here is based on [Ascher, Mattheij, and
c     Russell, Chapter 9]. The number of subintervals in the new mesh is
c     returned in `Nsub_star', while the points of the new mesh are returned
c     in the vector `mesh_new'. This routine will fail if the prescribed number
c     of new points is greater than the maximum number allowed, `MxNsub', in
c     which case we set `info = -1' to signal this.
c
c     `s_hat' is used here for temporary storage.
c-------------------------------------------------------------------------------
c***declaration of constants***
         integer Mxs
         parameter (Mxs=10)

c        `Mxs'    is the maximum value for the total number of stages
c                 of the interpolant.
c
c***declaration of factors controlling mesh redistribution***
         double precision     safety_factor
         double precision     rho
         double precision     upper_new_mesh
         double precision     lower_new_mesh
c
         parameter(safety_factor = 1.3d0)
         parameter(rho = 1.0d0)
         parameter(upper_new_mesh = 4.0d0)
         parameter(lower_new_mesh = 0.5d0)
c
c        `safety_factor' a factor to increase the prediction of the appropriate
c                        number of subintervals needed to achieve the given
c                        tolerance in the new mesh.
c
c        `rho'   a factor which determines how much bigger the maximum mesh
c                function value must be, compared to the average value, in
c                order for a mesh redistribution to occur.
c
c                Setting `rho=1' implies that a mesh redistribution will take
c                place everytime.
c
c        `upper_new_mesh', `lower_new_mesh'
c
c                are upper and lower bounds on the ratio between the number
c                of subintervals in current mesh and the number of subintervals
c                in the new mesh.
c
c***declaration of parameters***
c      imports:
         integer              neqns
         integer              Nsub
         double precision     mesh_current(0:Nsub)
         double precision     defect(neqns*Nsub)
         double precision     tol
         integer              MxNsub
c
c        `neqns'        the number of differential equations.
c
c        `Nsub'         the number of subintervals in the current mesh.
c
c        `mesh_current' a vector of points defining the current mesh.
c
c        `defect'       a vector giving an estimate of the blended
c                       relative/absolute defect, with `Nsub' components,
c                       each of size `neqns'.
c
c        `tol'          the user defined tolerance.
c
c        `MxNsub'       is the maximum number of subintervals.
c      exports:
         integer Nsub_star
         double precision mesh_new(0:MxNsub)
         integer info
c
c        `Nsub_star' the number of subintervals of the new mesh.
c
c        `mesh_new' the points defining the new mesh.
c
c        `info' communication flag
c
c                       info = -1    size of the new mesh would be too
c                                      large
c
c                       info =  0    successful generation of a new mesh
c
c***declaration of work array***
         double precision    s_hat(Nsub)
c
c        `s_hat' a vector of mesh function values, with one component per
c                subinterval. The mesh function value is based on the
c                value of the norm of the defect, the tolerance `tol',
c                the size of the subinterval, and the order of the Runge-Kutta
c                scheme, `p'.

c***declaration of local variables***
         integer             i
         double precision    hi
         double precision    r1, r2, r3
         integer             offset
         integer             Nsub_pred
         double precision    norm
c
c        `i'     loop index from 1 to `Nsub'.
c
c        `hi'    size of ith subinterval of current mesh.
c
c                                    1/(p+1)
c                s_hat(i) = (|| defect(i) ||/tol )
c                           ----------------------
c                                    hi
c
c        `r1'    maximum size of weighted mesh function values for any
c                subinterval of current mesh.
c
c        `r2'    weighted sum of mesh function values, `s_hat(i)', with
c                weights, `hi'. `r2' is the predicted number of points
c                that should be in the new mesh.
c
c        `r3'    average of weighted mesh function values.
c
c        `offset'     offset used in the vector `defect'
c                     for the `neqns' solution components of the ith
c                     subinterval.
c
c        `Nsub_pred'  predicted value for number of subintervals in
c                     new mesh if a redistribution is needed.
c
c        `norm'       max norm of the defect estimate for the ith subinterval.
c
c***declaration of variables in common block /interp_order/***
c      imports:
         integer p
         common /interp_order/ p
c
c        `p'      the order of the interpolant associated with the underlying
c                 discrete RK formula.
c
c***declaration of variables in common block /IOCNTRL/***
c      imports:
         integer   print_level_0, print_level_1, print_level_2
         parameter (print_level_0 = 0, print_level_1 = 1,
     *              print_level_2 = 2)
         integer                 profile
         common /IOCNTRL/  profile
c
c        `print_level_0' No output.
c        `print_level_1' Intermediate output.
c        `print_level_2' Full output.
c        `profile' Controls output, to standard output, of profiling infor-
c                  mation such as Newton iteration counts, mesh selection,
c                  relative defect estimate. This variable is identical to
c                  the `output_control' parameter.
c-----------------------------------------------------------------------------
c        called by: `mirkdc'
c        calls to : `half_mesh', `redistribute', `idamax'
                    integer idamax
c------------------------------------------------------------------------------
c     Compute weighted mesh function values as well as r1 and r2.
c     `r1' will be the maximum value of `s_hat(i)*hi'; `r2' will be the sum
c     of these same values. Initialize each to zero. Also set `info' to 0.
c
      info = 0
      r1 = 0.0d0
      r2 = 0.0d0
c
      do 10 i = 1 , Nsub
c
           hi = mesh_current(i)-mesh_current(i-1)
c
c          Compute offset to beginning of defect estimate vector for
c          ith subinterval. This vector is of length `neqns'.
c
           offset = (i-1)*neqns
c
c          Compute norm of the relative defect estimate vector for the
c          ith subinterval.
c
           norm = dabs(defect(offset+idamax(neqns,defect(offset+1),1)))
c
c          Compute corresponding mesh function value, `s_hat(i)'.
c
           s_hat(i) = (norm/tol)**(1.0d0/(p+1))/hi
c
c          Update maximum weighted mesh function value, if necessary.
c
           if ((s_hat(i)*hi) .GT. r1) r1 = s_hat(i)*hi
c
c          Accumulate ith weighted mesh function value in `r2'
c
           r2 = r2 + s_hat(i)*hi
c
 10   continue
c
c     `r3' is the average of the weighted mesh function values
c
      r3 = r2/Nsub
c
c     `r2' gives an estimate of the number of subintervals that the
c     new mesh should have. We employ a safety factor to increase
c     the number of points slightly to increase the probability
c     of achieving the tolerance.
c
      Nsub_pred = (safety_factor*r2)+1
c
c     To avoid cycling, make sure new mesh has a few more points than
c     the previous one.
c
      if ( dabs( (Nsub_pred - Nsub)/(1.0d0*Nsub)) .lt. 0.1d0)
     +                             Nsub_pred = 1.1d0 * Nsub
c
c     Decide whether to half the current mesh or produce a new mesh by
c     redistributing the points of the current mesh, while adding or deleting
c     a few. Redistribute if the maximum weighted mesh function value
c     is more than `rho' times the size of the average value.
c
      if (r1 .LE. rho*r3) then
c
c          The maximum value is not sufficiently large. That is, the
c          distribution of mesh function values is sufficiently uniform
c          and mesh redistribution is unwarrented. In this case, the current
c          mesh is halved.
c
           Nsub_star = 2*Nsub
c
           if (Nsub_star .GT. MxNsub) then
c
c              mesh halving fails; signal this by setting info = -1
               info = -1
c
               if (profile .GT. print_level_0) then
                 write(6,*)' New mesh would be too large'
               endif
c
           else
c
               if (profile .GT. print_level_0) then
                 write(6,*)' Half the current mesh'
               end if
c
               call half_mesh(Nsub,mesh_current,mesh_new)
c              A new mesh twice as fine as the current one is returned
c              in `mesh_new'.
c
           endif
c
      else
c
c          The distribution of mesh function values is sufficiently
c          non-uniform to warrent a mesh redistribution.
c
c          The value of `Nsub_pred' gives the required number of
c          subintervals for the new mesh. However, if this is
c          too big or too small we will not use it.
c
           Nsub_star = Nsub_pred
c
           if ( Nsub_star .GT. upper_new_mesh * Nsub )
     +               Nsub_star = upper_new_mesh * Nsub
           if ( Nsub_star .LT. lower_new_mesh * Nsub )
     +               Nsub_star = lower_new_mesh * Nsub
c
c          Also, check that `Nsub_star' is less than the maximum value.
c
           if (Nsub_star .GT. MxNsub) then
c
c              mesh redistribution fails; signal this by setting info = -1
               info = -1
c
               if (profile .GT. print_level_0) then
                  write(6,*)' New mesh would be too large'
               endif
c
           else
c
c               if (profile .GT. print_level_0) then
c                  write(6,*)'Call to redistibute with Nsub_star=',
c    +                                                Nsub_star
c               end if
c
               call redistribute(Nsub,mesh_current,s_hat,
     +                                       Nsub_star,mesh_new)
c              The new mesh, having `Nsub_star' subintervals, is returned in
c              `mesh_new'
c
           endif
c
      endif
c          `r1. LE. rho*r3'
c
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine NewIter(neqns,leftbc,Nsub,mesh,Y,newtol,maxiter,
     +             info,PHI,delta,top,bot,blocks,pivot,k_discrete,
     +             work,fsub,gsub,dfsub,dgsub)
c
c***overview***
c
c            This routine takes a mesh of `Nsub' subintervals, partitioning
c     the problem interval and a discrete initial approximation provided
c     at the meshpoints, (contained in the vectors `Y' and `mesh',
c     respectively) and uses a modified Newton method to produce a
c     solution approximation to within a tolerance of `newtol', for
c     the associated discrete nonlinear system `PHI(Y) = 0', where `PHI(Y)' is
c     the residual function. The number of differential equations
c     is `neqns', while the number of boundary conditions at the lefthand
c     endpoint of the problem interval is `leftbc'. The iteration completes
c     successfully if it produces a converged solution within `maxiter'
c     iterations. A successful return is indicated by `info'= 0. The
c     updated solution is returned in `Y'. The arrays `PHI', `delta', `top',
c     `bot', `blocks', `pivot', and `work' provide working space for the
c     `resid', `damped_newt', and `fixed_jacob' routines. This routine
c     provides, as a by-product, the stage values associated with the
c     application of a Runge-Kutta method on each subinterval of the mesh.
c     The array `k_discrete' returns these stages. The subroutines,
c     `fsub', `gsub', `dfsub', and `dgsub', are provided by the user to
c     define the boundary value ODE and the boundary conditions and their
c     derivatives.
c
c            The modified Newton iteration implemented here proceeds
c     through a sequence of steps, each of which must be either a damped
c     Newton step or a fixed Jacobian step. The code decides, depending
c     on how well the iterates appear to be converging, which type of step
c     to take each time. This routine is organized as follows. After some
c     initialization, we begin a loop which executes the modified Newton
c     iteration. The kind of step to be taken is controlled by the value of
c     the flag, `Fixed_Jacobian', which is initially set to FALSE since the
c     first step must be a full Newton step.
c
c            During a damped Newton step an iterative process is used to
c     choose a damping factor, `lambda', which is used to provide a damped
c     update to the current iterate `Y = Y - lambda*delta' where `delta' is
c     the Newton correction; `delta = inv(J)*PHI' where `PHI' is the residual
c     function evaluated at `Y' and J is the corresponding Jacobian also
c     evaluated at `Y'. The use of a damping factor can improve the
c     convergence properties of the Newton iteration as well as allowing for
c     the early detection of a divergent iteration. The update is accepted
c     if it produces a sufficient reduction in the value of the natural
c     criterion function (which is related to the size of the residual
c     function). At the end of a successful damped Newton step, the routine
c     takes a fixed Jacobian step if the step was a full Newton step
c     (i.e. `lambda' = 1).
c
c          The fixed Jacobian step consists of first computing a potential
c     new iterate, `Y_hat = Y - delta' where `delta' is the Newton correction
c     for this fixed Jacobian step. Here `delta = inv(J)*PHI', where `PHI' is
c     the residual function evaluated at `Y' and where J, the Jacobian, has
c     not been evaluated at `Y', but rather at some previous iterate.
c     The new iterate, `Y_hat' will be acceptable if it leads to a sufficient
c     reduction in the size of the natural criterion function. If the new
c     iterate value is accepted then another fixed Jacobian step is attempted.
c     Otherwise, the code returns to taking full Newton steps.
c
c            This routine can fail in three ways. If the Newton iteration takes
c     too many steps, the iteration halts with `info'=1. If a singular matrix is
c     encountered during the factorization of the Newton matrix (during the call
c     to `damped_newt'), the iteration halts with `info'=2. If it is impossible
c     for the `damped_newt' routine to obtain a suitable damping factor for the
c     Newton update (indicative of an ``effectively singular" Jacobian) the
c     iteration halts with `info'=3.
c-----------------------------------------------------------------------------
c***declaration of constants***
      integer  Mxs
      parameter(Mxs=10)
c
c     `Mxs' Maximum number of stages of the Runge-Kutta method.
c
c***declaration of parameters***
c   imports:
      integer                 neqns, leftbc, Nsub
      double precision        mesh(0:Nsub)
      double precision        Y(neqns*(Nsub+1))
      double precision        newtol
      integer                 maxiter
c
c     `neqns' Number of differential equations.
c     `leftbc' Number of boundary conditions at the left
c              end of the problem interval.
c     `Nsub' Number of subintervals into which the problem
c            interval is partitioned.
c     `mesh' Set of points which partition the problem interval.
c      `Y' On input, the initial approximation for the Newton iteration.
c     `newtol' Tolerance value to be applied to the Newton iteration,
c              i.e. the iteration stops if the Newton correction
c              is less than or equal to the user specified
c              tolerance. A blended relative/absolute tolerance
c              is applied. If `delta_j' is the Newton correction
c              for, `Y_j', the jth component of `Y', then we require
c                     `(|delta_j| / (|Y_j|+1)) <= newtol'.
c     `maxiter' Maximum number of damped Newton iterations
c               to be allowed during each call to this routine.
c   exports:
c     double precision      Y(neqns*(Nsub+1))
      integer               info
      double precision      k_discrete(Mxs*neqns*Nsub)
c
c     `Y' On output, when `info' = 0, the converged solution to
c         the Newton iteration.
c     `info' Internal communication flag:
c            0 - successful iteration converged to
c                within user specified tolerance.
c            1 - unsuccessful termination - more than `maxiter' iterations
c                were taken without convergence being achieved.
c            2 - unsuccessful termination - a singular coefficient matrix
c                was encountered during the attempted solution of the
c                Newton system.
c            3 - unsuccessful termination - it was impossible to obtain a
c                suitable damping factor for the Newton update (indicative
c                of an ``effectively singular" Jacobian) or evaluation
c                of natural criterion function overflowed.
c     `k_discrete' The discrete stages for all subintervals. The i-th set of
c                  `s*neqns' locations of `k_discrete' contains the `s
c                  stages, each of length `neqns' corresponding to the i-th
c                  subinterval.
c   work space:
      double precision      PHI(neqns*(Nsub+1))
      double precision      delta(neqns*(Nsub+1))
      double precision      top(neqns*neqns)
      double precision      bot(neqns*neqns)
      double precision      blocks(2*neqns**2 * Nsub)
      integer               pivot(neqns*(Nsub+1))
      double precision
     +            work(3*neqns*(Nsub+1)+neqns+2*neqns**2*(Mxs+1))
c
c     `PHI' Residual function; it has one set of `neqns' locations per
c          subinterval. The first `leftbc' locations correspond to the
c          left boundary conditions; the final `neqns-leftbc' locations
c          correspond to the right boundary conditions. Used as the
c          right hand side for the Newton systems.
c     `delta' Newton correction vector, equal to `inv(J)*PHI(Y)',
c            where J is the Jacobian evaluated at the current iterate,
c            `Y', when a damped Newton step is taken.
c     `top' Storage of derivative information for left boundary
c          conditions.
c     `bot' Storage of derivative information for right boundary
c          conditions.
c     `blocks' Storage of derivative information for residual function.
c     `pivot' Pivoting information associated with the factorization
c              of J. Together `top', `blocks', and `bot' contain
c              the almost block diagonal matrix representing the
c              coefficient matrix of the Newton system.
c     `work' Work array passed to `resid' and `fixed_jacob'.
c
c***user-supplied subroutines***
      external fsub,gsub,dfsub,dgsub
c
c     `fsub' Defines f(t,y) for the first order system
c            of differential equations, y' = f(t,y).
c     `gsub' Defines the boundary condition equations.
c     `dfsub' Defines the Jacobian of the system of
c             differential equations.
c     `dgsub' Defines the Jacobian of the boundary conditions.
c
c***declaration of local variables***
c
      logical            convrg
      integer            count
      double precision   delta_0_norm, lambda, g
      logical            Fixed_Jacobian
      integer            i_delta_0, i_work_d
      integer            i_Y_hat, i_PHI_hat, i_work_f
c
c     convrg' Used to check for convergence of the Newton iteration.
c     count'  Counts number of damped Newton iterations.
c     delta_0_norm' Norm of the Newton correction vector `delta_0'.
c     lambda' Damping factor used to perform a damped Newton update.
c     g' Value of natural criterion function.
c     Fixed_Jacobian' .TRUE. if the next step to be taken should be
c                     a fixed Jacobian step. .FALSE. if the next step
c                     to be taken should be a damped Newton step.
c     i_delta_0' Index into `work' for the start of the part of
c                `work' that will correspond to the temporary
c                storage array called `delta_0' within `damped_newt'.
c     i_work_d' Index to `work' for the start of the part of
c               `work' that will correspond to the work space used
c                within `damped_newt'.
c     i_Y_hat' Index to `work' for the start of the part of
c              `work' that will correspond to the temporary storage
c              array `Y_hat' used within `fixed_jacob'.
c     i_PHI_hat' Index to `work' for the start of the part of
c                `work' that will correspond to the temporary storage
c                array `PHI_hat' used within `fixed_jacob'.
c     i_work_f' Index to `work' for the start of the part of
c               `work' that will correspond to the work space used
c                within `fixed_jacob'.
c
c***declaration of variables for /IOCNTRL/ common block***
c   imports:
      integer   print_level_0, print_level_1, print_level_2
      parameter (print_level_0 = 0, print_level_1 = 1,
     *              print_level_2 = 2)
      integer                 profile
      common /IOCNTRL/        profile
c
c     `print_level_0' No output.
c     `print_level_1' Intermediate output.
c     `print_level_2' Full output.
c     `profile' Controls output, to standard output, of profiling infor-
c               mation such as Newton iteration counts, mesh selection,
c               and defect estimates.
c-----------------------------------------------------------------------------
c      Called by: `mirkdc'
c      Calls to: `fixed_jacob',`damped_newt', `resid'
c-----------------------------------------------------------------------------
c***initialization of work array indexes
       i_delta_0 = 1
       i_work_d    = neqns*(Nsub+1) + 1
       i_Y_hat = 1
       i_PHI_hat    = neqns*(Nsub+1) + 1
       i_work_f    = 2*neqns*(Nsub+1) + 1
c
c***initialization for Newton iteration***
c      The number of damped Newton iterations (equal to the number of
c      Jacobian evaluations and factorizations) will be counted using the
c      variable `count' which is initialized to 0. The variable `convrg',
c      monitors whether the iteration has converged or not and is
c      initialized to .FALSE..
c
       count = 0
       convrg = .FALSE.
c
c      The first iteration step must be a damped Newton step.
c      This is handled by setting `Fixed_Jacobian' to .FALSE.
c
       Fixed_Jacobian = .FALSE.
c
c      When a damped Newton step is used, it is assumed that the current
c      iterate, `Y', and also the value of, `PHI', the residual function
c      at `Y' will be available from the previous iteration step. For the
c      first step, which will be a damped Newton step, `Y' will be
c      available, but the residual function must be computed here in
c      preparation for the first time through the damped Newton iteration
c      step. This is done by calling `resid' which will return with the
c      residual function stored in the vector `PHI'.
c
       call resid(neqns, leftbc, Nsub, mesh, Y, PHI,
     +                                   k_discrete, work, fsub, gsub)
c
c      In order to predict the new damping factor, each damped Newton step
c      also assumes several other quantities are available from the previous
c      iteration step. For the first step, these will not be available. In
c      order for the prediction algorithm to detect this situation we set
c      `lambda=0' (The value of `lambda' can never reach 0 during the
c      usual execution of the algorithm).
c
       lambda = 0.0d0
c
c***Newton iteration loop***
c   REPEAT-UNTIL LOOP (Implemented here using an if-goto pair)
 1                                                        continue
       if (Fixed_Jacobian) then
c
              if (profile .GT. print_level_1) then
                write(6,*)'Fixed Jacobian step.'
              end if
c
              call fixed_jacob(neqns,leftbc,Nsub,mesh,Y,newtol,
     +             delta,g,PHI,top,bot,blocks,pivot,Fixed_Jacobian,
     +             lambda, convrg,info, k_discrete,work(i_Y_hat),
     +             work(i_PHI_hat),work(i_work_f), fsub,gsub)
c
       else
c
c             Update `count' to monitor the number of damped Newton iterations
              count = count + 1
c
              if (profile .GT. print_level_1) then
                write(6,100)'Damped Newton step, count =',count,'.'
 100            format(1x,a27,i3,a1)
              end if
c
              call damped_newt(neqns,leftbc,Nsub,mesh,Y,newtol,
     +               lambda,PHI,top,bot,blocks,pivot,Fixed_Jacobian,
     +               convrg,delta,delta_0_norm,g,info,k_discrete,
     +               work(i_delta_0),work(i_work_d),fsub,gsub,
     +               dfsub,dgsub)
       endif
c      End of `If (Fixed Jacobian)'
c
       if ((.NOT. convrg) .AND. (info .LE. 0)) then
c                Check to see that the maximum number of iterations has
c                not been exceeded.
                 if ( count .GT. maxiter) info = 1
       endif
c
c      Repeat the loop unless convergence has been obtained, the allowed
c      number of iterations has been exceeded, a singular matrix was
c      encountered during the Newton iteration, or a suitable damping
c      factor has not been found (the latter three signaled by `info' > 0).
c
c   UNTIL ( convrg .OR. (info .GT. 0))
                                    if ((.NOT. convrg) .AND.
     +                              (info.LE. 0)) go to 1
c***end of Newton iteration loop***
c
c***conclusion***
      if (profile .GT. print_level_0) then
           if (count .GT. 1) then
             write(6,200)'After',count,' Newton steps'
           else
             write(6,300)'After',count,' Newton step'
  200        format(1x,a5,i3,a13)
  300        format(1x,a5,i3,a12)
           endif
      end if
c
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine newmat(neqns,leftbc,Nsub,mesh,Y,top,blocks,bot,
     +                                      k_discrete,work,dfsub,dgsub)
c
c***overview***
c
c            This routine computes the Newton matrix for the discrete system
c     modeling the BVP.  The Newton matrix is the Jacobian of the residual
c     function.  It is evaluated at the current iterate `Y'. The Jacobian is
c     returned in the three vectors `top',`bot', and `blocks'; `top' contains
c     the elements of the Jacobian associated with the boundary  conditions
c     evaluated at the left endpoint - it consists of `leftbc' rows, each
c     of size `neqns', where `leftbc' is the number of boundary conditions
c     imposed at the left endpoint, and `neqns' is the number of differential
c     equations; `bot'  contains  the  corresponding information for the right
c     endpoint. It consists of `neqns-leftbc' rows, each of length `neqns'.
c     The vector `blocks' contains, for the i-th subinterval of the mesh which
c     partitions the problem interval, as prescribed by the vector, `mesh', the
c     blocks L_i and R_i, i=1,...,Nsub where `Nsub' is the number of subinter-
c     vals defined by the mesh.
c
c           L_i is the derivative of the i-th component of the residual function
c     with respect to y_i; R_i is the derivative of the i-th component of the
c     residual function with respect to y_i+1. These blocks, stored column
c     by column, are placed in `blocks'  one  after the  other  in  the order
c     L_i, R_i, i=1,...,Nsub. Each such block is `neqns' rows, each of length
c     `neqns'.
c
c            This routine accesses the discrete stages for each subinterval,
c     which are computed by the `subcom' routine, and saved in the array
c     `k_discrete' provided to this routine through the parameter list. Thus
c     a precondition for the use of this routine is that the `resid' routine
c     has already been called, with the same `Y' input.  `work' is a work
c     array provided to `jacblk' and also used here.
c-------------------------------------------------------------------------------
c***declaration of constants***
      integer Mxs
      parameter(Mxs=10)
c
c     `Mxs' Maximum value for the number of stages of the RK method.
c
c***declaration of parameters***
c   imports:
      integer                 neqns,leftbc,Nsub
      double precision        mesh(0:Nsub), Y(neqns*(Nsub+1))
      double precision        k_discrete(Mxs*neqns*Nsub)
c
c     `neqns' Number of differential equations.
c     `leftbc' Number of boundary conditions at the left end of the problem
c                                                                     interval.
c     `Nsub' Number of subintervals into which the problem
c                                                 interval is partitioned.
c     `mesh' Set of points which partition the problem interval.
c     `Y' Current discrete approximation to the solution,
c                                          evaluated at each meshpoint.
c     `k_discrete' Vector containing the discrete stages for
c                  all subintervals. The i-th set of `s*neqns'
c                  locations of `k_discrete' contains the `s'
c                  stages, each of length `neqns' corresponding
c                  to the i-th subinterval.
c   exports:
      double precision      top(leftbc*neqns)
      double precision      blocks(2*neqns**2*(Nsub+1))
      double precision      bot((neqns-leftbc)*neqns)
c
c     `top' Derivative information for left boundary conditions.
c     `blocks' Derivative information for residual function.
c     `bot' Derivative information for right boundary conditions.
c
c***declaration of work arrays***
      double precision         work(neqns+2*(Mxs+1)*neqns**2)
c
c     `work' Passed for working storage to the `jacblk' routine and
c            used here as the vector where the values of the
c            derivatives of the boundary conditions are returned from
c            the user defined routine `dgsub'.
c
c***user-supplied subroutines***
      external dfsub, dgsub
c
c     `dfsub' Defines the Jacobian of the system of differential equations.
c     `dgsub' Defines the Jacobian of the boundary condition equations.
c
c***declaration of local variables***
      integer                  i, j, k, neqnsq
      integer                  i1,i2,i3,i4,n1
      integer                  toploc, botloc, workloc
      double precision         h, t
c
c     `i,j,k' Loop variables; `i' ranges from 1 to `Nsub'; `j,k' range from 1
c             to `neqns'.
c     `i1',`i2',`i3',`i4',`n1' Indexes for the variables `Y', `k_discrete', and
c                              `blocks' .
c     `toploc', `botloc', `workloc' Offset locations in the vectors `top',
c                                   `bot', and `work'.
c     `neqnsq' Equals `neqns*neqns'.
c
c     `h' Size of the i-th subinterval, `mesh(i) - mesh(i-1)'.
c     `t' Value of the (i-1)st mesh point, `mesh(i-1)'.
c
c***declarations of indices for work array passing to `jacblk'***
      integer        i_yr, i_dfdy, i_ur, i_dkdyi, i_dkdyip1

c     `i_yr' Starting location in `work': work(1), for `y_r', the argument for
c            the evaluation of the `rth' stage.
c     `i_dfdy' Starting location in `work' :  work(neqns+1), for `dfdy', which
c              during the `rth' time through the main loop below, contains the
c              derivative of the right hand side of the system of ODEs,
c              evaluated at the arguments of the `rth' stage of the RK scheme.
c     `i_ur' Starting location in `work' : work(neqns+neqns^2+1), for `ur',
c            which is used to store the intermediate values arising in the
c            computation of each partial derivative.
c     `i_dkdyi' Starting location in `work' : work(neqns+2*neqns^2+1), for
c               `dkdyi'; the `rth' set of `neqns**2' locations of this vector
c               will contain the value of the `rth' partial derivative
c               with respect to `y_i'.
c     `i_dkdyip1' Starting location in `work'  : work(neqns+2*neqns^2+
c                                                              Mxs*neqns^2+1),
c                 for `dkdyip1'; the `rth' set of `neqns**2' locations of this
c                 vector will contain the value of the `rth' partial derivative
c                 with respect to `y_ip1'.
c
c***declaration of variables in common block /RK_s/***
      integer                 s
      common /RK_s/ s
c
c     `s' Number of discrete stages of the Runge-Kutta formula.
c----------------------------------------------------------------------------
c      Called by: `damped_newt'.
c      Calls to: `jacblk', `dgsub'.
c-----------------------------------------------------------------------------
c***computation of `blocks'***
c
      neqnsq     = neqns**2
      i_yr       =  1
      i_dfdy     =  neqns+1
      i_ur       =  neqns+neqnsq+1
      i_dkdyi    =  neqns+2*neqnsq+1
      i_dkdyip1  =  neqns+2*neqnsq+Mxs*neqnsq+1
c
c     We need to compute the Jacobian blocks L and R for each subinterval.
c
c     For each subinterval ...
c
      i1 = s*neqns
      n1= Nsub*neqns
c
      do 5 i=1, Nsub
           i2 = (i-1)*neqns
           i3 = (i-1)*i1
           i4 = 2*(i-1)*neqnsq
           h = mesh(i) - mesh(i-1)
           t = mesh(i-1)
c
c          y_i and y_i+1, the current solution approximations at the left and
c          right endpoints of the i-th subinterval, respectively, are required
c          for the computation of the i-th Jacobian blocks. Y((i-1)*neqns+1)
c          is y_i; Y(i*neqns+1) is y_i+1. We will also need the `s' stage
c          evaluations associated with the `ith' subinterval. These are already
c          available in the `ith' set of `s*neqns' locations of `k_discrete',
c          `k_discrete((i-1)*s*neqns+1)'.
c
c          The Jacobian blocks are passed to the appropriate places
c          in `blocks.' L is stored in the neqnsq locations starting at
c          location 2*(i-1)*neqnsq+1; R is stored in the neqnsq locations
c          (2*i-1)*neqnsq+1.
c
           call jacblk(neqns,h,t,Y(i2+1),Y(i2+neqns+1),
     +          k_discrete(i3+1), blocks(i4+1), blocks(i4+neqnsq+1),
     +          work(i_yr),work(i_dfdy),work(i_ur),
     +          work(i_dkdyi), work(i_dkdyip1),dfsub)
c
 5      continue
c
c***computation of `top' and `bot'***
c
c       We now compute the parts of the Jacobian of the residual
c       function associated with the boundary conditions.
c
c       `Y(1)' is y_0, `Y(Nsub*neqns+1)' is y_Nsub. We call `dgsub' to get
c       the boundary condition derivative information. The information is
c       returned in the vector `work'. Although, `work' is viewed inside
c       `dgsub' as a square matrix of size `neqns' by `neqns', it is more
c       convenient to view it here in the usual columnwise fashion, i.e.
c       with the columns of the square matrix laid out one after the
c       other in the column vector, `work'.
c
c       Set the Jacobian to zero then call `dgsub' to obtain the
c       non-zero entries.
c
        do 7 i = 1,neqns*neqns
          work(i) = 0.0d0
 7      continue
c
        call dgsub(neqns,Y(1),Y(n1+1),work)
c
c       The vector `top' must be assigned the values corresponding to
c       the first `leftbc' rows of the square version of the boundary condition
c       matrix. Since `work' contains the elements of the boundary condition
c       matrix stored in a columnwise fashion, this corresponds to selecting
c       the first `leftbc' elements from each column stored within `work', and
c       copying them into the appropriate locations of the vector `top'.
c
c       The vector `top' will also store the first `leftbc' rows of the
c       boundary condition Jacobian in a columnwise fashion. Thus it will
c       consist of `neqns' vectors, each of length `leftbc', representing the
c       columns of the `leftbc' by `neqns' matrix corresponding to the
c       derivatives of the boundary conditions at the left end of the problem
c       interval.
c
c       A similar situation holds for the vector `bot', which will
c       contain, in a columnwise fashion, the `neqns-leftbc' by `neqns'
c       matrix, corresponding to the derivatives of the boundary conditions
c       at the right end of the problem interval. The variables `toploc',
c       `workloc', and `botloc' are used to identify, within their respective
c       vectors, the appropriate locations for the jth column, for j=1 to
c       `neqns'.
c
        do 25 j = 1, neqns
c
           toploc = (j-1)*leftbc
           workloc = (j-1)*neqns
           botloc = (j-1)*(neqns-leftbc)
c
c          Copy the first leftbc locations of the jth "column" of `work'
c          into the leftbc locations corresponding to the jth "column" of `top'.
c
           do 10 k = 1, leftbc
            top(toploc+k) = work(workloc+k)
10         continue
c
c          Copy the remaining `neqns-leftbc' locations of the jth "column"
c          of `work' into the `neqns-leftbc' locations corresponding to the
c          jth "column" of `bot'.
c
           do 20 k = leftbc+1, neqns
            bot(botloc+k-leftbc) = work(workloc+k)
20         continue
c
25      continue
c
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine redistribute(Nsub,mesh_current,s_hat,Nsub_star,
     +                         mesh_new)
c
c***overview***
c
c         This routine takes as input a mesh of `Nsub' subintervals which
c      partitions the problem interval, stored in the vector `mesh_current',
c      and, for each subinterval of this mesh, the value of the mesh function,
c      stored in the vector, `s_hat'. The mesh function, which is piecewise
c      constant, i.e. has a constant value on each subinterval of the current
c      mesh, associated with some measure of the relative defect on each
c      subinterval (see `mesh_selector' routine for further details). The
c      routine returns a new mesh, of `Nsub_star' points, which equidistributes
c      the integral of the mesh function over the subintervals of the new mesh
c      in the vector `mesh_new'.
c
c         The equidistribution is performed as follows. The idea is to select
c      new mesh points so that the area under the mesh function is equally
c      distributed over the subintervals of the new mesh. The total area under
c      the mesh function is first computed and then divided by the total number
c      of subintervals in the new mesh, giving a quantity we call `zeta'.
c      The value of the integral of the mesh function over each subinterval of
c      the new mesh must equal `zeta'. The area under the mesh function for
c      each subinterval of the old mesh is computed and added to a sum,
c      proceeding from left to right. Whenever the sum reaches a value beyond
c      `zeta', only the area of that fraction of the subinterval needed to
c      reach the value of `zeta' is added. This fraction is used to decide
c      where to place the new mesh point. The procedure beginning with the new
c      mesh point is repeated and continued until every subinterval of the
c      current mesh is processed. See [Ascher, Mattheij, Russell, Chap.9] for
c      details
c
c         This routine assumes that `Nsub_star' is less than or equal to the
c      maximum allowable number of subintervals.
c--------------------------------------------------------------------------------
c***declaration of parameters***
c      imports:
          integer                Nsub
          double precision       mesh_current(0:Nsub)
          double precision       s_hat(Nsub)
          integer                Nsub_star
c
c         `Nsub'         the number of subintervals in the input mesh
c
c         `mesh_current' the vector containing the current mesh
c
c         `s_hat'        the vector containing the mesh function
c                        value for each subinterval.
c
c         `Nsub_star'    the number of points in the new mesh.
c      exports:
          double precision      mesh_new(0:Nsub_star)
c
c         `mesh_new'   the meshpoints making up the new mesh
c
c***declaration of local variables***
          integer i,k
          double precision sum
          double precision zeta
          double precision t
          double precision Integral
          double precision Next_Piece
c
c         `i'          loop index from 1 to Nsub_star.
c
c         `k'          loop index from 1 to Nsub.
c
c         `sum'        temporary variable used to compute `zeta'.
c
c         `zeta'       weighted sum of mesh function values, averaged over
c                      number of subintervals in the new mesh, `Nsub_star'.
c
c         `t'          used as marker to indicate the point within the
c                      problem interval to which the integration for the
c                      determination of new mesh points has proceeded so far.
c
c         `Integral'   used to store accumulation of contributions
c                      to the area under the mesh function for each
c                      subinterval of the new mesh.
c
c         `Next_Piece' used to store the value of each contribution
c                      to `Integral' as it arises during the determination
c                      of the position of each new mesh point.
c------------------------------------------------------------------------------
c         called by: `mesh_selector'
c------------------------------------------------------------------------------
c      Compute equidistribution constant, `zeta'. It equals the
c      weighted sum of the mesh function values divided by the desired
c      number of mesh points in the new mesh, `Nsub_star'.
c
       sum = 0.0d0
c
       do 10 k = 1, Nsub
            sum = sum + s_hat(k)*(mesh_current(k)-mesh_current(k-1))
 10    continue
c
       zeta = sum/Nsub_star
c
c      Loop to determine points of new mesh.
c
c      k gives index of the subinterval of the current mesh whose mesh
c      function contribution is being added to the integral for the ith
c      new mesh point.
c
       k = 1
       i = 0
       mesh_new(0) = mesh_current(0)
       t = mesh_current(0)
       Integral = 0.0d0
c
c      REPEAT
c                                                 --------
 1                                                continue
c                                                 --------
c
           Next_Piece = s_hat(k)*(mesh_current(k)-t)
c
           if ((Integral+Next_Piece) .GT. zeta) then
c
c              The next piece is too big. Set the next new meshpoint
c              so that it allows us to include only the appropriate
c              fraction of the current piece necessary to make the
c              integral for this new subinterval equal to `zeta'.
c
               mesh_new(i+1) = (zeta - Integral)/s_hat(k) + t
c
c              Set the "point of integration" marker, `t', to this new mesh
c              point value, set Integral to zero, and increment `i' in
c              preparation for the determination of the next new mesh point.
c
               t = mesh_new(i+1)
               i = i + 1
               Integral = 0
c
           else
c
c              The next piece to be added to the integral is still not
c              large enough to make the value of the integral `zeta'
c              or bigger, so just add it on to the current integral
c              value and adjust the "point of integration" marker, `t'.
c
               Integral = Integral + Next_Piece
               t = mesh_current(k)
c
c              The contribution from the kth subinterval of the current
c              mesh is complete. Increment k.
c
               k = k + 1
c
           endif
c
c      UNTIL ( k .GT. Nsub)
c                                         ----------------------
                                          if (k. LE. Nsub) goto1
c                                         ----------------------
c
c      Finish up the process by setting the last new mesh point value
c      to the value of the old last mesh point.
c
       mesh_new(Nsub_star) = mesh_current(Nsub)
c
       return
       end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine resid(neqns,leftbc,Nsub,mesh,Y,PHI,
     +                                    k_discrete,work,fsub,gsub)
c
c***overview***
c
c          This routine evaluates the residual function at the current
c     iterate, `Y', on the current mesh, `mesh'. The residual function value
c     is  returned  in  the vector `PHI'. The mesh has `Nsub' subintervals.
c
c          The evaluation  of the residual function is done in two phases:
c     one involving the calculation of the components  of `PHI' associated
c     with the boundary conditions (`leftbc' conditions at the left boundary
c     and `neqns-leftbc' at the right boundary, where `neqns' is the number of
c     differential equations), and the other associated  with  the
c     other `Nsub' components, each of size `neqns', corresponding  to each
c     subinterval of the current mesh. The  components  associated with
c     the boundary  conditions  are  obtained  directly  from  the  subroutine
c     `gsub'.  Each  of  the other components  is obtained through a call
c     to the `subcom' routine.
c
c          In  the  first  part of  this  routine we perform the  necessary
c     initializations.  We next treat the components of `PHI'  related  to the
c     subintervals, and lastly deal with the components  corresponding  to the
c     boundary conditions.
c
c          A by-product of the call to `subcom' routine is the return
c     of the `s' stages that arise as intermediate quantities during the
c     computation of the part of `PHI' associated with the i-th subinterval.
c     These stage values will be used later during the call to `Newmat' and the
c     defect estimation phases of this code, and are stored here in the
c     `k_discrete' vector which is returned from this routine through the
c     parameter list. `work' is a work array provided to the `subcom' routine.
c     The subroutines, `fsub' and `gsub' are provided by the user to define the
c     boundary value ODE and the boundary conditions.
c-------------------------------------------------------------------------------
c***declaration of constants***
      integer Mxs
      parameter (Mxs=10)
c
c     `Mxs'    is the maximum number of stages.
c
c***declaration of parameters***
c   imports:
      integer                 neqns,leftbc,Nsub
      double precision        mesh(0:Nsub), Y(neqns*(Nsub+1))
c
c     `neqns' Number of differential equations.
c     `leftbc' Number of boundary conditions at the left end of the problem
c                                                                    interval.
c     `Nsub' Number of subintervals into which the problem interval is
c                                                               partitioned.
c     `mesh' Set of points which partition the problem interval.
c     `Y' Current discrete approximation to the solution,
c                                        evaluated at each meshpoint.
c
c   exports:
      double precision     PHI(neqns*(Nsub+1))
      double precision     k_discrete(Mxs*neqns*Nsub)
c
c     `PHI' Value of the residual function on each subinterval.
c           It has one set of `neqns' locations per subinterval.
c           In addition, the first `leftbc' locations correspond to the
c           left boundary conditions; the final `neqns-leftbc' locations
c           correspond to the right boundary conditions.
c     `k_discrete' Vector containing the discrete stages for
c                  all subintervals. The i-th set of `s*neqns'
c                  locations of `k_discrete' contains the `s'
c                  stages, each of length `neqns' corresponding
c                  to the i-th subinterval.
c   work arrays:
      double precision      work(neqns)
c
c     `work' Temporary work array required by the `subcom' routine
c            and also used here to receive the values of the residual
c            function corresponding to the boundary conditions as
c            returned from the user defined routine `gsub'.
c
c***declaration of variables in common block /RK_s/***
      integer                 s
      common /RK_s/ s
c
c     `s' Number of discrete stages of the Runge-Kutta formula.
c
c***user-supplied subroutines***
      external fsub, gsub
c
c     `fsub' Defines f(t,y) for the first order
c                      system of differential equations, y' = f(t,y).
c     `gsub' Defines the boundary condition equations.
c
c***declaration of local variables***
      integer                 i, j, i1, n1
      double precision        h, t
c
c     `i,j'  Loop variables; `i' ranges from 1 to `Nsub';
c            `j' ranges from 1 to `neqns'.
c     `i1','n1' Indexes used for the variables `Y', `PHI', and `k_discrete'.
c     `h' Size of the i-th subinterval, `mesh(i) - mesh(i-1)'.
c     `t' Value of the i-th mesh point, `mesh(i-1)'.
c-------------------------------------------------------------------------------
c     Called by: `newiter','criterion', `fixed_jacob'.
c     Calls: `subcom', `gsub'.
c-------------------------------------------------------------------------------
c***computation of the main part of `PHI' associated with the subintervals***
c
c     For each subinterval ...
c
      do 5 i=1, Nsub
           i1 = (i-1)*neqns
           h = mesh(i) - mesh(i-1)
           t = mesh(i-1)
c
c          y_i and y_ip1, the current solution approximations at the endpoints
c          of the i-th subinterval, are required for the computation of the i-th
c          residual component. Y((i-1)*neqns+1) is y_i, Y(i*neqns+1) is y_ip1.
c
c          We now compute the component of `PHI' associated with the i-th sub-
c          interval. Note that `PHI(leftbc+(i-1)*neqns+1)' is simply
c          the place in `PHI' where this information is stored. The stages
c          of the discrete formula, for the i-th subinterval, are returned
c          by the `subcom' routine to the appropriate part of the `k_discrete'
c          vector.
c
           call subcom(neqns,h,t,Y(i1+1),Y(i1+neqns+1),
     +          PHI(leftbc+i1+1), k_discrete(i1*s+1),work,fsub)
c
 5    continue
c
c***boundary condition case***
c
      n1=Nsub*neqns
c
c     Y(1) is y_0, Y(Nsub*neqns+1) is y_Nsub. We call `gsub' to get the
c     boundary condition information.
c
      call gsub(neqns,Y(1),Y(n1+1),work)
c
c     We next assign the boundary condition info, stored in the array
c     `work' to the appropriate part of `PHI'.
c
      do 10 j = 1, leftbc
           PHI(j) = work(j)
 10   continue
c
      do 15 j = leftbc+1, neqns
           PHI(n1 + j) = work(j)
 15   continue
c
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         subroutine RK_tableau (method)
c
c***overview****
c
c          This routine sets up all coefficients related to the discrete
c      RK formula. `method' defines the MIRK method to be used. The possible
c      values for `method' and the corresponding MIRK scheme associated with
c      each of these values are listed below.
c
c             Value of `method' |   MIRK formula
c            -------------------------------------
c                  121          |   MIRK121 scheme
c                  221          |   MIRK221 scheme
c                  232          |   MIRK232 scheme (abscissa = 0,2/3)
c                  2321         |   MIRK232 scheme (abscissa = 1,1/3)
c                  343          |   MIRK343 scheme
c                  453          |   MIRK453 scheme
c                  563          |   MIRK563 scheme
c                  5631         |   An improved MIRK563 scheme
c
c      This routine assigns values to the parameters `v',`b',`c',`x' and `s';
c      thus defining the MIRK formula to be used. The tableau
c      of the method is
c
c            c(1) | v(1) |   0        0      ...     0
c            c(2) | v(2) | x(2,1)     0      ...     0
c            c(3) | v(3) | x(3,1)  x(3,2)    ...     0
c              .       .      .      .           .
c              .       .      .      .         .
c              .       .      .      .           .
c            c(s) | v(s) | x(s,1)  x(s,2)    ...     0
c            --------------------------------------------------------
c                          b(1)    b(2)      ...    b(s)
c
c          These coefficients are then available through the
c          /RK_s/,/RK_coeff/ common blocks.
c------------------------------------------------------------------------------
c***declaration of constants***
            integer Mxs
            parameter (Mxs=10)
c
c           `Mxs'    maximum number of stages in the discrete formula.

c***declaration of parameters***
c        imports:
            integer        method
c
c           `method'  defines the method to be used
c
c***declaration of variables in common block /RK_s/***
            integer                 s
            common /RK_s/ s
c
c           `s' is the number of discrete stages of the Runge-Kutta formula.

c***declaration of variables in common block /RK_coeff/***
            double precision        v(Mxs), b(Mxs)
            double precision        c(Mxs), x(Mxs**2)
            common /RK_coeff/ v, b, c, x
c
c           `v', `b', `c',and `x' are coefficients that define
c           the discrete Runge-Kutta formula.  See tableau above.
c
c           The first component of `c' contains the abscissa for the first
c           new stage. The second component of `c' contains the abscissa for
c           the second new stage, and so on.
c
c           The first component of `v' contains the blending coefficient
c           for the first new stage. The second component of `v' contains
c           the blending coefficient for the second new stage, and so on.
c
c           The matrix, X, is related to  the data structure, `x'
c           as follows (the matrix is stored columnwise in the vector):
c
c                 X(i,j) = x( (j-1)*s + i )
c
c           `x' contains the coupling coefficients for each of the new
c           stages
c-----------------------------------------------------------------------------
c     called by: `mirkdc'
c------------------------------------------------------------------------------
c       Translate short form of method values.
        if (method .EQ. 2) method = 221
        if (method .EQ. 4) method = 343
        if (method .EQ. 6) method = 563

c      Beginning method definition based on value of `method'.

       if (method .EQ. 121) then
c
c          Define the 1-stage, 2nd order, stage order 1 MIRK formula
c
           s = 1
c
           c(1) = 1.0d0/2.0d0
c
           v(1) = 1.0d0/2.0d0
c
           b(1) = 1.0d0
c
c          Column 1 of x
c
           x(1) = 0.0d0
c
       elseif (method .EQ. 221) then
c
c          Define the 2-stage, 2nd order, stage order 1 MIRK formula
c
           s = 2
c
           c(1) = 0.0d0
           c(2) = 1.0d0
c
           v(1) = 0.0d0
           v(2) = 1.0d0
c
           b(1) = 1.0d0/2.0d0
           b(2) = 1.0d0/2.0d0
c
c          Column 1 of x
c
           x(1) = 0.0d0
           x(2) = 0.0d0
c
c          Column 2 of x
c
           x(3) = 0.0d0
           x(4) = 0.0d0
c
       elseif (method .EQ. 232) then
c
c          Define the 2-stage, 3rd order, stage order 2 MIRK formula
c
           s = 2
c
           c(1) = 0.0d0
           c(2) = 2.0d0/3.0d0
c
           v(1) = 0.0d0
           v(2) = 4.0d0/9.0d0
c
           b(1) = 1.0d0/4.0d0
           b(2) = 3.0d0/4.0d0
c
c          Column 1 of x
c
           x(1) = 0.0d0
           x(2) = 2.0d0/9.0d0
c
c          Column 2 of x
c
           x(3) = 0.0d0
           x(4) = 0.0d0
c
       elseif (method. EQ. 2321) then
c
c          Define the reflection of the above 2-stage, 3rd order, stage order 2
c          MIRK formula
c
           s = 2
c
           c(1) = 1.0d0
           c(2) = 1.0d0/3.0d0

           v(1) = 1.0d0
           v(2) = 5.0d0/9.0d0

           b(1) = 1.0d0/4.0d0
           b(2) = 3.0d0/4.0d0
c
c          Column 1 of x
c
           x(1) = 0.0d0
           x(2) = -2.0d0/9.0d0
c
c          Column 2 of x
c
           x(3) = 0.0d0
           x(4) = 0.0d0
c
       elseif (method .EQ. 343) then
c
c          Define the 3-stage, 4th order, stage order 3 MIRK formula
c
           s = 3
c
           c(1) = 0.0d0
           c(2) = 1.0d0
           c(3) = 1.0d0/2.0d0
c
           v(1) = 0.0d0
           v(2) = 1.0d0
           v(3) = 1.0d0 / 2.0d0
c
           b(1) = 1.0d0 / 6.0d0
           b(2) = 1.0d0 / 6.0d0
           b(3) = 2.0d0 / 3.0d0
c
c          Column 1 of x
c
           x(1) = 0.0d0
           x(2) = 0.0d0
           x(3) = 1.0d0 / 8.0d0
c
c          Column 2 of x
c
           x(4) = 0.0d0
           x(5) = 0.0d0
           x(6) = -1.0d0 / 8.0d0
c
c          Column 3 of x
c
           x(7) = 0.0d0
           x(8) = 0.0d0
           x(9) = 0.0d0
c
       elseif (method. EQ. 453) then
c
c          Define the 4-stage, 5th order, stage order 3 MIRK formula
c
           s = 4
c
           c(1) = 0.0d0
           c(2) = 1.0d0
           c(3) = 3.0d0/4.0d0
           c(4) = 3.0d0/10.0d0
c
           v(1) = 0.0d0
           v(2) = 1.0d0
           v(3) = 27.0d0/32.0d0
           v(4) = 837.0d0/1250.0d0
c
           b(1) = 5.0d0/54.0d0
           b(2) = 1.0d0/14.0d0
           b(3) = 32.0d0/81.0d0
           b(4) = 250.0d0/567.0d0
c
c          Column 1 of x
c
           x(1) = 0.0d0
           x(2) = 0.0d0
           x(3) = 3.0d0/64.0d0
           x(4) = 21.0d0/1000.0d0
c
c          Column 2 of x
c
           x(5) = 0.0d0
           x(6) = 0.0d0
           x(7) = -9.0d0/64.0d0
           x(8) = 63.0d0/5000.0d0
c
c          Column 3 of x
c
           x(9)  = 0.0d0
           x(10) = 0.0d0
           x(11) = 0.0d0
           x(12) = -252.0d0/625.0d0
c
c          Column 4 of x
c
           x(13) = 0.0d0
           x(14) = 0.0d0
           x(15) = 0.0d0
           x(16) = 0.0d0
c
       elseif (method. EQ. 4531) then
c
           s = 4
c
           c(1) = 0.0d0
           c(2) = 1.0d0
           c(3) = 3.0d0/4.0d0
           c(4) = 3.0d0/10.0d0
c
           v(1) = 0.0d0
           v(2) = 1.0d0
           v(3) = 27.0d0/32.0d0
           v(4) = 837.0d0/1250.0d0
c
           b(1) = 5.0d0/54.0d0
           b(2) = 1.0d0/14.0d0
           b(3) = 32.0d0/81.0d0
           b(4) = 250.0d0/567.0d0
c
c          Column 1 of x
c
           x(1) = 0.0d0
           x(2) = 0.0d0
           x(3) = 3.0d0/64.0d0
           x(4) = 21.0d0/1000.0d0
c
c          Column 2 of x
c
           x(5) = 0.0d0
           x(6) = 0.0d0
           x(7) = -9.0d0/64.0d0
           x(8) = 63.0d0/5000.0d0
c
c          Column 3 of x
c
           x(9)  = 0.0d0
           x(10) = 0.0d0
           x(11) = 0.0d0
           x(12) = -252.0d0/625.0d0
c
c          Column 4 of x
c
           x(13) = 0.0d0
           x(14) = 0.0d0
           x(15) = 0.0d0
           x(16) = 0.0d0
c
       elseif (method. EQ. 563) then
c
c          Define the 5-stage, 6th order, stage order 3 MIRK formula
c
           s = 5
c
           c(1) = 0.0d0
           c(2) = 1.0d0
           c(3) = 1.0d0/4.0d0
           c(4) = 3.0d0/4.0d0
           c(5) = 1.0d0/2.0d0
c
           v(1) = 0.0d0
           v(2) = 1.0d0
           v(3) = 5.0d0 / 32.0d0
           v(4) = 27.0d0 / 32.0d0
           v(5) = 1.0d0 / 2.0d0
c
           b(1) = 7.0d0 / 90.0d0
           b(2) = 7.0d0 / 90.0d0
           b(3) = 32.0d0 / 90.0d0
           b(4) = 32.0d0 / 90.0d0
           b(5) = 12.0d0 / 90.0d0
c
c          Column 1 of x
c
           x(1) = 0.0d0
           x(2) = 0.0d0
           x(3) = 9.0d0 / 64.0d0
           x(4) = 3.0d0 / 64.0d0
           x(5) = -5.0d0 / 24.0d0
c
c          Column 2 of x
c
           x(6) = 0.0d0
           x(7) = 0.0d0
           x(8) = -3.0d0 / 64.0d0
           x(9) = -9.0d0 / 64.0d0
           x(10) = 5.0d0 / 24.0d0
c
c          Column 3 of x
c
           x(11) = 0.0d0
           x(12) = 0.0d0
           x(13) = 0.0d0
           x(14) = 0.0d0
           x(15) = 2.0d0 / 3.0d0
c
c          Column 4 of x
c
           x(16) = 0.0d0
           x(17) = 0.0d0
           x(18) = 0.0d0
           x(19) = 0.0d0
           x(20) = -2.0d0 / 3.0d0
c
c          Column 5 of x
c
           x(21) = 0.0d0
           x(22) = 0.0d0
           x(23) = 0.0d0
           x(24) = 0.0d0
           x(25) = 0.0d0
c
       elseif (method. EQ. 5631) then
c
c          Define the improved 5-stage, 6th order, stage order 3 MIRK formula
c
           s = 5
c
           c(1) = 0.0d0
           c(2) = 1.0d0
           c(3) = 1.0d0/2.0d0-sqrt(21.0d0)/14.0d0
           c(4) = 1.0d0/2.0d0+sqrt(21.0d0)/14.0d0
           c(5) = 1.0d0/2.0d0
c
           v(1) = 0.0d0
           v(2) = 1.0d0
           v(3) = 1.0d0/2.0d0-9.0d0*sqrt(21.0d0)/98.0d0
           v(4) = 1.0d0/2.0d0+9.0d0*sqrt(21.0d0)/98.0d0
           v(5) = 1.0d0 / 2.0d0
c
           b(1) = 1.0d0 / 20.0d0
           b(2) = 1.0d0 / 20.0d0
           b(3) = 49.0d0 / 180.0d0
           b(4) = 49.0d0 / 180.0d0
           b(5) = 16.0d0 / 45.0d0
c
c          Column 1 of x
c
           x(1) = 0.0d0
           x(2) = 0.0d0
           x(3) = 1.0d0/14.0d0+sqrt(21.0d0)/98.0d0
           x(4) = 1.0d0/14.0d0-sqrt(21.0d0)/98.0d0
           x(5) = -5.0d0 / 128.0d0
c
c          Column 2 of x
c
           x(6) = 0.0d0
           x(7) = 0.0d0
           x(8) = -1.0d0/14.0d0+sqrt(21.0d0)/98.0d0
           x(9) = -1.0d0/14.0d0-sqrt(21.0d0)/98.0d0
           x(10) = 5.0d0 / 128.0d0
c
c          Column 3 of x
c
           x(11) = 0.0d0
           x(12) = 0.0d0
           x(13) = 0.0d0
           x(14) = 0.0d0
           x(15) = 7.0d0*sqrt(21.0d0)/128.0d0
c
c          Column 4 of x
c
           x(16) = 0.0d0
           x(17) = 0.0d0
           x(18) = 0.0d0
           x(19) = 0.0d0
           x(20) = -7.0d0*sqrt(21.0d0)/128.0d0
c
c          Column 5 of x
c
           x(21) = 0.0d0
           x(22) = 0.0d0
           x(23) = 0.0d0
           x(24) = 0.0d0
           x(25) = 0.0d0
c
       endif
c
       return
       end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine sol_eval(t,z,z_prime,fail,iwork,work)
c
c***overview***
c
c          This routine uses the current solution approximation, `Y',
c       obtained on the mesh stored in `mesh', plus the stages of the
c       underlying discrete RK method plus some new stages associated with
c       the interpolant in order to evaluate the interpolant.  This interpolant
c       is evaluated at `t', with the value returned in `z', and the value of
c       the first derivative returned in `z_prime'.
c
c          The interpolant is computed as follows: Given `t', assumed to
c       be somewhere in the problem interval, the subinterval of `mesh' that
c       contains `t' is determined first, i.e. find i such that mesh(i-1)
c       =< t < mesh(i) for t < mesh(Nsub) or i = Nsub for t = mesh(Nsub).
c       The fraction of the way `t' is through the ith subinterval, a quantity
c       which we call `tau',  is also determined during this step. Then the
c       interpolant on this ith subinterval is evaluated - the interpolant has
c       the form z(mesh(i-1) + tau*hi) = yim1 + hi sum br(tau)*kr. The weights,
c       br(tau), can be computed with a call to `interp_weights'. `Nsub' is
c       the number of subintervals.
c
c          The remaining information required to define the interpolant
c       on the ith subinterval is the current solution approximation at the
c       (i-1)st mesh point, `yim1', which can be obtained from `Y', and the
c       stages associated with the ith subinterval. The first `s' stages,
c       contained in `k_discrete', come from the discrete formula and
c       remaining `s_star - s' stages, contained in `k_interp', are necessary
c       for the interpolant. Each stage is a vector of length `neqns' where
c       `neqns' is the number of differential equations.
c
c          `k_discrete' is the vector containing the discrete stages for
c       all subintervals. The ith set of `s*neqns' locations of `k_discrete'
c       contains the `s' stages, each of length `neqns' corresponding to the
c       ith subinterval. `k_interp' plays a similar role for the `s_star-s'
c       stages associated with the interpolant.
c------------------------------------------------------------------------------
c***declaration of constants***
        integer Mxs
        parameter(Mxs=10)
c
c       `Mxs'    is the maximum value for the total number of stages
c                of the interpolant.
c
c***declaration of parameters***
c      imports:
        double precision t
c       `t'     the point at which the interpolant is to be evaluated.
c      exports:
c       double precision z(neqns)
        double precision z(*)
c
c       double precision z_prime(neqns)
        double precision z_prime(*)
c
        integer fail
c
c       `z'       is the value of the interpolant at x.
c
c       `z_prime' is the value of the first derivative of the
c                 interpolant at x.
c
c       `fail'    if fail = 1 then the point `t', at which the interpolant
c                 is to be evaluated, is outside the problem interval and
c                 therefore `sol_eval' will not produce a solution
c                 approximation. If fail = 0 then `t' is a valid meshpoint.
c
c***LAYOUT OF WORK ARRAYS***
c
c        All of the input info to this routine except for the evaluation
c        point `t' is contained within the `iwork' and `work' arrays.
c
         integer iwork(5)
c
c        iwork(1) = neqns
         integer neqns
c        `neqns' is the number of differential equations.
c
c        iwork(2) = Nsub
         integer Nsub
c        `Nsub' is the number of subintervals.
c
c        iwork(3) = s
         integer s
c        `s' is the number of discrete stages.
c
c        iwork(4) = s_star
         integer s_star
c        `s_star'  is the total number of stages required to form the
c                  interpolant on each subinterval. It includes all the
c                  stages of the discrete formula plus the additional
c                  stages required for the interpolant.
c
c        iwork(5) = method
         integer method
c        `method' defines the Runge-Kutta method that was used to compute the
c                 solution. It is passed via the /rk_method/ common block to
c                 the `interp_weights' routine.
c
c
c        let MnN = Mxs*neqns*Nsub.
c
c        double precision work(2*Mxs*neqns*Nsub+(Nsub+1)+(Nsub+1)*neqns)
         double precision work(*)
c
c        work(1..MnN) = k_discrete(Mxs*neqns*Nsub)
c        `k_discrete' contains the `s' discrete stage values for each
c                                                                subinterval.
c        work(MnN+1..2*MnN) = k_interp(Mxs*neqns*Nsub)
c        `k_interp' contains the `s_star-s' extra stage values for each
c                                 subinterval need to form the interpolant.
c
c        work(2*MnN+1..2*MnN+(Nsub+1)) = mesh(0:Nsub)
c        `mesh'  is the list of points defining the subintervals
c
c        work(2*MnN+(Nsub+1)+1..2*MnN+(Nsub+1)+(Nsub+1)*neqns)) =
c                                                   Y(neqns*(Nsub+1))
c        `Y'     contains the discrete solution at the mesh points;
c                the ith set of `neqns' locations of `Y' contains
c                the solution approximation at the ith mesh point.
c
c***declaration of local variables***
        integer i
        double precision hi
        double precision tau
        double precision weights(Mxs)
        double precision weights_prime(Mxs)
        double precision a, b
c
c       `i'             index for subinterval containing `t', i ranges
c                       from 1 to `Nsub'
c
c       `hi'            the size of the ith subinterval.
c
c       `tau'           equal to (t-mesh(i-1))/hi
c
c       `weights'       contains the weight polynomials of the interpolant
c                       evaluated at `tau'.
c
c       `weights_prime' contains the first derivatives of the weight
c                       polynomials of the interpolant evaluated at `tau'.
c
c       `a'             the left hand endpoint of the interval.
c       `b'             the right hand endpoint of the interval.
c
c***declaration of work array indexes***
        integer i_k_discrete
        integer i_k_interp
        integer i_mesh, i_Y

c       `i_k_discrete' = 1; `i_k_interp' = Mxs*neqns*Nsub+1;
c       `i_mesh' = 2*Mxs*neqns*Nsub+1; `i_Y' = 2*Mxs*neqns*Nsub+(Nsub+1)+1

c***declaration of variables in /rk_method/ common block***
c      exports:
c       integer   method
        common   /rk_method/ method

c       `method'  defines the method to be used
c
c       This value is passed through the rk_method common block
c       to the `interp_weights' routine.
c-----------------------------------------------------------------------------
c      called by:`main',    user's main program subsequent to the
c                           completion of the call to mirkdc.
c
c      calls to: `interval', `interp_weights', `sum_stages'
c------------------------------------------------------------------------------
c         Setup access to info in `iwork' and `work'

          neqns = iwork(1)
          Nsub = iwork(2)
          s = iwork(3)
          s_star = iwork(4)
          method = iwork(5)

          i_k_discrete = 1
          i_k_interp = Mxs*neqns*Nsub+1
          i_mesh = 2*Mxs*neqns*Nsub+1
          i_Y = 2*Mxs*neqns*Nsub+(Nsub+1)+1

          a = work(i_mesh)
          b = work(i_mesh+(Nsub))
c
          fail = 0
          if ((t.LT.a).or.(t.GT.b)) then
           fail = 1
          else
c          Call `interval' to determine the subinterval containing t.
c          The index of the subinterval is returned in `i'.
c
           call interval(Nsub,work(i_mesh),t,i)
c
c          Compute `tau' for use in the evaluation of the weight polynomials.
c
           hi = work(i_mesh+i) - work(i_mesh+(i-1))
           tau = (t - work(i_mesh+(i-1)))/hi
c
c          Setup the weights, evaluated at `tau'.
c
           call interp_weights(s_star, tau, weights,weights_prime)
c
c          We need to evaluate the interpolant and its derivative at `tau'.
c
c          This is done by taking weighted sums of the stages, using the weights
c          we have just computed, and the stages corresponding to the ith
c          subinterval, and calling the `sum_stages' routine.
c
           call sum_stages(neqns,hi,work(i_Y+(i-1)*neqns),s,
     +               work(i_k_discrete+(i-1)*s*neqns),s_star,
     +               work(i_k_interp+(i-1)*(s_star-s)*neqns),
     +                     weights, weights_prime, z, z_prime)
c
          end if
          end

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       subroutine subcom(neqns,hi,ti,yi,yip1,phii,k,y,fsub)
c
c***overview***
c
c  ***Subcom Routine(Calculation of subcomponents of residual function)***
c
c            The purpose of this routine is to evaluate the component
c       of PHI, the residual function, associated with the i-th sub-
c       interval. This component of PHI is returned through the parameter
c       `phii'.
c                                          s
c                                          --
c                     phii = yip1-yi- hi * <  br*kr
c                                          --
c                                          r=1
c
c       The evaluation takes place for the current iterate `Y'.
c       However, the i-th component of PHI depends only on the solution
c       information directly associated with the i-th subinterval, namely
c       yi and yi+1, which have been assigned to `yi' and `yip1'.  Also
c       the size of the i-th subinterval is assigned to `hi' and the left
c       hand endpoint of the ith subinterval is assigned to `ti'.
c
c            The calculation uses the s-stage MIRK formula defined by the
c       values of `s', `c',`v', `b' and `x' (as passed in from the
c       /RK_s/,/RK_coeff/ commom blocks), The solution approximation vectors,
c       `yi' and `yip1', as well as `phii', are of size `neqns', the number
c       of differential equations.
c
c            During the computation of `phii', we will compute `s'
c       intermediate vectors, kr, r=1..s, each of size `neqns', and called
c       the stages of the formula. The stages are stored in the  vector `k'
c       in order, and are returned by this routine for later use during the
c       call to `newmat' and during the interpolation phase.
c       `y' provides a work array of size `neqns'.
c-----------------------------------------------------------------------------
c***declaration of constants***
         integer Mxs
         parameter (Mxs=10)
c
c         `Mxs' maximum number of stages.
c
c***declaration of parameters***
c        imports :
           integer                    neqns
           double precision           hi, ti, yi(neqns), yip1(neqns)
c
c         `neqns' is the number of differential equations
c
c         `hi'    is the size of the ith subinterval
c
c         `ti'    is the lefthand endpoint of the ith subinterval
c
c         `yi'    is the solution approximation value at the left endpoint
c
c         `yip1'  is the solution approxiamtion at the right endpoint
c        exports :
           double precision           phii(neqns)
           double precision           k(Mxs*neqns)
c
c         `phii' is the value of the component of the residual function
c                associated with the ith subinterval
c
c         `k'    stores the `s' stage values k1, k2, ..., ks, each of size
c                `neqns'
c
c***declaration of work array***
           double precision           y(neqns)

c          `y'     a work array for temporary storage.
c
c***user-supplied subroutines***
           external fsub
c
c          `fsub' - This routine defines f(t,y) for the first order system
c                   of differential equations, y' = f(t,y).
c
c***declaration of local variables***
           integer                    j, r, l
           integer                    i1,i2
           double precision           tr
c
c          `j,r,l' are loop index variables; `j,r' range from 1 to `s'
c                 `l' ranges from 1 to `neqns'.
c
c          `i1',`i2' are indexes used in calculating ur.
c
c          `tr'    point at which rth stage of Runge-Kutta method is
c                  evaluated.
c
c***declaration of variables in common block /RK_s/***
            integer                 s
            common /RK_s/ s
c
c           `s' is the number of discrete stages of the Runge-Kutta formula.

c***declaration of variables in common block /RK_coeff/***
            double precision        v(Mxs), b(Mxs)
            double precision        c(Mxs), x(Mxs**2)
            common /RK_coeff/ v, b, c, x
c
c           `v', `b', `c',and `x' are coefficients that define
c           the discrete Runge-Kutta formula.
c
c           The first component of `c' contains the abscissa for the first
c           new stage. The second component of `c' contains the abscissa for
c           the second new stage, and so on.
c
c           The first component of `v' contains the blending coefficient
c           for the first new stage. The second component of `v' contains
c           the blending coefficient for the second new stage, and so on.
c
c           The matrix, X, is related to  the data structure, `x'
c           as follows (the matrix is stored columnwise in the vector):
c
c                 X(i,j) = x( (j-1)*s + i )
c
c           `x' contains the coupling coefficients for each of the new
c           stages
c------------------------------------------------------------------------------
c      called by: `resid'
c      calls to: `fsub'
c------------------------------------------------------------------------------
c***COMPUTATION OF `phii'***
c
c       For each stage of the Runge-Kutta method, we compute the
c       argument of the function f, the right hand side of the ODE system,
c       and then evaluate f there by calling `fsub', with the results being
c       stored in the appropiate part of the vector `k'. The vector `y'
c       provides  storage for accumulation of the argument for f.
c
        do  7  r = 1, s
c
           do 2 l = 1, neqns
               y(l) = 0.0d0
 2         continue
c
c          Add the weighted sum of the previous stages
c
           do 4 j = 1, r-1
            i1 = (j-1)*s + r
            i2 = (j-1)*neqns
c
             do 3 l = 1, neqns
                 y(l) = y(l) + x(i1)*k(i2+l)
 3           continue
c
 4         continue
c
c         Multiply by hi
c
          do 5 l = 1, neqns
            y(l) = hi*y(l)
 5        continue

c          Add (1-v(r))*yi + v(r)*yip1
c
           do  6  l = 1, neqns
              y(l) = (1.0d0 - v(r))*yi(l) + v(r)*yip1(l) + y(l)
 6         continue
c
c             The second argument to f has been accumulated in the vector `y'.
c             Compute the abscissa for the argument to fsub.
c
              tr = ti + c(r)*hi
c
c            Evaluate the function f to get kr, the rth stage.
c            The user routine which evaluates the right hand side of the
c            system of ODE's must be called `fsub'
c
              call fsub(neqns, tr, y, k((r-1)*neqns+1))
c
c            The rth stage has now been computed and stored in the rth
c            vector location of k
c
 7      continue
c
c      All `s' stages of the formula have now been computed. We can now
c      take their weighted sum in order to complete the calculation of `phii'.
c
        do  8 l = 1, neqns
            phii(l) = 0.0d0
 8      continue
c
        do 9 r = 1, s
         i1 = (r-1)*neqns
         do 10 l = 1, neqns
            phii(l) = phii(l) - b(r)*k(i1+l)
 10         continue
 9      continue
c
c       Multiply phii by hi
c
        do 11 l = 1, neqns
          phii(l) = hi*phii(l)
 11     continue
c
c       Add  (yip1 - yi) to phii
c
        do  12  l = 1, neqns
           phii(l) = yip1(l) - yi(l) + phii(l)
 12     continue
c
        return
        end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine sum_stages(neqns,hi,yim1,s,ki_discrete,
     +            s_star,ki_interp,weights,weights_prime,z,z_prime)
c
c***overview***
c
c          Given the stages and corresponding weights, this routine
c      computes the value of the interpolant, `z', and its first derivative,
c      `z_prime', as follows:
c            z = yim1 + hi sum weights(r)*ki(r),
c            z_prime = sum weights_prime(r)*ki(r),
c      where the ki(r)'s are the stages. The first `s' stages, contained
c      in `ki_discrete', come from the discrete formula; the remaining
c      `s_star - s' stages, contained in `ki_interp', are necessary for
c      the interpolant only. Each stage is a vector of length `neqns' where
c      `neqns' is the number of differential equations. `hi' is the size of
c      the ith subinterval and `yim1'  is the discrete solution at the left
c      endpoint of the ith subinterval. `weights' contains the values of the
c      weight polynomials of the interpolant and `weights_prime' contains the
c      values of the first derivative of the weight polynomials.
c------------------------------------------------------------------------------
c***declaration of parameters***
c      imports:
            integer neqns
            double precision hi
            double precision yim1(neqns)
            integer s
            double precision ki_discrete(s*neqns)
            integer s_star
            double precision ki_interp((s_star-s)*neqns)
            double precision weights(s_star)
            double precision weights_prime(s_star)
c
c           `neqns'  is the number of differential equations.
c
c           `hi'     is the size of the current (ith) subinterval.
c
c           `yim1'   is the discrete solution at the left endpoint
c                    of the ith subinterval.
c
c           `s'      is the number of discrete stages.
c
c           `ki_discrete'   the set of `s' discrete stages associated with
c                           the ith subinterval.
c
c           `s_star'        is the total number of stages.
c
c           `ki_interp'     the set of `s_star-s' interpolation stages
c                           associated with the ith subinterval.
c
c           `weights'       the weights associated with the stages needed to
c                           evaluate the interpolant at the appropriate
c                           point. The weights have aleady been evaluated at
c                           that point.
c
c           `weights_prime' the weights associated with the stages needed to
c                           form the first derivative of the interpolant
c                           at the appropriate point. The weights have already
c                           been evaluated at this point.
c      exports:
            double precision z(neqns)
            double precision z_prime(neqns)
c
c           `z'       is the value of the interpolant.
c
c           `z_prime' is the value of the first derivative of the
c                     interpolant.
c
c***declaration of local variables***
            integer j
c
c           `j'       loop index from 1 to 's', or 1 to 's_star-s'
c---------------------------------------------------------------------------
c      Called by: `defect_estimate', `interp_eval'
c      Calls to : `dcopy', `daxpy'
c---------------------------------------------------------------------------
c          Initialize all vectors to be used in this computation to zero.
c
           call dcopy(neqns,0.0d0,0,z,1)
           call dcopy(neqns,0.0d0,0,z_prime,1)
c
c          We can now sum the appropriately weighted stages to compute
c          the interpolant `z' and the first derivative of the
c          interpolant, `z_prime'.
c
c          Accumulate weighted sum of discrete stages for `z' and `z_prime'.
c
           do 10 j = 1, s
c
                call daxpy(neqns,weights(j),
     +                      ki_discrete((j-1)*neqns+1),1,z,1)
                call daxpy(neqns,weights_prime(j),
     +                      ki_discrete((j-1)*neqns+1),1,z_prime,1)
 10        continue
c
c          Accumulate weighted sum of new interpolant stages for `z'
c          and `z_prime'.
c
           do 20 j = 1,(s_star-s)
c
            call daxpy(neqns,weights(s+j),
     +                      ki_interp((j-1)*neqns+1),1,z,1)

            call daxpy(neqns,weights_prime(s+j),
     +                      ki_interp((j-1)*neqns+1),1,z_prime,1)

 20        continue
c
c          The calculation of `z_prime' is now complete.
c          Scale the weighted sum of stages by hi, and then add
c          on yim1. This will complete the calculation of `z'.
c
           call dscal(neqns,hi,z,1)
           call daxpy(neqns,1.0d0,yim1,1,z,1)
c
           return
           end