c
c
c =========================================================
subroutine src(maxmx,maxmy,meqn,mbc,ibx,iby,mx,my,q,
& aux,maux,t,dt,ibnd)
c =========================================================
implicit double precision(a-h,o-z)
c
dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc)
c
c # source terms for cylindrical symmetry in 2d Euler equations
c # about x-axis (so y=radius)
c
dimension qstar(4)
common /param/ gamma,gamma1
c
c # 2-stage Runge-Kutta method
c
dt2 = dt/2.d0
press = 0.d0
ndim = 2
do 10 i=1-ibx*ibnd*mbc,mx+ibx*ibnd*mbc
do 10 j=1-iby*ibnd*mbc,my+iby*ibnd*mbc
y = aux(1,i,j)
rho = q(1,i,j)
u = q(2,i,j)/q(1,i,j)
v = q(3,i,j)/q(1,i,j)
press = gamma1*(q(4,i,j) - 0.5d0*rho*(u**2 + v**2))
if (y.eq.0.d0) write(6,*) 'y = 0 in src'
if (rho.eq.0.d0) write(6,*) 'rho = 0 in q'
qstar(1) = q(1,i,j) - dt2*(ndim-1)/y * q(3,i,j)
qstar(2) = q(2,i,j) - dt2*(ndim-1)/y *
& (rho*u*v)
qstar(3) = q(3,i,j) - dt2*(ndim-1)/y *
& (rho*v**2)
qstar(4) = q(4,i,j) - dt2*(ndim-1)/y *
& v*(q(4,i,j) + press)
c
c # second stage
c
rho = qstar(1)
u = qstar(2)/qstar(1)
v = qstar(3)/qstar(1)
press = gamma1*(qstar(4) - 0.5d0*rho*(u**2 + v**2))
if (rho.eq.0.d0) write(6,*) 'rho = 0 in qstar'
q(1,i,j) = q(1,i,j) - dt*(ndim-1)/y * qstar(3)
q(2,i,j) = q(2,i,j) - dt*(ndim-1)/y *
& (rho*u*v)
q(3,i,j) = q(3,i,j) - dt*(ndim-1)/y *
& (rho*v**2)
q(4,i,j) = q(4,i,j) - dt*(ndim-1)/y *
& v*(qstar(4) + press)
10 continue
return
end