c-----------------------------------------------------------------------
c Use preprocessor directives for certain aspects of this driver.
#ifndef NPDE
#  define NPDE 1
#endif
#ifndef NINTMX
#  define NINTMX 1000
#endif
#ifndef XA
#   define XA 0d0
#endif
#ifndef XB
#   define XB 1d0
#endif
#ifndef TOUT
#   define TOUT 1d0
#endif
c Define BACOL to use this driver with the old BACOL code, otherwise
c BACOLI will be used.
c Define DIRICHLET if the boundary conditions are such.
c burg1.f and burg2.f have Dirichlet boundary conditions, for example.
c-----------------------------------------------------------------------
        implicit none
c constants:
        integer                 kcol, kcolmx
c                               kcol is the number of collocation points
c                               to be used in each subinterval, which is
c                               equal to the degree of the piecewise
c                               polynomials minus one.
c                                  3 < kcol < 11 for the loi,
c                                  2 < kcol < 11 for the sci.
        parameter              (kcolmx = 10)
        integer                 npde
c                               number of pdes
        parameter              (npde = NPDE)
        integer                 nintmx
c                               maximum number of intervals required
        parameter              (nintmx = NINTMX)
        integer                 maxvec
c                               the dimension of the vector of
c                               B-spline coefficients
        parameter              (maxvec = npde*(nintmx*kcolmx+2))
        integer                 lrp
c                               uncomment the lines with a `-'
c                               continuation character
c                               and fix the parentheses if you
c                               only want to use the loi estimate.
#ifdef BACOL
        parameter              (lrp=134+nintmx*(35+35*kcolmx+31*npde
     +                               +38*npde*kcolmx+8*kcolmx*kcolmx)
     +                               +14*kcolmx+79*npde+npde*npde*(21
     +                               +4*nintmx*kcolmx*kcolmx+12*nintmx
     +                               *kcolmx+6*nintmx))
#else
        parameter              (lrp=113+57*npde+27*nintmx
     +                               +13*npde*npde+9*kcolmx+24*kcolmx
     +                               *nintmx+6*nintmx*kcolmx*kcolmx
     +                               +7*nintmx*npde+27*npde*nintmx
     +                               *kcolmx+2*npde*npde*nintmx*kcolmx
     +                               *kcolmx+4*npde*npde*kcolmx*nintmx)
c    -                               -15*nintmx+3*kcolmx+kcolmx*kcolmx
c    -                               -8*kcolmx*nintmx-3*nintmx*npde
c    -                               -nintmx*kcolmx*kcolmx+2
#endif
        integer                 lip
c                               see the comment for ipar
#ifdef BACOL
        parameter              (lip=115+npde*(nintmx*(2*kcolmx+1)+4))
#else
        parameter              (lip=101+npde*(nintmx*kcolmx)+2)

        external                f
        external                uinit
        external                derivf
        external                bndxa
        external                bndxb
        external                difbxa
        external                difbxb
#endif
        double precision        xa
c                               the left boundary point
        parameter              (xa = XA)
        double precision        xb
c                               the right boundary point
        parameter              (xb = XB)
c-----------------------------------------------------------------------
c bacoli input:
        double precision        t0
c                               t0 < tout is the initial time.
        double precision        tout
c                               tout is the desired final output time.
        double precision        atol(npde)
c                               atol is the absolute error tolerance
c                               request and is a scalar quantity if
c                               mflag(2) = 0.
        double precision        rtol(npde)
c                               rtol is the relative error tolerance
c                               request and is a scalar quantity if
c                               mflag(2) = 0.
        integer                 nint
c                               nint is the number of subintervals
c                               defined by the spatial mesh x.
        double precision        x(nintmx+1)
c                               x is the spatial mesh which divides the
c                               interval [x_a,x_b] as: x_a = x(1) <
c                               x(2) < x(3) < ... < x(nint+1) = x_b.
        integer                 mflag(12)
c                               this vector of user input determines
c                               the interaction of bacoli with dassl.
c                               PDE system definition subroutines.
c-----------------------------------------------------------------------
c bacoli work storage:
        double precision        rpar(lrp)
c                               rpar is a floating point work array
c                               of size lrp.
        integer                 ipar(lip)
c                               ipar is an integer work array
c                               of size lip.
c-----------------------------------------------------------------------
c bacoli output:
        double precision        y(maxvec)
c                               on successful return from bacoli, y is
c                               the vector of bspline
c                               coefficients at the current time t0.
        integer                 idid
c                               idid is the bacoli exit status flag
c                               which is based on the exit status from
c                               dassl on error checking performed by
c                               bacoli on initialization.
c-----------------------------------------------------------------------
        double precision        tstart
c                               stores initial variables that the
c                               subroutine bacoli will change, so they
c                               can be used in second call to bacoli
        integer                 est
c                               1 to use LOI, 2 to use SCI
        integer                 i
        character*5             eststr
c-----------------------------------------------------------------------
c problem specific variables:
        double precision        eps
        common /burger/         eps
c
        double precision        a1,a2,d1,d2,r,c,n,pe1,pe2
        common /rcdstuff/       a1,a2,d1,d2,r,c,n,pe1,pe2
c-----------------------------------------------------------------------
c additional counters
        integer ncrdmp, ncrslv
        common /counters/ ncrdmp, ncrslv
c-----------------------------------------------------------------------
c subroutines called:
c                               bacolri
c
c-----------------------------------------------------------------------
c  Initializations
      tstart = 0d0
      nint = 20
      t0 = tstart
      tout = TOUT

      a1 = 30.d0
      a2 = 30.d0
      d1 = 1.5d0
      d2 = 1.2d0
      r = 1000.d0
      c = 0.96d0
      n = 1.d0
      pe1 = 100.d0
      pe2 = 100.d0

      ncrdmp = 0
      ncrslv = 0
c-----------------------------------------------------------------------
c  Get some user input

      write(*,'(/a)') 'Pick a kcol'
      read(*,*) kcol

      write(*,'(/a)') 'Pick a stop time'
      read(*,*) tout

      write(*,'(/a)') 'Pick the error control mode'
      write(*,'(a)') '1. LOI/LE'
      write(*,'(a)') '2. SCI/ST'

      read(*,*) est

      write(*,'(/a)') 'Please input the tolerance to use'
      read(*,*) atol(1)
      rtol(1) = atol(1)

      write(*,'(/a)') 'Please input epsilon, too'
      read(*,*) eps

c-----------------------------------------------------------------------
c  Even more initializations

c  define the mesh based on a uniform step size.
      x(1) = xa
      do i = 2, nint
         x(i) = xa + ((i-1) * (xb - xa)) / nint
      end do
      x(nint+1) = xb

c  initialize the mflag vector.
      do i = 1, 8
         mflag(i) = 0
      end do
c  set mflag(5) = 1 for dirichlet boundary conditions, 0 otherwise
#ifdef DIRICHLET
      mflag(5) = 1
#endif

#ifndef BACOL
      if (est .eq. 1) then
         mflag(8) = 0
         eststr = 'LOI'
      elseif (est .eq. 2) then
         mflag(8) = 1
         eststr = 'SCI'
      endif
#else
      eststr = 'BAC'
#endif

c       This is to restrict the order used so that a solution is 
c        actually reached.      
#ifdef BACOLISCHROD
      mflag(7) = 1
      ipar(15) = 2
#endif

c-----------------------------------------------------------------------
c  Computation & output section

      write(*,'(/a)') 'THE INPUT IS'
      write(*,9001) 'kcol = ', kcol, ', nint0 = ', nint,
     &              ', npde = ', npde, ', tout = ', tout
 9001 format( a, i2, a, i4, a, i3, a, es7.1 )
      write(*,9002) 'atol = rtol = ', rtol(1),
     &              ',    ', eststr, ',    epsilon = ', eps
 9002 format( a, es7.1, 3a, es7.1 )

#ifndef BACOL
      call bacoli(t0, tout, atol, rtol, npde, kcol, nintmx, nint,
     &            x, mflag, rpar, lrp, ipar, lip, y, idid,
     &            f, derivf, bndxa, difbxa, bndxb, difbxb, uinit)
c       Function names aren't passed yet, uncomment when updated.
#else
      call bacol(t0, tout, atol, rtol, npde, kcol, nintmx, nint,
     &            x, mflag, rpar, lrp, ipar, lip, y, idid)
#endif

c     Check for errors running BACOLI.
      if (idid .le. 0) goto 700

c-----------------------------------------------------------------------
c  Sucessful return section

      write(*,'(a,i3)')     'IDID         = ', idid
      write(*,'(a,i10)')    'Final nint   = ', nint
      write(*,'(a,i10)')    'ipar(irshin) = ', ipar(9)
      write(*,'(a,i10)')    'ipar(isteps) = ', ipar(10)
      write(*,'(a,i10)')    'nremesh      = ', ipar(11)+ipar(9)
      write(*,'(a,i10)')    'ipar(icolds) = ', ipar(14)
      write(*,'(a,i10)')    'ncrdmp       = ', ncrdmp
      write(*,'(a,i10)')    'ncrlve       = ', ncrslv

      stop
c-----------------------------------------------------------------------
c  Unsucessful return section
  700 continue
      write(*,'(a)')        'CANNOT PROCEED DUE TO ERROR FROM BACOLI.'
      write(*,'(a,g25.18)') 'Failure at T = ', t0
      write(*,'(a,i3)')     'IDID = ', idid
      write(*,'(a,i10)')    'NINT = ', nint

      stop
c-----------------------------------------------------------------------
      end
