c  main program
      common /cstak/ ds
      double precision ds(4000)
      common /time/ t
      real t
      common /param/ vc, x
      real vc(4), x(3)
      external dee, handle, uofx, bc, af
      integer ndx, istkgt, k, immm, iu, is(1000)
      integer nu, nv, imesh, ilumb, nmesh
      real errpar(2), tstart, v(4), dt, xb(3), rs(1000)
      real 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 the hyperbolic problem
c      u sub t = - u sub x + g      on (-pi,+pi) * (0,pi)
c with a moving shock x(t) characterized by
c    u(x(t)+,t) = 0 and
c    u(x(t)+,t) - u(x(t)-,t) = x'(t)
c where g is chosen so that the solution is
c               sin(x+t)  for x < x(t)
c      u(x,t) = 
c               cos(x+t)  for x > x(t)
c with x(t) = pi/2 -t .
c v(1,2,3) gives the moving mesh and v(4) is the height of the jump.
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) = 0
c absolute error.
      errpar(2) = 1e-2
      tstart = 0
      tstop = 3.14
      dt = 0.4
      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 1 of multiplicity k-1.
      imesh = immm(imesh, nmesh, 1e0, k-1)
      x(1) = -3.14
      x(2) = 3.14/2.
      x(3) = 3.14
c initial values for v.
      call lplmg(3, x, vc)
c get u on the port stack.
      iu = istkgt(nmesh-k, 3)
c uofx needs time.
      t = tstart
c the initial height of the jump.
      vc(4) = 1
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 cos, sin, xxi(99), xtv(99), xvv(99), x(99)
      real xxiv(99), 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  4 i = 1, nx
         a(i, 1) = -u(i, 1)
         au(i, 1, 1) = -1
         f(i, 1) = ut(i, 1)
         fut(i, 1, 1) = 1
         if (xi(i) .gt. 1.) goto 2
            f(i, 1) = f(i, 1)-2.*cos(x(i)+t)
            fx(i) = 2.*sin(x(i)+t)
            goto  3
   2        f(i, 1) = f(i, 1)-vt(4)
            fvt(i, 1, 4) = -1
            f(i, 1) = f(i, 1)+2.*sin(x(i)+t)
            fx(i) = 2.*cos(x(i)+t)
   3     continue
   4     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 sin
      b(1, 1) = u(1, 1)-sin(t-3.14)
c u(-pi,t) = sin(-pi+t).
      bu(1, 1, 1) = 1
      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)
      integer intrvr, i, ileft
      real bx(10), xx(1), r1mach
      integer temp
      d(1) = v(1)+3.14
c x(0,v) = -pi.
      dv(1, 1) = 1
c xx(1) = 1 + a rounding error.
      xx(1) = r1mach(4)+1.
      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) = -v(4)
c u(x(t)+,t) - jump = 0.
      dv(2, 4) = -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)-3.14
c x(2,v) = +pi.
      dv(3, 3) = 1
c jump + d( x(1,v(t)) )/dt = 0.
      d(4) = vt(2)+v(4)
      dvt(4, 2) = 1
      dv(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
      real vc(4), xx(3)
      common /time/ tt
      real tt
      external uofx
      integer i1mach
      real eu, eesff, ev(2)
      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
   2  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 shock.
      ev(1) = v(2)-(3.14/2.-t)
c error in height of shock.
      ev(2) = v(4)-1.
      temp = i1mach(2)
      write (temp,  3) t, eu, ev
   3  format (14h error in u(x,, 1pe10.2, 4h ) =, 1pe10.2, 6h   v =, 2(
     1   1pe10.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
      real vc(4), x(3)
      common /time/ t
      real t
      integer ixv, ixx, istkgt, i, is(1000)
      real ewe, 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) = ewe(t, ws(temp-1), vc(2))
         if (xi(i) .gt. 1.) u(i) = u(i)+1.
   1     continue
      call leave
      return
      end
      real function ewe(t, x, xbreak)
      real t, x, xbreak
      real cos, sin
      if (x .ge. xbreak) goto 1
         ewe = sin(x+t)
         return
   1     if (x .le. xbreak) goto 2
            ewe = cos(x+t)
            return
   2        call seterr(17hewe - x == xbreak, 17, 1, 2)
   3  continue
   4  stop 
      end