c c =================================================================== subroutine rec3(ixyz,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) 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 mu = 3 mv = 4 mw = 5 mE = 6 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, 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, the sum of the c # 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=mu,mE 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 c # ---------------------------------------------------------------- else 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) - q(i ,2)*q0 - & 0.5d0*rho *(u**2 +v**2 +w**2)) p1 = gamma1*(q(i-1,mE) - q(i-1,2)*q0 - & 0.5d0*rho1*(u1**2+v1**2+w1**2)) p2 = gamma1*(q(i+1,mE) - q(i+1,2)*q0 - & 0.5d0*rho2*(u2**2+v2**2+w2**2)) c 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,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 + ql(i,2)*q0 + & 0.5d0*rhol*(ul**2+vl**2+wl**2) qr(i,mE) = pr/gamma1 + qr(i,2)*q0 + & 0.5d0*rhor*(ur**2+vr**2+wr**2) endif c 110 continue c return end c c