c
c
c     =====================================================
      subroutine rpt2(ixy,maxm,meqn,mwaves,mbc,mx,
     &                  ql,qr,aux1,aux2,aux3,
     &                  ilr,asdq,bmasdq,bpasdq)
c     =====================================================
      implicit double precision (a-h,o-z)
c
c     # Riemann solver in the transverse direction for the Euler equations
c     #  with a tracer variable.
c     # Split asdq (= A^* \Delta q, where * = + or -)
c     # into down-going flux difference bmasdq (= B^- A^* \Delta q)
c     #    and up-going flux difference bpasdq (= B^+ A^* \Delta q)
c
c     # Uses Roe averages and other quantities which were 
c     # computed in rpn2eu and stored in the common block comroe.
c
      dimension     ql(1-mbc:maxm+mbc, meqn)
      dimension     qr(1-mbc:maxm+mbc, meqn)
      dimension   asdq(1-mbc:maxm+mbc, meqn)
      dimension bmasdq(1-mbc:maxm+mbc, meqn)
      dimension bpasdq(1-mbc:maxm+mbc, meqn)
c
      common /param/  gamma,gamma1
      dimension waveb(5,4),sb(4)
      parameter (maxm2 = 802)  !# assumes at most 600x600 grid with mbc=2
      common /comroe/ u2v2(-1:maxm2),
     &       u(-1:maxm2),v(-1:maxm2),enth(-1:maxm2),a(-1:maxm2),
     &       g1a2(-1:maxm2),euv(-1:maxm2) 
c
      if (-1.gt.1-mbc .or. maxm2 .lt. maxm+mbc) then
         write(6,*) 'need to increase maxm2 in rpt'
         stop
         endif
c
      if (ixy.eq.1) then
          mu = 2
          mv = 3
        else
          mu = 3
          mv = 2
        endif
c
         do 20 i = 2-mbc, mx+mbc
            a3 = g1a2(i) * (euv(i)*asdq(i,1) 
     &             + u(i)*asdq(i,mu) + v(i)*asdq(i,mv) - asdq(i,4))
            a2 = asdq(i,mu) - u(i)*asdq(i,1)
            a4 = (asdq(i,mv) + (a(i)-v(i))*asdq(i,1) - a(i)*a3)
     &              / (2.d0*a(i))
            a1 = asdq(i,1) - a3 - a4
c
            waveb(1,1) = a1
            waveb(mu,1) = a1*u(i)
            waveb(mv,1) = a1*(v(i)-a(i))
            waveb(4,1) = a1*(enth(i) - v(i)*a(i))
            waveb(5,1) = 0.d0
            sb(1) = v(i) - a(i)
c
            waveb(1,2) = a3
            waveb(mu,2) = a3*u(i) + a2
            waveb(mv,2) = a3*v(i)
            waveb(4,2) = a3*0.5d0*u2v2(i) + a2*u(i)
            waveb(5,2) = 0.d0
            sb(2) = v(i)
c
            waveb(1,3) = a4
            waveb(mu,3) = a4*u(i)
            waveb(mv,3) = a4*(v(i)+a(i))
            waveb(4,3) = a4*(enth(i)+v(i)*a(i))
            waveb(5,3) = 0.d0
            sb(3) = v(i) + a(i)
c
            waveb(1,4) = 0.d0
            waveb(mu,4) = 0.d0
            waveb(mv,4) = 0.d0
            waveb(4,4) = 0.d0
            waveb(5,4) = asdq(i,5)
            sb(4) = v(i)
c
c           # compute the flux differences bmasdq and bpasdq
c
            do 10 m=1,meqn
               bmasdq(i,m) = 0.d0
               bpasdq(i,m) = 0.d0
               do 10 mw=1,4
                  bmasdq(i,m) = bmasdq(i,m) 
     &                         + dmin1(sb(mw), 0.d0) * waveb(m,mw)
                  bpasdq(i,m) = bpasdq(i,m)
     &                         + dmax1(sb(mw), 0.d0) * waveb(m,mw)
   10             continue
c                 
   20       continue
c
      return
      end