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@cacr.caltech.edu
c
implicit double precision (a-h,o-z)
include "ck.i"
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(LeNsp), Yk2(LeNsp), Ykl(LeNsp), Ykr(LeNsp)
c
mu = Nsp+1
mE = Nsp+2
mT = Nsp+3
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 = 0.d0
rho1 = 0.d0
rho2 = 0.d0
do k = 1, Nsp
rho = rho + q(i ,k)
rho1 = rho1 + q(i-1,k)
rho2 = rho2 + q(i+1,k)
enddo
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, Nsp
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, Nsp
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=Nsp+1,mE
call reclim(q(i,m),q(i-1,m),q(i+1,m),
& mlim,om,ql(i,m),qr(i,m))
enddo
c # No reconstruction for temperature - is updated in flx1eurhok
ql(i,mT) = q(i,mT)
qr(i,mT) = q(i,mT)
c
c # Reconstruction in primitive variables
c # ----------------------------------------------------------------
else if (method(2).eq.5) then
rhoW = 0.d0
rhoW1 = 0.d0
rhoW2 = 0.d0
do k = 1, Nsp
rhoW = rhoW + q(i ,k)/Wk(k)
rhoW1 = rhoW1 + q(i-1,k)/Wk(k)
rhoW2 = rhoW2 + q(i+1,k)/Wk(k)
enddo
c
call reclim(q(i,mu)/rho,q(i-1,mu)/rho1,
& q(i+1,mu)/rho2,mlim,om,ul,ur)
call reclim(rhoW*RU*q(i,mT),rhoW1*RU*q(i-1,mT),
& rhoW2*RU*q(i+1,mT),mlim,om,pl,pr)
ql(i,mu) = ul*rhol
qr(i,mu) = ur*rhor
c
rhoWl = 0.d0
rhoWr = 0.d0
do k = 1, Nsp
rhoWl = rhoWl + ql(i,k)/Wk(k)
rhoWr = rhoWr + qr(i,k)/Wk(k)
enddo
ql(i,mT) = pl/(rhoWl*RU)
qr(i,mT) = pr/(rhoWr*RU)
ql(i,mE) = rhol * (0.5d0*ul**2 +
& avgtabip(ql(i,mT),Ykl,hms,Nsp)) - pl
qr(i,mE) = rhor * (0.5d0*ur**2 +
& avgtabip(qr(i,mT),Ykr,hms,Nsp)) - pr
else if (method(2).eq.6) then
call reclim(q(i,Nsp+1),q(i-1,Nsp+1),q(i+1,Nsp+1),
& mlim,om,ql(i,Nsp+1),qr(i,Nsp+1))
c
rhoW = 0.d0
rhoW1 = 0.d0
rhoW2 = 0.d0
do k = 1, Nsp
rhoW = rhoW + q(i ,k)/Wk(k)
rhoW1 = rhoW1 + q(i-1,k)/Wk(k)
rhoW2 = rhoW2 + q(i+1,k)/Wk(k)
enddo
c
call reclim(rhoW*RU*q(i,mT),rhoW1*RU*q(i-1,mT),
& rhoW2*RU*q(i+1,mT),mlim,om,pl,pr)
c
rhoWl = 0.d0
rhoWr = 0.d0
do k = 1, Nsp
rhoWl = rhoWl + ql(i,k)/Wk(k)
rhoWr = rhoWr + qr(i,k)/Wk(k)
enddo
ql(i,mT) = pl/(rhoWl*RU)
qr(i,mT) = pr/(rhoWr*RU)
ql(i,mE) = rhol * avgtabip(ql(i,mT),Ykl,hms,Nsp) +
& 0.5d0*ql(i,Nsp+1)**2/rhol - pl
qr(i,mE) = rhor * avgtabip(qr(i,mT),Ykr,hms,Nsp) +
& 0.5d0*qr(i,Nsp+1)**2/rhor - pr
else if (method(2).eq.7) then
call reclim(q(i,Nsp+1),q(i-1,Nsp+1),q(i+1,Nsp+1),
& mlim,om,ql(i,Nsp+1),qr(i,Nsp+1))
c
rhoW = 0.d0
rhoW1 = 0.d0
rhoW2 = 0.d0
do k = 1, Nsp
rhoW = rhoW + q(i ,k)/Wk(k)
rhoW1 = rhoW1 + q(i-1,k)/Wk(k)
rhoW2 = rhoW2 + q(i+1,k)/Wk(k)
enddo
c
call reclim(q(i,mT),q(i-1,mT),q(i+1,mT),
& mlim,om,ql(i,mT),qr(i,mT))
c
rhoWl = 0.d0
rhoWr = 0.d0
do k = 1, Nsp
rhoWl = rhoWl + ql(i,k)/Wk(k)
rhoWr = rhoWr + qr(i,k)/Wk(k)
enddo
ql(i,mE) = rhol * avgtabip(ql(i,mT),Ykl,hms,Nsp) +
& 0.5d0*ql(i,Nsp+1)**2/rhol - ql(i,mT)*(rhoWl*RU)
qr(i,mE) = rhor * avgtabip(qr(i,mT),Ykr,hms,Nsp) +
& 0.5d0*qr(i,Nsp+1)**2/rhor - qr(i,mT)*(rhoWr*RU)
endif
c
110 continue
c
return
end
c