3d/equations/euler/rp/rpn3euswg.f

c
c
c     ==================================================================
      subroutine rpn3eu(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,
     &			maux,auxl,auxr,wave,s,fl,fr)
c     ==================================================================
c
c     # solve Riemann problems for the 3D 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     # 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 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)
      dimension auxl(1-mbc:maxm+mbc, maux, 3)
      dimension auxr(1-mbc:maxm+mbc, maux, 3)
      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 and mw to the 
c     # orthogonal momentum:
c
      if(ixyz .eq. 1)then
	  mu = 2
	  mv = 3
          mw = 4
      else if(ixyz .eq. 2)then
	  mu = 3
	  mv = 4
          mw = 2
      else
          mu = 4
          mv = 2
          mw = 3
      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
         wl = qr(i-1,mw)/rhol
         wr = ql(i  ,mw)/rhor
	 pl = gamma1*(qr(i-1,5) - 0.5d0*(ul**2+vl**2+wl**2)*rhol)
	 pr = gamma1*(ql(i  ,5) - 0.5d0*(ur**2+vr**2+wr**2)*rhor)
         Hl = (qr(i-1,5)+pl)/rhol
         Hr = (ql(i  ,5)+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,mw) = wl*taul + wr*taur
         fl(i,5)  = 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 mws=1,mwaves
            s(i,mws) = dmax1(dabs(el(mws)),dabs(er(mws)))
            do 10 m=1,meqn
               wave(i,m,mws) = 0.d0
 10   continue
c
      return
      end