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

c
c
c     =====================================================
      subroutine rpn3eu(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr,
     &                  maux,auxl,auxr,wave,s,fl,fr)
c     =====================================================
c
c     # FORCE scheme for the 3D Euler equations. The flux of the FORCE
c     # scheme is the arithmetic mean of the fluxes of the finite difference
c     # schemes of Richtmyer and Lax-Friedrichs. Use parameters
c     # richtmyer, laxfriedrich to switch to the original schemes.
c
c     # Eleuterio F. Toro, "Riemann solvers and numerical methods
c     # for fluid dynamics", Springer-Verlag, Berlin 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     # 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 step1, 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)
      common /param/  gamma,gamma1
      include "call.i"
c
c     # local storage
c     ---------------
      parameter (max2 = 10000)  !# assumes at most 10000 grid points with mbc=2
      dimension qint(-1:max2+2,5), fint(-1:max2+2,5), 
     &     auxint(-1:max2+2,0,3)
      logical richtmyer, laxfriedrich
c
      data richtmyer    /.true./     
      data laxfriedrich /.true./     
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
      dxdt = 0.5d0*dxcom/dtcom
      dtdx = 0.5d0*dtcom/dxcom
c
      call flx3(ixyz,maxm,meqn,mbc,mx,ql,maux,auxl,fl)
      call flx3(ixyz,maxm,meqn,mbc,mx,qr,maux,auxr,fr)
c
      do 50 i = 2-mbc, mx+mbc
         do 50 m=1,meqn
            qint(i,m) = 0.5d0*(qr(i-1,m) + ql(i,m)) + 
     &           dtdx*(fr(i-1,m) - fl(i,m))
 50   continue
      do 60 i = 2-mbc, mx+mbc
         do 60 m=1,maux
            auxint(i,m,2) = 0.5d0*(auxl(i,m,2) + auxr(i,m,2)) 
 60   continue
      call flx3(ixyz,max2,meqn,mbc,mx,qint,maux,auxint,fint)
c
      do 100 i = 2-mbc, mx+mbc
         ul = 0.5d0*qr(i-1,mu)/qr(i-1,1)
         ur = 0.5d0*ql(i  ,mu)/ql(i  ,1)
         pl = gamma1*(qr(i-1,5) - 0.5d0*(qr(i-1,mu)**2+
     &        qr(i-1,mv)**2+qr(i-1,mw)**2)/qr(i-1,1))
         pr = gamma1*(ql(i  ,5) - 0.5d0*(ql(i  ,mu)**2+
     &        ql(i  ,mv)**2+ql(i  ,mw)**2)/ql(i  ,1))
         al = dsqrt(gamma*pl/qr(i-1,1))
         ar = dsqrt(gamma*pr/ql(i  ,1))
         s(i,1) = dmax1(dabs(ul-al),dabs(ur-ar))
         s(i,2) = dmax1(dabs(ul   ),dabs(ur   ))
         s(i,3) = dmax1(dabs(ul+al),dabs(ur+ar))
         do 110 mws=1,mwaves
            do 110 m=1,meqn
               wave(i,m,mws) = 0.d0
 110     continue
         do 100 m=1,meqn
            if (richtmyer) 
     &           fl(i,m) = fint(i,m)
            if (laxfriedrich) 
     &           fl(i,m) = dxdt*(qr(i-1,m) - ql(i,m)) + 
     &           0.5d0*(fr(i-1,m) + fl(i,m))
            if (richtmyer.and.laxfriedrich)
     &           fl(i,m) = 0.5d0*(fl(i,m) + fint(i,m))
 100  continue
c
      do 120 i = 2-mbc, mx+mbc
         do 120 m=1,meqn
            fr(i,m) = -fl(i,m)
 120  continue
c
      return
      end
c