2d/equations/euler/rp/rpn2euausmdvg.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     # an improved version of the Liou-Steffen Flux-Vector-Splitting 
c   
c     # Yasuhiro Wada, Meng-Sing Liou "An accurate and robust flux
c     # splitting scheme for shock and contact discontinuities",
c     # SIAM J. Sci. Comput., Vol. 18, No.2, pp 633-657, May 1997.
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, 
c     #            s the speeds, 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, rp is called with ql = qr = q.
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 l(4), r(4)
      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     # AUSM 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*(qr(i-1,mu)**2+
     &      qr(i-1,mv)**2)/rhol)
	 pr = gamma1*(ql(i  ,4) - 0.5d0*(ql(i  ,mu)**2+
     &      ql(i  ,mv)**2)/rhor)
         Hl = (qr(i-1,4)+pl)/rhol
         Hr = (ql(i  ,4)+pr)/rhor
         al = dsqrt(gamma*pl/rhol)
         ar = dsqrt(gamma*pr/rhor)
c
         am = dmax1(al,ar)
         alphal = 2.d0*(pl/rhol)/(pl/rhol+pr/rhor)
         alphar = 2.d0*(pr/rhor)/(pl/rhol+pr/rhor)
c
         ulp = 0.5d0*(ul+dabs(ul))
         plp = pl*ulp/ul
         if (dabs(ul).le.am) then
            ulp = 0.25d0*alphal*(ul+al)**2/am + (1.d0-alphal)*ulp
            plp = 0.25d0*pl*(ul+al)**2/am**2*(2.d0-ul/am)
         endif
c
         urm = 0.5d0*(ur-dabs(ur))
         prm = pr*urm/ur
         if (dabs(ur).le.am) then
            urm = -0.25d0*alphar*(ur-ar)**2/am + (1.d0-alphar)*urm
            prm =  0.25d0*pr*(ur-ar)**2/am**2*(2.d0+ur/am)
         endif
c
c     #  Blending between AUSMV and AUSMD 
c     #  sf=1.d0 gives AUSMV, sf=-1.d0 gives AUSMD
         sf = dmin1(1.d0, 10.d0*dabs(pr-pl)/dmin1(pl,pr))
c
         l(1)  = 0.5d0*(ulp*rhol+dabs(ulp*rhol))
         l(mu) = 0.5d0*((1.d0+sf)*ulp*rhol*ul + (1.d0-sf)*l(1)*ul) + plp
         l(mv) = l(1)*vl
         l(4)  = l(1)*Hl
c
         r(1)  = 0.5d0*(urm*rhor-dabs(urm*rhor))
         r(mu) = 0.5d0*((1.d0+sf)*urm*rhor*ur + (1.d0-sf)*r(1)*ur) + prm
         r(mv) = r(1)*vr
         r(4)  = r(1)*Hr
c
         do 20 m = 1,meqn
            fl(i,m) = l(m) + r(m)
            fr(i,m) = -fl(i,m)
 20      continue
c
         s(i,1) = 0.5d0*(ul-al + ur-ar)
         s(i,2) = 0.5d0*(ul    + ur)
         s(i,3) = 0.5d0*(ul+al + ur+ar)
         do 10 mw=1,mwaves
            do 10 m=1,meqn
               wave(i,m,mw) = 0.d0
 10   continue
c
      return
      end
c
c