1d/equations/euler/rp/rp1euswg.f

c
c
c =========================================================
      subroutine rp1eu(maxmx,meqn,mwaves,mbc,mx,ql,qr,maux,
     &     auxl,auxr,wave,s,fl,fr)
c =========================================================
c
c     # solve Riemann problems for the 1D 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     # 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 routine step1, rp is called with ql = qr = q.
c
c     Author:  Ralf Deiterding
c
      implicit double precision (a-h,l,o-z)
      dimension   ql(1-mbc:maxmx+mbc, meqn)
      dimension   qr(1-mbc:maxmx+mbc, meqn)
      dimension    s(1-mbc:maxmx+mbc, mwaves)
      dimension wave(1-mbc:maxmx+mbc, meqn, mwaves)
      dimension   fl(1-mbc:maxmx+mbc, meqn)
      dimension   fr(1-mbc:maxmx+mbc, meqn)
      double precision el(3), er(3)
      common /param/  gamma,gamma1
c
c     # Method returns fluxes
c     ------------
      common /rpnflx/ mrpnflx
      mrpnflx = 1
c
      do 10 i=2-mbc,mx+mbc
         rhol = qr(i-1,1)
         rhor = ql(i  ,1)
         ul = qr(i-1,2)/rhol
         ur = ql(i  ,2)/rhor
	 pl = gamma1*(qr(i-1,3) - 0.5d0*ul**2*rhol)
	 pr = gamma1*(ql(i  ,3) - 0.5d0*ur**2*rhor)
         Hl = (qr(i-1,3)+pl)/rhol
         Hr = (ql(i  ,3)+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*rhol/gamma
         facr = 0.5d0*rhor/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,2) = ul*taul - zetal + ur*taur - zetar
         fl(i,3) = 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
c