2d/equations/euler/rp/rpn2euswg.f

c
c
c     =====================================================
      subroutine rpn2eu(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr,maux,
     &     auxl,auxr,wave,s,fl,fr)
c     =====================================================
c
c     # solve Riemann problems for the 2D Euler equations using 
c     # Steger & Warming - Flux Vector Splitting 
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     # This data is along a slice in the x-direction if ixy=1 
c     #                            or the y-direction if ixy=2.
c     # On output, wave contains the waves, s the speeds, 
c     # fl and fr the positive and negative flux.
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 routines, this routine is called with ql = qr
c
c     Author: Ralf Deiterding
c
      implicit double precision (a-h,o-z)
      dimension wave(1-mbc:maxm+mbc, meqn, mwaves)
      dimension    s(1-mbc:maxm+mbc, mwaves)
      dimension   ql(1-mbc:maxm+mbc, meqn)
      dimension   qr(1-mbc:maxm+mbc, meqn)
      dimension   fl(1-mbc:maxm+mbc, meqn)
      dimension   fr(1-mbc:maxm+mbc, meqn)
      double precision el(3), er(3)
      common /param/  gamma,gamma1
c
c     # Method returns fluxes
c     ------------
      common /rpnflx/ mrpnflx
      mrpnflx = 1
c
c     # set mu to point to  the component of the system that corresponds
c     # to momentum in the direction of this slice, mv to the orthogonal
c     # momentum:
c
      if (ixy.eq.1) then
         mu = 2
         mv = 3
      else
         mu = 3
         mv = 2
      endif
c
c     #  Steger & Warming - Flux Vector Splitting 
c
      do 10 i=2-mbc,mx+mbc
         rhol = qr(i-1,1)
         rhor = ql(i  ,1)
         ul = qr(i-1,mu)/rhol
         ur = ql(i  ,mu)/rhor
         vl = qr(i-1,mv)/rhol
         vr = ql(i  ,mv)/rhor
	 pl = gamma1*(qr(i-1,4) - 0.5d0*(ul**2+vl**2)*rhol)
	 pr = gamma1*(ql(i  ,4) - 0.5d0*(ur**2+vr**2)*rhor)
         Hl = (qr(i-1,4)+pl)/rhol
         Hr = (ql(i  ,4)+pr)/rhor
c
         al2 = gamma*pl/rhol
         al  = dsqrt(al2)
         ar2 = gamma*pr/rhor
         ar  = dsqrt(ar2)
c
         el(1) = 0.5d0*(ul-al + dabs(ul-al))
         el(2) = 0.5d0*(ul    + dabs(ul)   )
         el(3) = 0.5d0*(ul+al + dabs(ul+al))
         er(1) = 0.5d0*(ur-ar - dabs(ur-ar))
         er(2) = 0.5d0*(ur    - dabs(ur)   )
         er(3) = 0.5d0*(ur+ar - dabs(ur+ar))
c
         facl = 0.5d0*qr(i-1,1)/gamma
         facr = 0.5d0*ql(i  ,1)/gamma
c
         taul  = facl*(el(1) + 2.d0*gamma1*el(2) + el(3))
         taur  = facr*(er(1) + 2.d0*gamma1*er(2) + er(3))
         zetal = al*facl*(el(1)-el(3)) 
         zetar = ar*facr*(er(1)-er(3)) 
c
         fl(i,1)  = taul + taur
         fl(i,mu) = ul*taul - zetal + ur*taur - zetar
         fl(i,mv) = vl*taul + vr*taur
         fl(i,4)  = Hl*taul - ul*zetal - 2.d0*el(2)*facl*al2 + 
     &              Hr*taur - ur*zetar - 2.d0*er(2)*facr*ar2
c
         do 20 m = 1, meqn
            fr(i,m) = -fl(i,m)
 20      continue
c
         do 10 mw=1,mwaves
            s(i,mw) = dmax1(dabs(el(mw)),dabs(er(mw)))
            do 10 m=1,meqn
               wave(i,m,mw) = 0.d0
 10   continue
c
      return
      end