c
c ===================================================================
subroutine fmod1(maxmx,mvar,meqn,maux,mbc,mx,q,ql,qr,
& aux,dx,dt,method,cfl,fm,fp,q1,aux1,
& amdq,apdq)
c ===================================================================
c
c # This function applies Larrouturou's positivity fix to the fluctuations
c # of the two species of ZND-Euler equations.
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)
dimension q(mvar, 1-mbc:maxmx+mbc)
dimension ql(mvar, 1-mbc:maxmx+mbc)
dimension qr(mvar, 1-mbc:maxmx+mbc)
dimension fm(mvar, 1-mbc:maxmx+mbc)
dimension fp(mvar, 1-mbc:maxmx+mbc)
dimension aux(maux, 1-mbc:maxmx+mbc)
dimension q1(1-mbc:maxmx+mbc, meqn)
dimension aux1(1-mbc:maxmx+mbc, maux)
dimension amdq(1-mbc:maxmx+mbc, meqn)
dimension apdq(1-mbc:maxmx+mbc, meqn)
dimension method(7)
common /rpnflx/ mrpnflx
c
if (mrpnflx.ne.0) then
write(6,*) '*** Riemann solver returns fluxes.'
write(6,*) '*** Correction is implemented for fluctuations.'
write(6,*) '*** Set method(1)=0.'
stop
endif
c
mu = 3
do 100 i=2-mbc,mx+mbc
c
c # Apply the fix to the Godunov-fluxes of the RECONSTRUCTED values
c
fmrho = fm(1,i)+fm(2,i)+q(mu,i-1)
fprho = fp(1,i)+fp(2,i)+q(mu,i )
rhol = qr(1,i-1)+qr(2,i-1)
rhor = ql(1,i )+ql(2,i )
do 110 m=1,2
if (fmrho.gt.0.d0) then
Y = qr(m,i-1)/rhol
else
Y = ql(m,i )/rhor
endif
fm(m,i) = Y*fmrho
fp(m,i) = Y*fprho
110 continue
c
c # Now subtract flux of the original cell values to obtain the waves
c
rhol = q(1,i-1)+q(2,i-1)
rhor = q(1,i )+q(2,i )
do 120 m=1,2
fm(m,i) = fm(m,i) - q(m,i-1)/rhol*q(mu,i-1)
fp(m,i) = fp(m,i) - q(m,i )/rhor*q(mu,i )
120 continue
100 continue
c
return
end