c
c
c ==========================================================
subroutine unsp2(maxm,maxmx,maxmy,mvar,meqn,maux,mwaves,mbc,
& mx,my,qold,aux,dx,dy,dt,method,mthlim,cflgrid,
& fm,fp,gm,gp,faddm,faddp,gaddm,gaddp,
& q1d,dtdx1d,dtdy1d,aux1,aux2,aux3,
& work,mwork,rpn2,rpt2)
c ==========================================================
c
c # Take one time step, updating q with the wave propagation method
c # of Randall J. LeVeque.
c
c # fm, fp are fluxes to left and right of single cell edge,
c # gm, gp are fluxes at bottom and top.
c # fadd and gadd are used to return flux increments from flux2.
c # See the flux2 documentation for more information.
c
implicit double precision (a-h,o-z)
include "call.i"
c
dimension qold(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
dimension fm(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
dimension fp(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
dimension gm(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
dimension gp(mvar, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
dimension q1d(1-mbc:maxm+mbc, meqn)
dimension faddm(1-mbc:maxm+mbc, meqn)
dimension faddp(1-mbc:maxm+mbc, meqn)
dimension gaddm(1-mbc:maxm+mbc, meqn, 2)
dimension gaddp(1-mbc:maxm+mbc, meqn, 2)
dimension aux(maux, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
dimension aux1(1-mbc:maxm+mbc, maux)
dimension aux2(1-mbc:maxm+mbc, maux)
dimension aux3(1-mbc:maxm+mbc, maux)
dimension dtdx1d(1-mbc:maxm+mbc)
dimension dtdy1d(1-mbc:maxm+mbc)
dimension method(7),mthlim(mwaves)
dimension work(mwork)
external rpn2,rpt2
c
c # partition work array into pieces needed for local storage in
c # flux2 routine. Find starting index of each piece:
c
i0wave = 1
i0s = i0wave + (maxm+2*mbc)*meqn*mwaves
i0amdq = i0s + (maxm+2*mbc)*mwaves
i0apdq = i0amdq + (maxm+2*mbc)*meqn
i0cqxx = i0apdq + (maxm+2*mbc)*meqn
i0bmadq = i0cqxx + (maxm+2*mbc)*meqn
i0bpadq = i0bmadq + (maxm+2*mbc)*meqn
if (method(1).ge.1) then
i0qls = i0bpadq + (maxm+2*mbc)*meqn
i0qrs = i0qls + (maxmx+2*mbc)*(maxmy+2*mbc)*mvar
i0qbs = i0qrs + (maxmx+2*mbc)*(maxmy+2*mbc)*mvar
i0qts = i0qbs + (maxmx+2*mbc)*(maxmy+2*mbc)*mvar
i0ql = i0qts + (maxmx+2*mbc)*(maxmy+2*mbc)*mvar
i0qr = i0ql + (maxm+2*mbc)*meqn
i0slope = i0ql
else
i0slope = i0bpadq + (maxm+2*mbc)*meqn
endif
iused = i0slope + (maxm+2*mbc)*meqn*4 - 1
mslope = iused-i0slope+1
c
if (iused.gt.mwork) then
write(6,*) '*** not enough work space in unsp2'
write(6,*) '*** iused = ', iused, ' mwork =',mwork
stop
endif
c
mcapa = method(6)
c
cflgrid = 0.d0
dtdx = dt/dx
dtdy = dt/dy
c
do 10 j=1-mbc,my+mbc
do 10 i=1-mbc,mx+mbc
do 10 m=1,mvar
fm(m,i,j) = 0.d0
fp(m,i,j) = 0.d0
gm(m,i,j) = 0.d0
gp(m,i,j) = 0.d0
10 continue
c
if (mcapa.eq.0) then
c # no capa array:
do 5 i=1-mbc,maxm+mbc
dtdx1d(i) = dtdx
dtdy1d(i) = dtdy
5 continue
endif
c
c
c # perform x-sweeps
c ==================
c
do 50 j = 0,my+1
c
c # copy data along a slice into 1d arrays:
do 20 i = 1-mbc, mx+mbc
do 20 m=1,meqn
q1d(i,m) = qold(m,i,j)
20 continue
c
if (mcapa.gt.0) then
do 21 i = 1-mbc, mx+mbc
dtdx1d(i) = dtdx / aux(mcapa,i,j)
21 continue
endif
c
if (maux .gt. 0) then
do 22 i = 1-mbc, mx+mbc
do 22 ma=1,maux
aux1(i,ma) = aux(ma,i,j-1)
aux2(i,ma) = aux(ma,i,j)
aux3(i,ma) = aux(ma,i,j+1)
22 continue
endif
c
c
c # Store the value of j along this slice in the common block
c # comxyt in case it is needed in the Riemann solver (for
c # variable coefficient problems)
jcom = j
c
c # compute modifications fadd and gadd to fluxes along this slice:
call flux2(1,maxm,meqn,maux,mwaves,mbc,mx,
& q1d,dtdx1d,aux1,aux2,aux3,method,mthlim,
& faddm,faddp,gaddm,gaddp,cfl1d,
& work(i0wave),work(i0s),work(i0amdq),work(i0apdq),
& work(i0cqxx),work(i0bmadq),work(i0bpadq),
& work(i0slope),mslope,rpn2,rpt2)
c
cflgrid = dmax1(cflgrid,cfl1d)
c
c # update fluxes for use in AMR:
c
do 25 i=1,mx+1
do 25 m=1,meqn
fm(m,i,j) = fm(m,i,j) + faddm(i,m)
fp(m,i,j) = fp(m,i,j) + faddp(i,m)
if (method(3).gt.0) then
gm(m,i,j) = gm(m,i,j) + gaddm(i,m,1)
gp(m,i,j) = gp(m,i,j) + gaddp(i,m,1)
gm(m,i,j+1) = gm(m,i,j+1) + gaddm(i,m,2)
gp(m,i,j+1) = gp(m,i,j+1) + gaddp(i,m,2)
endif
25 continue
c
if (method(1).ge.1.and.method(2).ge.3)
& call saverec2(1,maxm,maxmx,maxmy,mvar,meqn,mbc,mx,my,
& qold,work(i0qls),work(i0qrs),work(i0qbs),
& work(i0qts),work(i0ql),work(i0qr))
c
50 continue
c
c
c
c # perform y sweeps
c ==================
c
c
do 100 i = 0, mx+1
c
c # copy data along a slice into 1d arrays:
do 70 j = 1-mbc, my+mbc
do 70 m=1,meqn
q1d(j,m) = qold(m,i,j)
70 continue
c
if (mcapa.gt.0) then
do 71 j = 1-mbc, my+mbc
dtdy1d(j) = dtdy / aux(mcapa,i,j)
71 continue
endif
c
if (maux .gt. 0) then
do 72 j = 1-mbc, my+mbc
do 72 ma=1,maux
aux1(j,ma) = aux(ma,i-1,j)
aux2(j,ma) = aux(ma,i, j)
aux3(j,ma) = aux(ma,i+1,j)
72 continue
endif
c
c
c # Store the value of i along this slice in the common block
c # comxyt in case it is needed in the Riemann solver (for
c # variable coefficient problems)
icom = i
c
c # compute modifications fadd and gadd to fluxes along this slice:
call flux2(2,maxm,meqn,maux,mwaves,mbc,my,
& q1d,dtdy1d,aux1,aux2,aux3,method,mthlim,
& faddm,faddp,gaddm,gaddp,cfl1d,
& work(i0wave),work(i0s),work(i0amdq),work(i0apdq),
& work(i0cqxx),work(i0bmadq),work(i0bpadq),
& work(i0slope),mslope,rpn2,rpt2)
c
cflgrid = dmax1(cflgrid,cfl1d)
c
c #
c # update fluxes for use in AMR:
c
do 75 j=1,my+1
do 75 m=1,meqn
gm(m,i,j) = gm(m,i,j) + faddm(j,m)
gp(m,i,j) = gp(m,i,j) + faddp(j,m)
if (method(3).gt.0) then
fm(m,i,j) = fm(m,i,j) + gaddm(j,m,1)
fp(m,i,j) = fp(m,i,j) + gaddp(j,m,1)
fm(m,i+1,j) = fm(m,i+1,j) + gaddm(j,m,2)
fp(m,i+1,j) = fp(m,i+1,j) + gaddp(j,m,2)
endif
75 continue
c
if (method(1).ge.1.and.method(2).ge.3)
& call saverec2(2,maxm,maxmx,maxmy,mvar,meqn,mbc,mx,my,
& qold,work(i0qls),work(i0qrs),work(i0qbs),
& work(i0qts),work(i0ql),work(i0qr))
c
100 continue
c
if (method(1).ge.1) then
if (method(2).le.2) then
call fmod2(maxm,maxmx,maxmy,mvar,meqn,maux,mbc,mx,my,
& qold,qold,qold,qold,qold,aux,dx,dy,dt,method,
& cflgrid,fm,fp,gm,gp,q1d,aux2,faddm,faddp,0)
else
call fmod2(maxm,maxmx,maxmy,mvar,meqn,maux,mbc,mx,my,
& qold,work(i0qls),work(i0qrs),work(i0qbs),
& work(i0qts),aux,dx,dy,dt,method,cflgrid,
& fm,fp,gm,gp,q1d,aux2,faddm,faddp,0)
endif
endif
c
do 110 j=1,my
do 110 i=1,mx
do 110 m=1,meqn
if (mcapa.gt.0) then
c # with capa array.
qold(m,i,j) = qold(m,i,j)
& - (dtdx * (fm(m,i+1,j) - fp(m,i,j))
& + dtdy * (gm(m,i,j+1) - gp(m,i,j))) /
& aux(mcapa,i,j)
c # no capa array. Standard flux differencing:
else
qold(m,i,j) = qold(m,i,j)
& - dtdx * (fm(m,i+1,j) - fp(m,i,j))
& - dtdy * (gm(m,i,j+1) - gp(m,i,j))
endif
110 continue
c
return
end
c