c  main program
      common /cstak/ ds
      double precision ds(4000)
      common /time/ t
      real t
      common /param/ vc, x, xi0
      real vc(4), x(3), xi0
      external dee, handle, uofx, bc, af
      integer ndx, istkgt, k, immm, iu, is(1000)
      integer nu, nv, imesh, ilumb, nmesh
      real errpar(2), tstart, d, v(4), dt, xb(3)
      real rs(1000), ws(1000), tstop
      logical ls(1000)
      complex cs(500)
      equivalence (ds(1), cs(1), ws(1), rs(1), is(1), ls(1))
c to test  post 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 = 1e+6
      dt = 1e-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 = ilumb(xb, 3, ndx, k, nmesh)
c make 1d0 of multiplicity k-1.
      imesh = immm(imesh, nmesh, 1e0, k-1)
      x(1) = 0
      x(2) = 2.3*d/tstart
      x(3) = 1
c initial values for v.
      call lplmg(3, x, vc)
c get u on port stack.
      iu = istkgt(nmesh-k, 3)
c uofx needs time.
      t = tstart
      vc(4) = tstart
c uofx needs v for mapping.
      call movefr(nv, vc, v)
c initial conditions for u.
      call l2sff(uofx, k, ws(imesh), nmesh, ws(iu))
c output ics.
      call handle(t-1., ws(iu), v, t, ws(iu), v, nu, nmesh-k, nv, k, ws(
     1   imesh), nmesh, dt, tstop)
      call post(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
      real t, xi(nx), u(nx, nu), ux(nx, nu), ut(nx, nu), utx(nx, nu)
      real v(nv), vt(nv), a(nx, nu), au(nx, nu, nu), aux(nx, nu, nu), 
     1   aut(nx, nu, nu)
      real autx(nx, nu, nu), av(nx, nu, nv), avt(nx, nu, nv), f(nx, nu),
     1   fu(nx, nu, nu), fux(nx, nu, nu)
      real fut(nx, nu, nu), futx(nx, nu, nu), fv(nx, nu, nv), fvt(nx, 
     1   nu, nv)
      common /postf/ failed
      logical failed
      integer i
      real xxi(99), xtv(99), xvv(99), x(99), expl, xxiv(99)
      real ax(99), fx(99), xt(99), xv(99)
      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 lplm(xi, nx, v, 3, x, xxi, xxiv, xv, xvv, xt, xtv)
c map u into x system.
      call postu(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))-expl((-x(i))*v(4))*(x(i)+v(4)**2)
         fut(i, 1, 1) = -1
         fv(i, 1, 4) = (-expl((-x(i))*v(4)))*(2.*v(4)+(x(i)+v(4)**2)*(-x
     1      (i)))
         fx(i) = (-expl((-x(i))*v(4)))*(1.-v(4)*x(i)-v(4)**3)
   2     continue
c map a and f into xi system.
      call posti(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
      real t, l, r, u(nu, 2), ux(nu, 2), ut(nu, 2)
      real utx(nu, 2), v(nv), vt(nv), b(nu, 2), bu(nu, nu, 2), bux(nu, 
     1   nu, 2)
      real but(nu, nu, 2), butx(nu, nu, 2), bv(nu, nv, 2), bvt(nu, nv, 2
     1   )
      real expl
c u(0,t) = 1
      b(1, 1) = u(1, 1)-1.
c u(1,t) = exp(-t)
      b(1, 2) = u(1, 2)-expl(-v(4))
      bu(1, 1, 1) = 1
      bu(1, 1, 2) = 1
      bv(1, 4, 2) = expl(-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
      real t, x(nx), u(nxmk, nu), ut(nxmk, nu), v(nv), vt(nv)
      real d(nv), du(nv, nxmk, nu), dut(nv, nxmk, nu), dv(nv, nv), dvt(
     1   nv, nv)
      common /param/ vc, xc, xi0
      real vc(4), xc(3), xi0
      integer intrvr, i, ileft
      real exp, bx(10), xx(1)
      integer temp
      d(1) = v(1)
c x(0,v) = 0.
      dv(1, 1) = 1
      xx(1) = xi0
      ileft = intrvr(nx, x, xx(1))
c get the b-spline basis at xx.
      call bspln(k, x, nx, xx, 1, ileft, bx)
      d(2) = -exp(-1e0)
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)-1.
c x(2,v) = 1.
      dv(3, 3) = 1
      d(4) = vt(4)-1.
      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
      real t0, u0(nxmk, nu), v0(nv), t, u(nxmk, nu), v(nv)
      real x(nx), dt, tstop
      common /param/ vc, xx, xi0
      real vc(4), xx(3), xi0
      common /time/ tt
      real tt
      external uofx
      integer i1mach
      real eu, eesff, ev, lplmt
      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 = lplmt(t, v, nv, t0, v0, 1e-1, dt)
      tt = t
c uofx needs v for mapping.
      call movefr(nv, v, vc)
      eu = eesff(k, x, nx, u, uofx)
c error in position of boundary layer.
      ev = v(2)-1./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
      real xi(nx), u(nx), w(nx)
      common /cstak/ ds
      double precision ds(500)
      common /param/ vc, x, xi0
      real vc(4), x(3), xi0
      common /time/ t
      real t
      integer ixv, ixx, istkgt, i, is(1000)
      real expl, rs(1000), ws(1000)
      logical ls(1000)
      complex cs(500)
      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, 3)
c space for x and xv.
      ixv = istkgt(3*nx, 3)
c map into user system.
      call lplmx(xi, nx, vc, 3, ws(ixx), ws(ixv))
      do  1 i = 1, nx
         temp = ixx+i
         u(i) = expl((-ws(temp-1))*t)
   1     continue
      call leave
      return
      end