c     ==================================================================
      subroutine rpn3  (ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,
     &			auxl,auxr,maux,wave,s,amdq,apdq)
c     ==================================================================
c
c     # Roe-solver for the Euler equations
c      -----------------------------------------------------------
c
c     # solve Riemann problems along one slice of data.
c     # This data is along a slice in the x-direction if ixyz=1
c     #                               the y-direction if ixyz=2.
c     #                               the z-direction if ixyz=3.
c
c     # On input, ql contains the state vector at the left edge of each cell
c     #           qr contains the state vector at the right edge of each cell
c
c     # On output, wave contains the waves, s the speeds,
c     # and amdq, apdq the left-going and right-going flux differences,
c     # respectively.
c
c     # Note that the i'th Riemann problem has left state qr(i-1,:)
c     #                                    and right state ql(i,:)
c     # From the basic clawpack routines, this routine is called with ql = qr
c
c
      implicit none

      integer ixyz, maxm, meqn, mwaves, mbc, mx, maux
c
      double precision wave(1-mbc:maxm+mbc, meqn, mwaves)
      double precision    s(1-mbc:maxm+mbc, mwaves)
      double precision   ql(1-mbc:maxm+mbc, meqn)
      double precision   qr(1-mbc:maxm+mbc, meqn)
      double precision amdq(1-mbc:maxm+mbc, meqn)
      double precision apdq(1-mbc:maxm+mbc, meqn)
      double precision auxl(1-mbc:maxm+mbc, maux)
      double precision auxr(1-mbc:maxm+mbc, maux)

      double precision ql_proj(5),qr_proj(5), s_rot(3), wave_rot(5,3)
      double precision amdq_rot(5), apdq_rot(5), rot(9)
      double precision uvw(3), delta(5), rhol, rhor, ul, ur, vl, vr
      double precision wl, wr, el, er,pl, pr, u2v2w2l, u2v2w2r
      double precision rhsqrtl, rhsqrtr, rhsq2, enth
      double precision speeds(3,2), area, cl, cr

      integer mcapa, locrot, locarea, i, j,m,mws

      logical efix

      double precision dtcom, dxcom, dycom, dzcom, tcom
      double precision gamma, gamma1
      integer icom, jcom, kcom

      common /comxyzt/ dtcom,dxcom,dycom,dzcom,tcom,icom,jcom,kcom
      common /param/ gamma, gamma1
c
      data efix /.false./    !# use entropy fix for transonic rarefactions
c
      if (mwaves .ne. 3) then
         write(6,*) '*** Should have mwaves=3 for this Riemann solver'
         stop
      endif

      call get_aux_locations_n(ixyz,mcapa,locrot,locarea)

      do i = 2-mbc, mx+mbc

c        # Get  two stored vectors; third is computed when rotations are done
         do j = 1,6
            rot(j) = auxl(i,locrot+j-1)
         enddo
         call compute_binormal(rot)

         do m = 1,meqn
            qr_proj(m) = ql(i,m)
            ql_proj(m) = qr(i-1,m)
         enddo

         call rotate3(rot,qr_proj(2))
         call rotate3(rot,ql_proj(2))

         rhol = ql_proj(1)
         rhor = qr_proj(1)

         ul = ql_proj(2)/rhol
         ur = qr_proj(2)/rhor

         vl = ql_proj(3)/rhol
         vr = qr_proj(3)/rhor

         wl = ql_proj(4)/rhol
         wr = qr_proj(4)/rhor

         el = ql_proj(5)
         er = qr_proj(5)

         u2v2w2l = ul*ul + vl*vl + wl*wl
         u2v2w2r = ur*ur + vr*vr + wr*wr
         pl = gamma1*(el - 0.5d0*rhol*u2v2w2l)
         pr = gamma1*(er - 0.5d0*rhor*u2v2w2r)


c        # Use Roe averaged values for normal solves
c        # Get Roe averaged values
         rhsqrtl = sqrt(rhol)
         rhsqrtr = sqrt(rhor)
         rhsq2 = rhsqrtl + rhsqrtr

         uvw(1) = (ul*rhsqrtl + ur*rhsqrtr) / rhsq2
         uvw(2) = (vl*rhsqrtl + vr*rhsqrtr) / rhsq2
         uvw(3) = (wl*rhsqrtl + wr*rhsqrtr) / rhsq2
         enth = ((el + pl)/rhsqrtl + (er + pr)/rhsqrtr)/rhsq2

         do j = 1,meqn
            delta(j) = qr_proj(j) - ql_proj(j)
         enddo

c        # Solve normal Riemann problem
         call solve_riemann(uvw, enth, delta, wave_rot, s_rot)

         do mws = 1,mwaves
            speeds(mws,1) = min(s_rot(mws),0.d0)
            speeds(mws,2) = max(s_rot(mws),0.d0)
         enddo

         if (efix) then
c           # This modifies the speeds, but we will still have
c           # s(mw) = speeds(mw,1) + speeds(mw,2)
            cl = sqrt(gamma*pl/rhol)
            cr = sqrt(gamma*pr/rhor)
            call apply_entropy_fix(ql_proj,qr_proj,cl,cr,
     &            wave_rot, speeds)
         endif

         area = auxl(i,locarea)
         do mws = 1,mwaves
            call rotate3_tr(rot,wave_rot(2,mws))
            speeds(mws,1) = area*speeds(mws,1)
            speeds(mws,2) = area*speeds(mws,2)
            s(i,mws) = speeds(mws,1) + speeds(mws,2)
         enddo

         do m = 1,meqn
            amdq(i,m) = 0.d0
            apdq(i,m) = 0.d0
            do mws = 1,mwaves
               wave(i,m,mws) = wave_rot(m,mws)
               amdq(i,m) = amdq(i,m) + speeds(mws,1)*wave(i,m,mws)
               apdq(i,m) = apdq(i,m) + speeds(mws,2)*wave(i,m,mws)
            enddo
         enddo

      enddo  !! end of i loop over 1d sweep array

      return
      end


      subroutine solve_riemann(uvw,enth,delta,wave_rot,s)
      implicit none

      double precision enth, uvw(3), delta(5)
      double precision wave_rot(5,3), s(3)
      double precision u2v2w2, a2, a, g1a2, euv
      double precision a1, a3, a4, a5,u,v,w
      double precision gamma, gamma1
      logical debug

      common /param/ gamma, gamma1


      u = uvw(1)
      v = uvw(2)
      w = uvw(3)

      u2v2w2 = u**2 + v**2 + w**2
      a2 = gamma1*(enth - 0.5d0*u2v2w2)
      if (a2 .lt. 0.d0) then
         write(6,*) 'a2 .lt. 0; ', a2
c         stop
      endif
      a = dsqrt(a2)
      g1a2 = gamma1 / a2
      euv = enth - u2v2w2

      a4 = g1a2 * (euv*delta(1)
     &      + u*delta(2) + v*delta(3) + w*delta(4)
     &      - delta(5))
      a2 = delta(3) - v*delta(1)
      a3 = delta(4) - w*delta(1)
      a5 = (delta(2) + (a-u)*delta(1) - a*a4) / (2.d0*a)
      a1 = delta(1) - a4 - a5
c
c     # Compute the waves.
c     # Note that the 2-wave, 3-wave and 4-wave travel at the same speed
c     # and are lumped together in wave(.,.,2).  The 5-wave is then stored
c     # in wave(.,.,3).
c
      wave_rot(1,1)  = a1
      wave_rot(2,1) = a1*(u-a)
      wave_rot(3,1) = a1*v
      wave_rot(4,1) = a1*w
      wave_rot(5,1)  = a1*(enth - u*a)
      s(1) = u - a
c
      wave_rot(1,2)  = a4
      wave_rot(2,2) = a4*u
      wave_rot(3,2) = a4*v + a2
      wave_rot(4,2) = a4*w + a3
      wave_rot(5,2)  = a4*0.5d0*u2v2w2  + a2*v + a3*w
      s(2) = u
c
      wave_rot(1,3)  = a5
      wave_rot(2,3) = a5*(u+a)
      wave_rot(3,3) = a5*v
      wave_rot(4,3) = a5*w
      wave_rot(5,3)  = a5*(enth+u*a)
      s(3) = u + a

      end

      subroutine apply_entropy_fix(ql,qr,cl, cr, wave_local,speeds)
      implicit none

      double precision ql(5), qr(5),wave_local(5,3)
      double precision speeds(3,2), s1, s2, s3
      double precision sl, sml, smr, sr, qml(5),qmr(5)
      double precision ul, cl, ur, cr, pml, pmr, cml,cmr
      double precision sfract
      double precision gamma, gamma1
      logical trans1

      integer m

      common /param/ gamma, gamma1

      s1 = speeds(1,1) + speeds(1,2)
      s2 = speeds(2,1) + speeds(2,2)
      s3 = speeds(3,1) + speeds(3,2)

      do m = 1,5
         qml(m) = ql(m) + wave_local(m,1)
      enddo
      sl = ql(2)/ql(1) - cl
      pml = gamma1*(qml(5) - 0.5d0*(qml(2)**2/qml(1) +
     &      qml(3)**2/qml(1) + qml(4)**2/qml(1)))
      sml = qml(2)/qml(1) - sqrt(gamma*pml/qml(1))

c     # Check the 1-wave
      trans1 = .false.
      if (sl .lt. 0 .and. sml .gt. 0) then
c        # apply transonic entropy fix
         trans1 = .true.
         s1 = (sl + sml)/2.d0
         sfract = (sml - s1)/(sml - sl)
         speeds(1,1) = sfract*sl
         speeds(1,2) = (1-sfract)*sml
      endif

c     # If the 1-wave is transonic,then we are done...
c     # Otherwise, we have to check the 3-wave
      if (.not. trans1) then
         do m = 1,5
            qmr(m) = qr(m) - wave_local(m,3)
         enddo
         sr = qr(2)/qr(1) + cr
         pmr = gamma1*(qmr(5) - 0.5d0*(qmr(2)**2/qmr(1) +
     &         qmr(3)**2/qmr(1) + qmr(4)**2/qmr(1)))
         smr = qmr(2)/qmr(1) + sqrt(gamma*pmr/qmr(1))
         if (smr .lt. 0 .and. sr .gt. 0) then
c           # apply transonic entropy fix
            s3 = (smr + sr)/2.d0
            sfract = (sr - s3)/(sr - smr)
            speeds(3,1) = sfract*smr
            speeds(3,2) = (1-sfract)*sr
         endif
      endif

      end
