c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
c
c =========================================================
subroutine src(maxmx,meqn,mbc,ibx,mx,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)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc)
c
parameter ( mtab = 10000, mmeqn = 4 )
dimension ql(mmeqn, -1:mtab+2), qr(mmeqn, -1:mtab+2)
dimension Yk1(2), Yk2(2), Ykl(2), Ykr(2)
c
mlim = 1
c
if (ibnd*ibx.gt.0) then
do 100 m=1,meqn
ql(m,-1) = q(m,-1)
qr(m,-1) = q(m,-1)
ql(m,mx+mbc) = q(m,mx+mbc)
qr(m,mx+mbc) = q(m,mx+mbc)
100 continue
endif
c
om = 0.d0
do 110 i=2-ibnd*ibx*mbc,mx-1+ibnd*ibx*mbc
c
c # Reconstruction
c
do m=1,4
call reclim(q(m,i),q(m,i-1),q(m,i+1),
& mlim,om,ql(m,i),qr(m,i))
enddo
c
110 continue
c
call src2(mtab,meqn,mbc,ibx,mx,ql,aux,maux,t,dt,ibnd)
call src2(mtab,meqn,mbc,ibx,mx,qr,aux,maux,t,dt,ibnd)
c
do 300 i=1-ibnd*ibx*mbc, mx+ibnd*ibx*mbc
do 300 m=1,meqn
q(m,i) = 0.5d0*(ql(m,i) + qr(m,i))
300 continue
c
return
end
c
c
c =========================================================
subroutine src2(maxmx,meqn,mbc,ibx,mx,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)
dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc)
c
dimension y(1),c(24),w(1,9)
common /euparam/ e
external fcn3
c
n=1
tol=1.d-8
nw=1
pmax = 0.d0
do 10 i=1-ibnd*ibx*mbc,mx+ibnd*ibx*mbc
rho=q(1,i)+q(2,i)
u=q(3,i)/rho
e=q(4,i)/rho-0.5d0*u**2
temp=gamma1*(e-q(2,i)/rho*q0)
if (temp.gt.1.d0) then
ind=1
z=t
zend=t+dt
y(1)=q(2,i)/rho
c if (y(1).lt.0.d0.or.y(1).gt.1.d0) write (6,*) 'rp'
c if (y(2).lt.0.d0.or.y(2).gt.1.d0) write (6,*) 'rp'
call dverk(n,fcn3,z,y,zend,tol,ind,c,nw,w)
c if (y(1).lt.0.d0.or.y(1).gt.1.d0) write (6,*) 'src'
c if (y(2).lt.0.d0.or.y(2).gt.1.d0) write (6,*) 'src'
q(2,i)=rho*y(1)
q(1,i)=rho-q(2,i)
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)
common /euparam/ e
c
temp=gamma1*(e-y(1)*q0)
yprime(1)=-dhalbe*y(1)*exp(-E0/temp)
return
end