3d/equations/euler/rp/rpn3euvijag.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 ZND-Euler equations using 
c     # the Flux-Vector-Splitting of Vijayasundaram
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 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
      do 10 i=2-mbc,mx+mbc
         rhol  = qr(i-1,1)
         rhor  = ql(i  ,1)
         rhoul = qr(i-1,mu)
         rhour = ql(i  ,mu)
         rhovl = qr(i-1,mv)
         rhovr = ql(i  ,mv)
         rhowl = qr(i-1,mw)
         rhowr = ql(i  ,mw)
         rhoEl = qr(i-1,5)
         rhoEr = ql(i  ,5)
c
         rho  = 0.5d0*(rhol  + rhor )
         rhou = 0.5d0*(rhoul + rhour)
         rhov = 0.5d0*(rhovl + rhovr)
         rhow = 0.5d0*(rhowl + rhowr)
         rhoE = 0.5d0*(rhoEl + rhoEr)
c
         u = rhou/rho
         v = rhov/rho
         w = rhow/rho
	 p = gamma1*(rhoE - 0.5d0*rho*(u**2+v**2+w**2))
         H = (rhoE+p)/rho
         if (p.le.0.d0.or.rho.le.0.d0) 
     &        write (6,*) 'Error in middle state in',i,p,pl,pr,
     &        rho,rhol,rhor,a,al,ar
         a = dsqrt(gamma*p/rho)
         f = 0.5d0/a**2
c
         el(1) = 0.5d0*(u-a + dabs(u-a))
         el(2) = 0.5d0*(u   + dabs(u)  )
         el(3) = 0.5d0*(u+a + dabs(u+a))
         er(1) = 0.5d0*(u-a - dabs(u-a))
         er(2) = 0.5d0*(u   - dabs(u)  )
         er(3) = 0.5d0*(u+a - dabs(u+a))
c
         zl = el(1)-el(3)
         zr = er(1)-er(3)
         ol = el(1)-2.d0*el(2)+el(3)
         or = er(1)-2.d0*er(2)+er(3)
         dul = a*(rhol*u-rhoul)
         dur = a*(rhor*u-rhour)
         dEl = gamma1*(rhoEl+0.5d0*rhol*(u**2+v**2+w**2)-
     &                 rhoul*u-rhovl*v-rhowl*w)
         dEr = gamma1*(rhoEr+0.5d0*rhor*(u**2+v**2+w**2)-
     &                 rhour*u-rhovr*v-rhowr*w)
         f1 =   f*(zl*dul + ol*dEl + zr*dur + or*dEr)
         f2 = a*f*(ol*dul + zl*dEl + or*dur + zr*dEr)
c
         fl(i,1)  = rhol *el(2) + rhor *er(2) + f1
         fl(i,mu) = rhoul*el(2) + rhour*er(2) + u*f1 -   f2
         fl(i,mv) = rhovl*el(2) + rhovr*er(2) + v*f1
         fl(i,mw) = rhowl*el(2) + rhowr*er(2) + w*f1
         fl(i,5)  = rhoEl*el(2) + rhoEr*er(2) + H*f1 - u*f2
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
c