2d/equations/euler/rp/rec2eu.f

c
c ===================================================================
      subroutine rec2(ixy,maxm,meqn,mwaves,mbc,mx,q,method,mthlim,ql,qr)
c ===================================================================
c
c     # Reconstruction in primitive variables for the two-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
      mE = 4
c
      mlim = 0
      do 90 mw=1,mwaves
	 if (mthlim(mw) .gt. 0) then
            mlim = mthlim(mw)
            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
         p  = gamma1*(q(i  ,mE) - 0.5d0*rho *(u**2 +v**2 ))
         p1 = gamma1*(q(i-1,mE) - 0.5d0*rho1*(u1**2+v1**2))
         p2 = gamma1*(q(i+1,mE) - 0.5d0*rho2*(u2**2+v2**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(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,mE) = pl/gamma1+0.5d0*rhol*(ul**2+vl**2)
         qr(i,mE) = pr/gamma1+0.5d0*rhor*(ur**2+vr**2)
c     
 110  continue
c     
      return
      end
c
c