c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@amroc.net
c
c =========================================================
subroutine src(maxmx,maxmy,maxmz,meqn,mbc,ibx,iby,ibz,
& mx,my,mz,q,aux,maux,t,dt,ibnd)
c =========================================================
implicit double precision(a-h,o-z)
include "cuser.i"
c
dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc,1-ibz*mbc:maxmz+ibz*mbc)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc,
& 1-iby*mbc:maxmy+iby*mbc,1-ibz*mbc:maxmz+ibz*mbc)
c
dimension y(1),c(24),w(1,9)
external fcn3
c
n=1
nw=1
tol=1.d-12
do 10 k=1-ibnd*ibz*mbc,mz+ibnd*ibz*mbc
do 10 j=1-ibnd*iby*mbc,my+ibnd*iby*mbc
do 10 i=1-ibnd*ibx*mbc,mx+ibnd*ibx*mbc
rho=q(1,i,j,k)+q(2,i,j,k)
p=gamma1*(q(6,i,j,k)-q(2,i,j,k)*q0-
& 0.5d0*(q(3,i,j,k)**2+q(4,i,j,k)**2+
& q(5,i,j,k)**2)/rho)
y(1)=q(2,i,j,k)/rho
if (p.gt.Pact.and.y(1).gt.0.d0) then
ind=1
z=t
zend=t+dt
call dverk(n,fcn3,z,y,zend,tol,ind,c,nw,w)
if (y(1).lt.0.d0) y(1)=0.d0
if (y(1).gt.1.d0) y(1)=1.d0
q(2,i,j,k)=rho*y(1)
q(1,i,j,k)=rho-q(2,i,j,k)
endif
10 continue
c
return
end
c
c =========================================================
subroutine fcn3(n,x,y,yprime)
c =========================================================
implicit double precision (a-h,o-z)
include "cuser.i"
c
dimension y(n),yprime(n)
c
yprime(1) = 0.d0
if (y(1).le.1.d0.and.y(1).gt.0.d0)
& yprime(1)=-2.d0/treac*dsqrt(y(1))
c
return
end