c
c
c ===================================================================
subroutine step1(maxmx,mvar,meqn,maux,mwaves,mbc,mx,
& qold,aux,dx,t,dt,method,mthlim,cflgrid,
& fm,fp,wave,s,
& amdq,apdq,dtdx,aux1,q1,work,mwork,rp1)
c ===================================================================
c
c Author: Randall J. LeVeque
c Modified for AMROC: Ralf Deiterding
c
c # Main entry from AMROC.
c # Take one time step, updating qold.
c
c # fm, fp are fluxes to left and right of single cell edge
c
c method(2) = 1 ==> Godunov method
c method(2) = 2 ==> Wave limiter method
c mthlim(p) controls what limiter is used in the pth family
c method(2) = 3 ==> Slope-limiting for the conserved quantities
c method(2) > 3 ==> User defined slope-limiter method
c
c amdq, apdq, wave, s, and f are used locally:
c
c amdq(1-mbc:maxmx+mbc, meqn) = left-going flux-differences
c apdq(1-mbc:maxmx+mbc, meqn) = right-going flux-differences
c e.g. amdq(i,m) = m'th component of A^- \Delta q from i'th Riemann
c problem (between cells i-1 and i).
c
c wave(1-mbc:maxmx+mbc, meqn, mwaves) = waves from solution of
c Riemann problems,
c wave(i,m,mw) = mth component of jump in q across
c wave in family mw in Riemann problem between
c states i-1 and i.
c
c s(1-mbc:maxmx+mbc, mwaves) = wave speeds,
c s(i,mw) = speed of wave in family mw in Riemann problem between
c states i-1 and i.
c
c f(1-mbc:maxmx+mbc, meqn) = correction fluxes for second order method
c f(i,m) = mth component of flux at left edge of ith cell
c --------------------------------------------------------------------
c
implicit double precision (a-h,o-z)
include "call.i"
c
dimension qold(mvar, 1-mbc:maxmx+mbc)
dimension fm(mvar, 1-mbc:maxmx+mbc)
dimension fp(mvar, 1-mbc:maxmx+mbc)
dimension wave(1-mbc:maxmx+mbc, meqn, mwaves)
dimension s(1-mbc:maxmx+mbc, mwaves)
dimension amdq(1-mbc:maxmx+mbc, meqn)
dimension apdq(1-mbc:maxmx+mbc, meqn)
dimension aux(maux, 1-mbc:maxmx+mbc)
dimension dtdx(1-mbc:maxmx+mbc)
dimension aux1(1-mbc:maxmx+mbc, maux)
dimension q1(1-mbc:maxmx+mbc, meqn)
dimension method(7),mthlim(mwaves)
dimension work(mwork)
external rp1
common /rpnflx/ mrpnflx
c
i0ql = 1
i0qr = i0ql + (maxmx+2*mbc)*meqn
i0fl = i0qr + (maxmx+2*mbc)*meqn
i0fr = i0fl + (maxmx+2*mbc)*meqn
if (method(1).eq.1) then
i0qls = i0fr + (maxmx+2*mbc)*meqn
i0qrs = i0qls + (maxmx+2*mbc)*mvar
iused = i0qrs + (maxmx+2*mbc)*mvar - 1
else
iused = i0fr + (maxmx+2*mbc)*meqn - 1
endif
c
if (iused.gt.mwork) then
write(6,*) '*** not enough work space in step1'
write(6,*) '*** iused = ', iused, ' mwork =',mwork
stop
endif
c
tcom = t
dtcom = dt
dxcom = dx
c
mcapa = method(6)
c
c # check if any limiters are used:
limit = 0
do 5 mw=1,mwaves
if (mthlim(mw) .gt. 0) limit = 1
5 continue
c
do 10 i=1-mbc,mx+mbc
if (mcapa.gt.0) then
if (aux(mcapa,i) .le. 0.d0) then
write(6,*) 'Error -- capa must be positive'
stop
endif
dtdx(i) = dt / (dx*aux(mcapa,i))
else
dtdx(i) = dt/dx
endif
10 continue
do 20 i = 1-mbc, mx+mbc
do 20 m=1,meqn
q1(i,m) = qold(m,i)
20 continue
if (maux .gt. 0) then
do 22 i = 1-mbc, mx+mbc
do 22 ma=1,maux
aux1(i,ma) = aux(ma,i)
22 continue
endif
do 23 i = 1-mbc, mx+mbc
do 23 m=1,mvar
fm(m,i) = 0.d0
fp(m,i) = 0.d0
23 continue
c
c
c
c # solve Riemann problem at each interface
c -----------------------------------------
c
if (method(2).le.2) then
call rp1(maxmx,meqn,mwaves,mbc,mx,q1,q1,maux,aux1,aux1,
& wave,s,amdq,apdq)
else
call slope1(maxmx,meqn,maux,mwaves,mbc,mx,q1,aux1,dx,dt,
& method,mthlim,wave,s,amdq,apdq,dtdx,rp1,
& work(i0ql),work(i0qr),work(i0fl),work(i0fr))
if (method(1).eq.1)
& call saverec1(maxmx,mvar,meqn,mbc,mx,qold,work(i0qls),
& work(i0qrs),work(i0ql),work(i0qr))
endif
c
do 40 i=1,mx+1
do 40 m=1,meqn
fp(m,i) = fp(m,i) - apdq(i,m)
fm(m,i) = fm(m,i) + amdq(i,m)
40 continue
c
c # compute maximum wave speed:
cflgrid = 0.d0
do 50 mw=1,mwaves
do 50 i=1,mx+1
cflgrid = dmax1(cflgrid, dtdx(i)*dabs(s(i,mw)))
50 continue
c
if (method(2).eq.1.or.method(2).ge.3) go to 900
c
if (mrpnflx.ne.0) then
write(6,*) '*** Riemann solver returns fluxes.'
write(6,*) '*** Wave limiting not possible.'
write(6,*) '*** Set method(2)>=3 for slope limiting.'
stop
endif
c
c # compute correction fluxes for second order q_{xx} terms:
c ----------------------------------------------------------
c
c # apply limiter to waves:
if (limit.eq.1) call limiter(maxmx,meqn,mwaves,mbc,mx,
& wave,s,mthlim)
c
do 120 i=1,mx+1
dtdxave = 0.5d0 * (dtdx(i-1) + dtdx(i))
do 120 m=1,meqn
cqxx = 0.d0
do 110 mw=1,mwaves
c # second order corrections:
cqxx = cqxx + dabs(s(i,mw))
& * (1.d0 - dabs(s(i,mw))*dtdxave) * wave(i,m,mw)
110 continue
fm(m,i) = fm(m,i) + 0.5d0 * cqxx
fp(m,i) = fp(m,i) + 0.5d0 * cqxx
c
120 continue
c
900 continue
c
if (method(1).eq.1) then
if (method(2).le.2) then
call fmod1(maxmx,mvar,meqn,maux,mbc,mx,qold,qold,qold,
& aux,dx,dt,method,cflgrid,fm,fp,q1,aux1,
& amdq,apdq)
else
call fmod1(maxmx,mvar,meqn,maux,mbc,mx,qold,work(i0qls),
& work(i0qrs),aux,dx,dt,method,cflgrid,fm,fp,
& q1,aux1,amdq,apdq)
endif
endif
c
c # update q
c ============================================
c
do 150 i=1,mx
do 150 m=1,meqn
qold(m,i) = qold(m,i) - dtdx(i) * (fm(m,i+1) - fp(m,i))
150 continue
c
c ============================================
c
return
end
c
c
c ===================================================================
subroutine slope1(maxm,meqn,maux,mwaves,mbc,mx,
& q,aux,dx,dt,method,mthlim,
& wave,s,amdq,apdq,dtdx,rp1,
& ql,qr,fl,fr)
c ===================================================================
c
c # Implements the standard MUSCL-Hancock method, which gives 2nd
c # order accuracy for conventional finite-volume schemes that
c # return the numerical fluxes instead of fluctuations.
c
c # Author: Ralf Deiterding, ralf@amroc.net
c
implicit double precision (a-h,o-z)
common /rpnflx/ mrpnflx
c
dimension q(1-mbc:maxm+mbc, meqn)
dimension ql(1-mbc:maxm+mbc, meqn)
dimension qr(1-mbc:maxm+mbc, meqn)
dimension aux(1-mbc:maxm+mbc, maux)
dimension fl(1-mbc:maxm+mbc, meqn)
dimension fr(1-mbc:maxm+mbc, meqn)
dimension wave(1-mbc:maxm+mbc, meqn, mwaves)
dimension s(1-mbc:maxm+mbc, mwaves)
dimension amdq(1-mbc:maxm+mbc, meqn)
dimension apdq(1-mbc:maxm+mbc, meqn)
dimension dtdx(1-mbc:maxm+mbc)
dimension method(7),mthlim(mwaves)
external rp1
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
do 100 m=1,meqn
ql(1-mbc,m) = q(1-mbc,m)
qr(1-mbc,m) = q(1-mbc,m)
ql(mx+mbc,m) = q(mx+mbc,m)
qr(mx+mbc,m) = q(mx+mbc,m)
100 continue
c
if (method(2).gt.3) then
call rec1(maxm,meqn,mwaves,mbc,mx,q,method,mthlim,ql,qr)
c
c # MUSCL reconstruction with slope-limiting for conserved
c # variables. Linear reconstruction: om=0.d0
c # Quadratic spatial reconstuction: om!=0.d0
c # 2nd order spatial reconstruction: om=1.d0/3.d0
c # Reconstructions with om!=0.d0 are not conservative!
c
else
om = 0.d0
do 110 i=2-mbc,mx+mbc-1
do 110 m=1,meqn
call reclim(q(i,m),q(i-1,m),q(i+1,m),
& mlim,om,ql(i,m),qr(i,m))
110 continue
endif
c
call flx1(maxm,meqn,mbc,mx,ql,maux,aux,fl)
call flx1(maxm,meqn,mbc,mx,qr,maux,aux,fr)
c
do 200 i=2-mbc,mx+mbc-1
do 200 m=1,meqn
ql(i,m) = ql(i,m) + 0.5d0*dtdx(i)*(fl(i,m)-fr(i,m))
qr(i,m) = qr(i,m) + 0.5d0*dtdx(i)*(fl(i,m)-fr(i,m))
200 continue
c
call rp1(maxm,meqn,mwaves,mbc,mx,ql,qr,maux,aux,aux,
& wave,s,amdq,apdq)
c
c # Add differences of fluxes between original and reconstructed values
c # to fluctuations, if a wave propagation scheme is applied.
c
if (mrpnflx.eq.0) then
call flx1(maxm,meqn,mbc,mx,ql,maux,aux,fl)
call flx1(maxm,meqn,mbc,mx,qr,maux,aux,fr)
do 300 i=2-mbc,mx+mbc-1
do 300 m=1,meqn
amdq(i,m) = amdq(i,m) + fr(i-1,m)
apdq(i,m) = apdq(i,m) - fl(i ,m)
300 continue
call flx1(maxm,meqn,mbc,mx,q,maux,aux,fl)
do 310 i=2-mbc,mx+mbc-1
do 310 m=1,meqn
amdq(i,m) = amdq(i,m) - fl(i-1,m)
apdq(i,m) = apdq(i,m) + fl(i ,m)
310 continue
endif
c
return
end
c
c
c
c ===================================================================
subroutine saverec1(maxm,mvar,meqn,mbc,mx,q,qls,qrs,ql,qr)
c ===================================================================
c
c # Store reconstructed values for later use.
c
implicit double precision (a-h,o-z)
dimension q(mvar, 1-mbc:maxm+mbc)
dimension qls(mvar, 1-mbc:maxm+mbc)
dimension qrs(mvar, 1-mbc:maxm+mbc)
dimension ql(1-mbc:maxm+mbc, meqn)
dimension qr(1-mbc:maxm+mbc, meqn)
c
do 10 i = 1-mbc, mx+mbc
do m=1,meqn
qls(m,i) = ql(i,m)
qrs(m,i) = qr(i,m)
enddo
do m=meqn+1,mvar
qls(m,i) = q(m,i)
qrs(m,i) = q(m,i)
enddo
10 continue
c
return
end
c