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

c
c ===================================================================
      subroutine rec3(ixyz,maxm,meqn,mwaves,mbc,mx,
     &     q,method,mthlim,ql,qr)
c ===================================================================
c
c     # Reconstruction in primitive variables for the three-dimensional
c     # Euler equations. This reconstruction is not conservative!
c 
      implicit double precision (a-h,o-z)
c
      dimension    q(1-mbc:maxm+mbc, meqn)
      dimension   ql(1-mbc:maxm+mbc, meqn)
      dimension   qr(1-mbc:maxm+mbc, meqn)
      dimension method(7),mthlim(mwaves)
      common /param/  gamma,gamma1
c
      mu = 2
      mv = 3
      mw = 4
      mE = 5
c
      mlim = 0
      do 90 mws=1,mwaves
	 if (mthlim(mws) .gt. 0) then
            mlim = mthlim(mws)
            goto 95
         endif
 90   continue
 95   continue
c
c     Linear interpolation: om=0.d0, 2nd order: om!=0.d0,
c     3rd order: om=1.d0/3.d0
c
      om = 1.d0/3.d0     
      do 110 i=2-mbc,mx+mbc-1
c
c     # Reconstruction
         rho  = q(i  ,1)
         rho1 = q(i-1,1)
         rho2 = q(i+1,1)   
         u  = q(i  ,mu)/rho
         u1 = q(i-1,mu)/rho1
         u2 = q(i+1,mu)/rho2
         v  = q(i  ,mv)/rho
         v1 = q(i-1,mv)/rho1
         v2 = q(i+1,mv)/rho2
         w  = q(i  ,mw)/rho
         w1 = q(i-1,mw)/rho1
         w2 = q(i+1,mw)/rho2
         p  = gamma1*(q(i  ,mE) - 0.5d0*rho *(u**2 +v**2 +w**2 ))
         p1 = gamma1*(q(i-1,mE) - 0.5d0*rho1*(u1**2+v1**2+w1**2))
         p2 = gamma1*(q(i+1,mE) - 0.5d0*rho2*(u2**2+v2**2+w2**2))
c    
         call reclim(rho,rho1,rho2,mlim,om,rhol,rhor)  
         call reclim(u,u1,u2,mlim,om,ul,ur)  
         call reclim(v,v1,v2,mlim,om,vl,vr)  
         call reclim(w,w1,w2,mlim,om,wl,wr)  
         call reclim(p,p1,p2,mlim,om,pl,pr)  
c     
         ql(i,1)  = rhol
         qr(i,1)  = rhor
         ql(i,mu) = ul*rhol
         qr(i,mu) = ur*rhor
         ql(i,mv) = vl*rhol
         qr(i,mv) = vr*rhor
         ql(i,mw) = wl*rhol
         qr(i,mw) = wr*rhor
         ql(i,mE) = pl/gamma1+0.5d0*rhol*(ul**2+vl**2+wl**2)
         qr(i,mE) = pr/gamma1+0.5d0*rhor*(ur**2+vr**2+wr**2)
c     
 110  continue
c     
      return
      end
c
c