c
c ===================================================================
subroutine rec1(maxm,meqn,mwaves,mbc,mx,q,method,mthlim,ql,qr)
c ===================================================================
c
c # Copyright (C) 2002 Ralf Deiterding
c # Brandenburgische Universitaet Cottbus
c
c # Copyright (C) 2003-2007 California Institute of Technology
c # Ralf Deiterding, ralf@amroc.net
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),
& Yk1(2), Yk2(2), Ykl(2), Ykr(2)
common /param/ gamma,gamma1,q0
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, quadratic interpolation: om!=0.d0,
c # Second order accuracte reconstruction for om=1.d0/3.d0
c
om = 0.d0
do 110 i=2-mbc,mx+mbc-1
c
c # Reconstruction of total density
c
rho = q(i ,1) + q(i ,2)
rho1 = q(i-1,1) + q(i-1,2)
rho2 = q(i+1,1) + q(i+1,2)
call reclim(rho,rho1,rho2,mlim,om,rhol,rhor)
c
c # Reconstruction of mass fractions - if the same limiter value
c # is choosen for left and right side and for all mass fractions,
c # the sum of the reconstructed mass fractions is 1.
c
sl = 1.d0
do k = 1, 2
Yk1(k) = q(i ,k)/rho - q(i-1,k)/rho1
Yk2(k) = q(i+1,k)/rho2 - q(i ,k)/rho
sl = dmin1(sl,slopelim(Yk1(k),Yk2(k),mlim))
sl = dmin1(sl,slopelim(Yk2(k),Yk1(k),mlim))
enddo
c
do k = 1, 2
Ykl(k) = q(i,k)/rho - 0.25d0*((1.d0+om)*sl*Yk1(k) +
& (1.d0-om)*sl*Yk2(k))
Ykr(k) = q(i,k)/rho + 0.25d0*((1.d0-om)*sl*Yk1(k) +
& (1.d0+om)*sl*Yk2(k))
ql(i,k) = Ykl(k)*rhol
qr(i,k) = Ykr(k)*rhor
enddo
c
c # Reconstruction in conservative variables
c # ----------------------------------------------------------------
if (method(2).eq.4) then
do m=3,4
call reclim(q(i,m),q(i-1,m),q(i+1,m),
& mlim,om,ql(i,m),qr(i,m))
enddo
c
c # Reconstruction in primitive variables - not conservative
c # gives worse results thanmethod(2)=4 !
c # ----------------------------------------------------------------
else if (method(2).eq.5) then
u = q(i ,3)/rho
u1 = q(i-1,3)/rho1
u2 = q(i+1,3)/rho2
p = gamma1*(q(i ,4) - q(i ,2)*q0 - 0.5d0*rho*u**2)
p1 = gamma1*(q(i-1,4) - q(i-1,2)*q0 - 0.5d0*rho1*u1**2)
p2 = gamma1*(q(i+1,4) - q(i+1,2)*q0 - 0.5d0*rho2*u2**2)
c
call reclim(u,u1,u2,mlim,om,ul,ur)
call reclim(p,p1,p2,mlim,om,pl,pr)
c
ql(i,3) = ul*rhol
qr(i,3) = ur*rhor
ql(i,4) = pl/gamma1+ql(i,2)*q0+0.5d0*rhol*ul**2
qr(i,4) = pr/gamma1+qr(i,2)*q0+0.5d0*rhor*ur**2
else if (method(2).eq.6) then
u = q(i ,3)/rho
u1 = q(i-1,3)/rho1
u2 = q(i+1,3)/rho2
p = gamma1*(q(i ,4) - q(i ,2)*q0 - 0.5d0*rho*u**2)
p1 = gamma1*(q(i-1,4) - q(i-1,2)*q0 - 0.5d0*rho1*u1**2)
p2 = gamma1*(q(i+1,4) - q(i+1,2)*q0 - 0.5d0*rho2*u2**2)
c
call reclim(p,p1,p2,mlim,om,pl,pr)
c
call reclim(q(i,3),q(i-1,3),q(i+1,3),
& mlim,om,ql(i,3),qr(i,3))
ql(i,4) = pl/gamma1+ql(i,2)*q0+0.5d0*ql(i,3)**2/rhol
qr(i,4) = pr/gamma1+qr(i,2)*q0+0.5d0*qr(i,3)**2/rhor
else
u = q(i ,3)/rho
u1 = q(i-1,3)/rho1
u2 = q(i+1,3)/rho2
p = gamma1*(q(i ,4) - q(i ,2)*q0 - 0.5d0*rho*u**2)
p1 = gamma1*(q(i-1,4) - q(i-1,2)*q0 - 0.5d0*rho1*u1**2)
p2 = gamma1*(q(i+1,4) - q(i+1,2)*q0 - 0.5d0*rho2*u2**2)
c
call reclim(p/rho,p1/rho1,p2/rho2,mlim,om,Tl,Tr)
c
call reclim(q(i,3),q(i-1,3),q(i+1,3),
& mlim,om,ql(i,3),qr(i,3))
ql(i,4) = Tl*rhol/gamma1+ql(i,2)*q0+0.5d0*ql(i,3)**2/rhol
qr(i,4) = Tr*rhor/gamma1+qr(i,2)*q0+0.5d0*qr(i,3)**2/rhor
endif
c if ((ql(i,4).lt.q(i-1,4).and.ql(i,4).lt.q(i+1,4)
c & .and.ql(i,4).lt.q(i,4)).or.
c & (ql(i,4).gt.q(i-1,4).and.ql(i,4).gt.q(i+1,4)
c & .and.ql(i,4).gt.q(i,4)).or.
c & (qr(i,4).lt.q(i-1,4).and.qr(i,4).lt.q(i+1,4)
c & .and.qr(i,4).lt.q(i,4)).or.
c & (qr(i,4).gt.q(i-1,4).and.qr(i,4).gt.q(i+1,4)
c & .and.qr(i,4).gt.q(i,4))) write (6,*) 'TVD viol.', i
c
110 continue
c
return
end
c
c