3d/equations/euler/rp/rpn3euvlg.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     # van Leer's 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 routine, 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 Ml, Mr, sl(3), sr(3), fvl(5), fvr(5)
      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 momentums:
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     # Van Leer's Flux Vector Splitting
c
      gamma2 = gamma**2-1
      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
         El = qr(i-1,5)/rhol
         Er = ql(i  ,5)/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)
         al = dsqrt(gamma*pl/rhol)
         ar = dsqrt(gamma*pr/rhor)
c
         Ml = ul/al
         Mr = ur/ar
c
         sl(1) = ul-al
         sl(2) = ul
         sl(3) = ul+al
         sr(1) = ur-ar
         sr(2) = ur
         sr(3) = ur+ar
c
         if (Ml.gt.1d0) then
            fvl(1)  = rhol*ul
            fvl(mu) = fvl(1)*ul+pl
            fvl(mv) = fvl(1)*vl
            fvl(mw) = fvl(1)*wl
            fvl(5)  = ul*(rhol*El+pl)
         else if (Ml.lt.-1.d0) then
            do m = 1,meqn
               fvl(m) = 0.d0
            enddo
         else
            fvl(1)  = 0.25d0*rhol*al*(Ml+1.d0)**2
            fvl(mu) = fvl(1)*2.d0*al/gamma*(0.5d0*gamma1*Ml+1.d0)
            fvl(mv) = fvl(1)*vl
            fvl(mw) = fvl(1)*wl
            fvl(5)  = fvl(1)*(0.5d0*(vl**2+wl**2) + 2.d0*al**2/gamma2*
     &                        (0.5d0*gamma1*Ml+1.d0)**2)
         endif
c
         if (Mr.lt.-1.d0) then
            fvr(1)  = rhor*ur
            fvr(mu) = fvr(1)*ur+pr
            fvr(mv) = fvr(1)*vr
            fvr(mw) = fvr(1)*wr
            fvr(5)  = ur*(rhor*Er+pr)
         else if (Mr.gt.1.d0) then
            do m = 1,meqn
               fvr(m) = 0.d0
            enddo
         else
            fvr(1)  = -0.25d0*rhor*ar*(Mr-1.d0)**2
            fvr(mu) = fvr(1)*2.d0*ar/gamma*(0.5d0*gamma1*Mr-1.d0)
            fvr(mv) = fvr(1)*vr
            fvr(mw) = fvr(1)*wr
            fvr(5)  = fvr(1)*(0.5d0*(vr**2+wr**2) + 2.d0*ar**2/gamma2* 
     &                        (0.5d0*gamma1*Mr-1.d0)**2)
         endif
c
         do 20 m = 1,meqn
            fl(i,m) = fvl(m) + fvr(m)
            fr(i,m) = -fl(i,m)
 20      continue
c
         if (dabs(Ml).lt.1.d0) then
            facl = (gamma+3.d0)/(2.d0*gamma+dabs(Ml)*(3.d0-gamma))
         else
            facl = 1.d0
         endif
         if (dabs(Mr).lt.1.d0) then
            facr = (gamma+3.d0)/(2.d0*gamma+dabs(Mr)*(3.d0-gamma))
         else
            facr = 1.d0
         endif
c
         do 10 mws=1,mwaves
            s(i,mws) = dmax1(dabs(facl*sl(mws)),dabs(facr*sr(mws)))
            do 10 m=1,meqn
               wave(i,m,mws) = 0.d0
 10   continue
c
      return
      end
c