C$TEST  DPOST7
c  main program
      common /cstak/ ds
      double precision ds(4000)
      common /time/ t
      double precision t
      common /param/ vc, x, xi0
      double precision vc(4), x(3), xi0
      external dee, handle, uofx, bc, af
      integer ndx, idlumb, istkgt, k, iu, is(1000)
      integer nu, nv, immmd, imesh, nmesh
      real errpar(2), rs(1000)
      logical ls(1000)
      complex cs(500)
      double precision tstart, d, v(4), dt, xb(3), ws(500)
      double precision tstop
      equivalence (ds(1), cs(1), ws(1), rs(1), is(1), ls(1))
c to test dpost on
c      u sub t = u sub xx + f      on (20,10**6)
c where f is chosen so that the solution is
c      u(x,t) = exp(-x*t),
c and x(1,t) is chosen so that the boundary-layer is tracked
c implicitly by forcing u(x(1,t)/2.3/d,t) = 1/e.
c this is the same as requiring the exact solution to have
c u(x(1,t),t) = 10 ** -d.
c v(1,2,3) gives the moving mesh, v(4) is time.
c the port library stack and its aliases.
c initialize the port library stack length.
      call istkin(4000, 4)
      call enter(1)
      nu = 1
      nv = 4
      errpar(1) = 1e-2
c mixed relative and absolute error.
      errpar(2) = 1e-2
      d = 3
c w(xi0,t) = 1/e.
      xi0 = 1./2.3/d
      tstart = 20
      tstop = 1d+6
      dt = 1d-2
      k = 4
c ndx uniform mesh points on each interval of xb.
      ndx = 6
      xb(1) = 0
      xb(2) = 1
      xb(3) = 2
c get mesh on port stack.
      imesh = idlumb(xb, 3, ndx, k, nmesh)
c make 1d0 of multiplicity k-1.
      imesh = immmd(imesh, nmesh, 1d0, k-1)
      x(1) = 0
      x(2) = 2.3*d/tstart
      x(3) = 1
c initial values for v.
      call dlplmg(3, x, vc)
c get u on port stack.
      iu = istkgt(nmesh-k, 4)
c uofx needs time.
      t = tstart
      vc(4) = tstart
c uofx needs v for mapping.
      call movefd(nv, vc, v)
c initial conditions for u.
      call dl2sff(uofx, k, ws(imesh), nmesh, ws(iu))
c output ics.
      call handle(t-1d0, ws(iu), v, t, ws(iu), v, nu, nmesh-k, nv, k, 
     1   ws(imesh), nmesh, dt, tstop)
      call dpost(ws(iu), nu, k, ws(imesh), nmesh, v, nv, tstart, tstop
     1   , dt, af, bc, dee, errpar, handle)
      call leave
      call wrapup
      stop 
      end
      subroutine af(t, xi, nx, u, ux, ut, utx, nu, v, vt, nv, a, 
     1   au, aux, aut, autx, av, avt, f, fu, fux, fut, futx, fv, fvt)
      integer nu, nv, nx
      double precision t, xi(nx), u(nx, nu), ux(nx, nu), ut(nx, nu), 
     1   utx(nx, nu)
      double precision v(nv), vt(nv), a(nx, nu), au(nx, nu, nu), aux(nx,
     1   nu, nu), aut(nx, nu, nu)
      double precision autx(nx, nu, nu), av(nx, nu, nv), avt(nx, nu, nv)
     1   , f(nx, nu), fu(nx, nu, nu), fux(nx, nu, nu)
      double precision fut(nx, nu, nu), futx(nx, nu, nu), fv(nx, nu, nv)
     1   , fvt(nx, nu, nv)
      common /dpostf/ failed
      logical failed
      integer i
      double precision xxi(99), xtv(99), xvv(99), x(99), xxiv(99), ax(
     1   99)
      double precision fx(99), xt(99), xv(99), dexpl
      logical temp
      temp = v(2) .le. v(1)
      if (.not. temp) temp = v(2) .ge. v(3)
      if (.not. temp) goto 1
         failed = .true.
         return
c map xi into x.
   1  call dlplm(xi, nx, v, 3, x, xxi, xxiv, xv, xvv, xt, xtv)
c map u into x system.
      call dpostu(xi, x, xt, xxi, xv, vt, nx, 3, ux, ut, nu, ax, fx)
      do  2 i = 1, nx
         a(i, 1) = -ux(i, 1)
         aux(i, 1, 1) = -1
         f(i, 1) = (-ut(i, 1))-dexpl((-x(i))*v(4))*(x(i)+v(4)**2)
         fut(i, 1, 1) = -1
         fv(i, 1, 4) = (-dexpl((-x(i))*v(4)))*(2d0*v(4)+(x(i)+v(4)**2)*(
     1      -x(i)))
         fx(i) = (-dexpl((-x(i))*v(4)))*(1d0-v(4)*x(i)-v(4)**3)
   2     continue
c map a and f into xi system.
      call dposti(xi, x, xt, xxi, xv, xtv, xxiv, xvv, nx, ux, ut, nu, v,
     1   vt, nv, 1, 3, a, ax, au, aux, aut, autx, av, avt, f, fx, fu, 
     2   fux, fut, futx, fv, fvt)
      return
      end
      subroutine bc(t, l, r, u, ux, ut, utx, nu, v, vt, nv, b, bu,
     1   bux, but, butx, bv, bvt)
      integer nu, nv
      double precision t, l, r, u(nu, 2), ux(nu, 2), ut(nu, 2)
      double precision utx(nu, 2), v(nv), vt(nv), b(nu, 2), bu(nu, nu, 2
     1   ), bux(nu, nu, 2)
      double precision but(nu, nu, 2), butx(nu, nu, 2), bv(nu, nv, 2), 
     1   bvt(nu, nv, 2)
      double precision dexpl
c u(0,t) = 1
      b(1, 1) = u(1, 1)-1d0
c u(1,t) = exp(-t)
      b(1, 2) = u(1, 2)-dexpl(-v(4))
      bu(1, 1, 1) = 1
      bu(1, 1, 2) = 1
      bv(1, 4, 2) = dexpl(-v(4))
      return
      end
      subroutine dee(t, k, x, nx, u, ut, nu, nxmk, v, vt, nv, d, 
     1   du, dut, dv, dvt)
      integer nxmk, nu, nv, nx
      integer k
      double precision t, x(nx), u(nxmk, nu), ut(nxmk, nu), v(nv), vt(
     1   nv)
      double precision d(nv), du(nv, nxmk, nu), dut(nv, nxmk, nu), dv(
     1   nv, nv), dvt(nv, nv)
      common /param/ vc, xc, xi0
      double precision vc(4), xc(3), xi0
      integer intrvd, i, ileft
      double precision dexp, bx(10), xx(1)
      integer temp
      d(1) = v(1)
c x(0,v) = 0.
      dv(1, 1) = 1
      xx(1) = xi0
      ileft = intrvd(nx, x, xx(1))
c get the b-spline basis at xx.
      call dbspln(k, x, nx, xx, 1, ileft, bx)
      d(2) = -dexp(-1d0)
c d(2) = w(xi0,t) - exp(-1).
      do  1 i = 1, k
         temp = ileft+i-k
         d(2) = d(2)+u(temp, 1)*bx(i)
         temp = ileft+i-k
         du(2, temp, 1) = bx(i)
   1     continue
      d(3) = v(3)-1d0
c x(2,v) = 1.
      dv(3, 3) = 1
      d(4) = vt(4)-1d0
      dvt(4, 4) = 1
      return
      end
      subroutine handle(t0, u0, v0, t, u, v, nu, nxmk, nv, k, x, 
     1   nx, dt, tstop)
      integer nxmk, nu, nv, nx
      integer k
      double precision t0, u0(nxmk, nu), v0(nv), t, u(nxmk, nu), v(nv)
      double precision x(nx), dt, tstop
      common /param/ vc, xx, xi0
      double precision vc(4), xx(3), xi0
      common /time/ tt
      double precision tt
      external uofx
      integer i1mach
      double precision deesff, dlplmt, eu, ev
      integer temp
c output and checking routine.
      if (t0 .ne. t) goto 2
         temp = i1mach(2)
         write (temp,  1) t, dt
   1     format (16h restart for t =, 1pe10.2, 7h   dt =, 1pe10.2)
         return
c let dt carry v(2) down by no more than a factor of 10.
   2  dt = dlplmt(t, v, nv, t0, v0, 1d-1, dt)
      tt = t
c uofx needs v for mapping.
      call movefd(nv, v, vc)
      eu = deesff(k, x, nx, u, uofx)
c error in position of boundary layer.
      ev = v(2)-1d0/xi0/t
      temp = i1mach(2)
      write (temp,  3) t, eu, ev
   3  format (14h error in u(x,, 1pe10.2, 4h ) =, 1pe10.2, 6h   v =, 1p
     1   e10.2)
      return
      end
      subroutine uofx(xi, nx, u, w)
      integer nx
      double precision xi(nx), u(nx), w(nx)
      common /cstak/ ds
      double precision ds(500)
      common /param/ vc, x, xi0
      double precision vc(4), x(3), xi0
      common /time/ t
      double precision t
      integer ixv, ixx, istkgt, i, is(1000)
      real rs(1000)
      logical ls(1000)
      complex cs(500)
      double precision ws(500), dexpl
      integer temp
      equivalence (ds(1), cs(1), ws(1), rs(1), is(1), ls(1))
c the port library stack and its aliases.
      call enter(1)
      ixx = istkgt(nx, 4)
c space for x and xv.
      ixv = istkgt(3*nx, 4)
c map into user system.
      call dlplmx(xi, nx, vc, 3, ws(ixx), ws(ixv))
      do  1 i = 1, nx
         temp = ixx+i
         u(i) = dexpl((-ws(temp-1))*t)
   1     continue
      call leave
      return
      end