c-----------------------------------------------------------------------
c This driver file writes output for each PDE using BACOLR's mesh points
c at ntout uniform points along the temporal domain. tstart is 0, tstop
c prompted for and the spatial domain is [0,1]. A number of parameters
c are hard coded and should be edited as needed: npde, xa, xb, tstart,
c nintmx, kcol and nintinit.
c
c The output format used here allows the user to take advantage of
c BACOLR's spatial adaptivity in their data visualisation, to some
c extent, but the trade-off is that data points are slightly more
c difficult to work with than those over a uniform grid.
c
c If returning after uniform time steps is not desirable,
c it is possible to modify this code to have bacolr return after
c a specific number of time steps (as performed by DASSL).
c See the documentation inside bacolr.f for details.
c
c The subroutine values is also able to return spatial partial
c derivatives of the approximate solution (e.g., U_x, U_xx,)
c but that is not demonstrated here. Refer to the documentation
c of the values subroutine.
c
c When prompted to input a tolerance, it is recommended to
c use engineering notation (e.g., '1e-6').
c
c After a successful computation, this program will have written files
c named Points001, ..., PointsNPDE containing data points as ordered
c triples, one per line. The coordinates are ordered (X, T, U).
c
c This driver should be linked with with bacolr.f and a problem
c definition file, eg, burg1.f.
c
c Example Python plotting code:
c
c   import matplotlib as mpl
c   mpl.use('AGG')  # for systems not running a GUI
c   from mpl_toolkits.mplot3d import Axes3D
c   from matplotlib import cm
c   import matplotlib.pyplot as plt
c   import numpy as np
c
c   styling = {
c     'cmap': cm.coolwarm,
c     'linewidth': 0,
c     'antialiased': True
c   }
c
c   x, t, u = np.loadtxt('Points001', unpack=True)
c   fig = plt.figure()
c   ax = fig.add_subplot(111, projection='3d')
c   ax.plot_trisurf(x, t, u, **styling)
c
c   ax.set_xlabel('$x$')
c   ax.set_ylabel('$t$')
c   ax.set_zlabel('$u(t,x)$')
c
c   plt.savefig('trimesh.png')
c
c axes3d.plot_trisurf() was added to matplotlib in version 1.2.0.
c It drapes a surface over a triangularization of the data.
c-----------------------------------------------------------------------
        implicit none
c constants:
        integer                 npde
        parameter              (npde = 1)
c                               number of PDEs in the system.
        double precision        xa
        parameter              (xa = 0d0)
c                               the left spatial boundary point.
        double precision        xb
        parameter              (xb = 1d0)
c                               the right spatial boundary point.
        double precision        tstart
        parameter              (tstart = 0d0)
c                               the start of the time domain, used to
c                               initialise t0.
        integer                 nintinit
        parameter              (nintinit = 20)
c                               the number of subintervals in the
c                               initial uniform mesh.
        integer                 nintmx
        parameter              (nintmx = 1000)
c                               maximum number of subintervals that
c                               bacolr may use
        integer                 kcol
        parameter              (kcol = 3)
c                               kcol is the number of collocation
c                               points to be used in each subinterval.
c                                  2 <= kcol <= 10.
c                               A value of 3 is most efficient in
c                               typical cases, and higher when
c                               a sharp tolerance is enforced.
        integer                 maxvec
        parameter              (maxvec = npde*(nintmx*kcol+2))
c                               the maximum size of the vector of
c                               B-spline coefficients, y.
         integer                 lrp
        parameter              (lrp = 74+24*npde*npde*nintmx*kcol
     +                            + 8*npde*npde*nintmx*kcol*kcol
     +                            + 29*npde*nintmx*kcol+61*npde+14*kcol
     +                            + 35*nintmx*kcol
     +                            + 35*nintmx+21*nintmx*npde+8*nintmx
     +                                * kcol*kcol
     +                            + 37*npde*npde+12*npde*npde*nintmx)
c                               length of bacolr's floating point
c                               work array, rpar.
        integer                 lip
        parameter              (lip = 100+3*npde*(nintmx*(2*kcol+1)+4))
c                               length of bacolr's integer
c                               work array, ipar.
        integer                 lcp
        parameter              (lcp=npde*(4+(2*kcol+1)*nintmx)+npde*npde
     +                            * (8+nintmx*(2*kcol*(kcol+3)+3)))
c                               length of bacolr's complex
c                               work array, cpar.
        integer                 lenwrk
        parameter              (lenwrk = (kcol+2)+kcol*(nintmx+1)+4)
c                               length of the work array used by VALUES
c                               after a successful return from bacolr.
c-----------------------------------------------------------------------
c bacolr input:
        double precision        t0
c                               t0 < tout is the initial time for a
c                               given call to bacolr.
        double precision        tout
c                               tout is the next desired output time.
        double precision        tstop
c                               tout is the final desired 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 [xa,xb] as:
c                               xa = x(1) < x(2) < .. < x(nint+1) = xb
        integer                 mflag(7)
c                               this vector is used for more specific
c                               settings. see documentation in bacolr.f
c-----------------------------------------------------------------------
c bacolr 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.
        double complex          cpar(lcp)
c                               cpar is an complex work array
c                               of size lip.
c-----------------------------------------------------------------------
c bacolr output:
        double precision        y(maxvec)
c                               on successful return from bacolr, y is
c                               the vector of time dependent B-spline
c                               coefficients at the current time t0.
        integer                 idid
c                               idid is the bacolr exit status flag
c-----------------------------------------------------------------------
c miscellany:
        integer                 ntout
c                               number of evenly spaced time steps to
c                               give the computed solution values.
        double precision        uout(npde,nintmx)
c                               contains solution approximations,
c                               which are converted from their internal
c                               representation as B-spline coefficients
c                               by the subroutine values.
        double precision        valwrk(lenwrk)
c                               valwrk is the work array for values.
        integer                 i, j, k
c                               loop indices
        character*3             eststr
c                               eststr is printed in the input summary.
        character*32            fname
c                               fname is a string formatted as
c                               Points001, .., PointsNPDE.
c-----------------------------------------------------------------------
c problem specific variables:
        double precision        eps
        common /burger/         eps
c                               eps is used in the burg1.f, burg2.f and
c                               Cahn_Allen.f problem definition files.
c
        double precision        a1,a2,d1,d2,r,c,n,pe1,pe2
        common /rcdstuff/       a1,a2,d1,d2,r,c,n,pe1,pe2
c                               used in RCDsys.f.
c-----------------------------------------------------------------------
c subroutines called:
c                               bacolr
c                               values
        external                uinit
c-----------------------------------------------------------------------
c  Common block initializations

      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 = 10000.d0
      pe2 = 10000.d0

c-----------------------------------------------------------------------
c  Get some user input

      write(*,*) 'Enter tstop, the end of the temporal domain.'
      read(*,*,err=999) tstop

      write(*,*) 'At how many equally-spaced points along the '
     &        // 'time domain is output desired?'
      read(*,*,err=999) ntout

      write(*,*) 'Please input the error tolerance.'
      read(*,*,err=999) atol(1)
      rtol(1) = atol(1)

c-----------------------------------------------------------------------
c  BACOLR parameter initializations

      nint = nintinit
      t0 = tstart

c  set the initial mesh to be uniformly spaced
      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, 7
         mflag(i) = 0
      end do
c  set mflag(5) = 1 for Dirichlet boundary conditions, 0 otherwise
c  recall that Dirichlet boundary conditions involve only function
c  values, i.e., U_x is not involved in bndxa or bndxb.
      mflag(5) = 1

      do i = 1, lip
         ipar(i) = 0
      end do

      do i = 1, lrp
         rpar(i) = 0d0
      end do

      eststr = 'BAC'

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

      write(*,*)     'THE INPUT IS'
      write(*,9000) 'kcol = ', kcol, ', nint0 = ', nint,
     &              ', npde = ', npde, ', tout = ', tstop
      write(*,9001) 'atol = ', atol(1), ', rtol = ', rtol(1),
     &              ',                 ', eststr

c     open files to receive output.
      do k = 1, npde
c        String manipulation sure is awkward in F77.
c        F90/95 intrinsics and/or format specifiers are recommended.
         if (k .le. 9) then
            write(fname,'("Points00",i1)') k
         else if (k .le. 99) then
            write(fname,'("Points0",i2)') k
         else if (k .le. 999) then
            write(fname,'("Points",i3)') k
         end if
         open(unit=10+k,file=fname)
      end do

c     Call BACOLR at each tout and save output to files.
      do j = 2, ntout

c        idid must be set to 1 for continuation calls.
c        On the initial call, it does not matter if it is 1
         idid = 1

c        tout is the return time for this bacolr call
         tout = tstart + (j-1)*(tstop-tstart)/(ntout-1)

         call bacolr(t0, tout, atol, rtol, npde, kcol, nintmx,
     &               nint, x, mflag, rpar, lrp, ipar, lip, cpar, lcp,
     &               y, idid)
c        Check for errors running BACOLR.
         if (idid .le. 0) goto 700

c        After the first call to bacolr, the spatial mesh will
c        have been adapted to the PDE system.
c        Use this mesh, and the user subroutine uinit, to output
c        data points for tstart.
         if (j .eq. 2) then
            do i = 1, nint+1
               call uinit(x(i), uout(1, i), npde)
            end do
            do k = 1, npde
               do i = 1, nint+1
                  write(10+k,*) x(i), tstart, uout(k, i)
               end do
            end do
         end if

c        Write computed solution values to files for each
c        PDE system component. Evaluation is done at mesh points
c        so that the output benefits from bacolr's spatial adaptivity.
         call values(kcol, x, nint, x, npde, nint+1, 0, uout, y,
     &               valwrk)
         do k = 1, npde
            do i = 1, nint+1
               write(10+k,*) x(i), tout, uout(k, i)
            end do
         end do

c        mflag(1) must be 0 on the initial call and 1 on continuation.
         mflag(1) = 1
      end do

c-----------------------------------------------------------------------
c  Sucessful return section
      write(*,9002) 'IDID       = ', idid
      write(*,9002) 'Final nint = ', nint
      write(*,9002) 'nsteps     = ', ipar(10)

  600 continue
c     close files
      do k = 1, npde
         close(10+k)
      end do

      stop
c-----------------------------------------------------------------------
c  Unsuccessful return section
  700 continue
      write(*,*) 'CANNOT PROCEED DUE TO ERROR FROM BACOLR.'
      write(*,9003) 'Failure at T = ', t0
      write(*,9002) 'IDID = ', idid
      write(*,9002) 'NINT = ', nint

      goto 600
c-----------------------------------------------------------------------
c  Invalid input error section
  999 continue
      write(*,*) 'Error: Improperly formatted input'

      stop
c-----------------------------------------------------------------------
c  Formats!
 9000 format( a, i2, a, i4, a, i3, a, e7.1 )
 9001 format( 2(a, e7.1), 2a )
 9002 format( a, i10 )
 9003 format( a, f15.8 )
c-----------------------------------------------------------------------
      end
