c-----------------------------------------------------------------------
c driver-L2.f
c With the true solution defined in TRUU, this driver uses DQAGE from
c quadpack to approximate the L2-norm of the error in the approximate
c solution after BACOLI has returned.
c Some bits need to be modified to suit the current test being done,
c such as the parameters npde, nintmx, xa and xb,
c or common block declarations and initializations for problem-specific
c parameters.
c Also the output given here probably isn't what you want, so do
c something about that.
c
c When prompted to input the tolerance or epsilon, it is recommended to
c use engineering notation (like `1e-7'.)
c
c This driver should be compiled with:
c   some version of BACOLI, eg, bacoli.f
c   a problem definition file, eg, burg1.f
c   quadpack routines
c
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                 npde
c                               number of pdes
        parameter              (npde = 4)
        integer                 nintmx
c                               maximum number of intervals required
        parameter              (nintmx = NINTMX)

        integer                 kcol
        parameter              (kcol=10)

        integer                 maxvec
c                               the dimension of the vector of
c                               B-spline coefficients
        parameter              (maxvec = npde*(nintmx*kcol+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.
        parameter              (lrp=134+nintmx*(35+35*kcol+31*npde
     +                               +38*npde*kcol+8*kcol*kcol)
     +                               +14*kcol+79*npde+npde*npde*(21
     +                               +4*nintmx*kcol*kcol+12*nintmx
     +                               *kcol+6*nintmx))

        integer                 lip
c                               see the comment for ipar

        parameter              (lip=115+npde*(nintmx*(2*kcol+1)+4))

        integer                 lenwrk
c                               the dimension of array work when we
c                               call values
        parameter              (lenwrk =(kcol+2)+kcol*(nintmx+1)+4)
        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                 nintinit
        integer                 nint
c                               nint is the number of subintervals
c                               defined by the spatial mesh x.
        double precision        x(nintmx+1)
        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.

        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                 i, m, k
c                               loop counters
        character*3             eststr
        integer                 nderiv
        parameter              (nderiv=0)
c-----------------------------------------------------------------------
c L2-norm error calculation:
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

      tstart = 0d0
      nintinit = 100
      tout = TOUT
      eps = 1d-3

      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

      atol(1) = 1d-20
      rtol(1) = atol(1)

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

c  vectors.
      do i = 1, 8
         mflag(i) = 0
      end do

c     Set error estimation scheme for BACOLI. For BACOL, set different
c     kcol value for use in call to values in order to extract higher
c     or lower-order interpolant.
       eststr = 'BAC'

c-----------------------------------------------------------------------
c  Computation & output section
      write(*,*) 'kcol = ', kcol, ', nint0 = ', nintinit,
     &              ', npde = ', npde, ', tout = ', tout
      write(*,*) 'atol = rtol = ', rtol(1),
     &              ',    ', eststr, ',    epsilon = ', eps

c      some per-run initialisations
       t0 = tstart
c         rpar(1) = tout part of a no longer existing flag
       nint = nintinit
       x(1) = xa
       do i = 2, nint
          x(i) = xa + ((i-1) * (xb - xa)) / nint
       end do
       x(nint+1) = xb
       mflag(1) = 0

      call bacol(t0, tout, atol, rtol, npde, kcol, nintmx, nint,
     &            x, mflag, rpar, lrp, ipar, lip, y, idid)

      open(100, file='reference_coefs.txt', status='NEW')
      open(101, file='reference_x.txt', status='NEW')
      open(102, file='parameters.txt', status = 'NEW')

      do i = 1, maxvec
        write(100,*) y(i)
      end do
      do i = 1, nint+1
        write(101,*) x(i)
      end do
      write(102, *) kcol
      write(102, *) nint

      close(100)
      close(101)
      close(102)

      stop
      end
