c
c =====================================================
subroutine ic(maxmx,maxmy,maxmz,meqn,mbc,mx,my,mz,
& xc,yc,zc,dx,dy,dz,q)
c =====================================================
c
c Copyright (C) 2003-2007 California Institute of Technology
c Ralf Deiterding, ralf@amroc.net
c
implicit double precision (a-h,l,o-z)
include "cuser.i"
c
dimension q(meqn, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension xc(1-mbc:maxmx+mbc),yc(1-mbc:maxmy+mbc),
& zc(1-mbc:maxmz+mbc)
dimension y(1),c(24),w(1,9)
external fcn
c
n=1
tol=1.d-10
ind=1
nw=1
y(1)=0.d0
lambda=0.d0
if (zc(mz) .lt. sloc) then
z=0.d0
zend=sloc-zc(mz)-0.5d0*dz
tol=1.d-8
ind=1
call dverk(n,fcn,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
lambda = y(1)
endif
c
mcount = 1
dmc = 1.d0/mcount
do k=mz,1,-1
c
do j = 1, my
do i = 1, mx
q(3,i,j,k)=0.d0
q(4,i,j,k)=0.d0
enddo
enddo
c
if ( zc(k) .gt. sloc ) then
do j = 1, my
do i = 1, mx
q(1,i,j,k)=0.d0
q(2,i,j,k)=rho0
q(5,i,j,k)=-D*rho0*U0
if (moving.eq.1) q(5,i,j,k) = 0.d0
q(6,i,j,k)=(1.d0/gamma1+qn)*P0+
& 0.5d0*q(5,i,j,k)**2/q(2,i,j,k)
enddo
enddo
else if ( lambda.lt.1.d0 .and.
& zc(k)-0.5d0*dz.lt.sloc ) then
z=sloc-zc(k)+0.5d0*dz
do j = 1, my
do i = 1, mx
q(1,i,j,k) = 0.d0
q(2,i,j,k) = 0.d0
q(5,i,j,k) = 0.d0
q(6,i,j,k) = 0.d0
enddo
enddo
do m=1,mcount
zend=z+dmc*dz
tol=1.d-8
ind=1
call dverk(n,fcn,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
lambda=y(1)
a=0.d0
if (lambda/clambda.lt.1.d0)
& a=cfact*dsqrt(1.d0-lambda/clambda)
V=Vj*(1.d0-a/gamma)
P=Pj*(1.d0+a)
DU=-D*V
if (moving.eq.1) DU = DU + D
do j = 1, my
do i = 1, mx
q(1,i,j,k)=q(1,i,j,k)+dmc*lambda/V*rho0
q(2,i,j,k)=q(2,i,j,k)+
& dmc*(1.d0-lambda)/V*rho0
q(5,i,j,k)=q(5,i,j,k)+dmc*DU/V*rho0*U0
q(6,i,j,k)=q(6,i,j,k)+
& dmc*(P/gamma1+(1.d0-lambda)/
& V*qn)*P0+0.5d0*q(5,i,j,k)**2/
& (q(1,i,j,k)+q(2,i,j,k))
enddo
enddo
z=zend
enddo
else
a=0.d0
if (1.d0/clambda.lt.1.d0)
& a=cfact*dsqrt(1.d0-1.d0/clambda)
V=Vj*(1.d0-a/gamma)
P=Pj*(1.d0+a)
DU=-D*V
if (moving.eq.1) DU = DU + D
do j = 1, my
do i = 1, mx
q(1,i,j,k) = 1.d0/V*rho0
q(2,i,j,k) = 0.d0
q(5,i,j,k) = DU/V*rho0*U0
q(6,i,j,k) = P/gamma1*P0+
& 0.5d0*q(5,i,j,k)**2/q(1,i,j,k)
enddo
enddo
endif
enddo
c
do k = 1, mz
do j = 1, my
do i = 1, mx
if (dsqrt(xc(i)**2+yc(j)**2).gt.rfi) then
q(1,i,j,k)=0.d0
q(3,i,j,k)=0.d0
q(4,i,j,k)=0.d0
q(2,i,j,k)=rho0
q(5,i,j,k)=-D*rho0*U0
if (moving.eq.1) q(5,i,j,k) = 0.d0
q(6,i,j,k)=(1.d0/gamma1+qn)*P0+
& 0.5d0*q(5,i,j,k)**2/q(2,i,j,k)
endif
enddo
enddo
enddo
return
end
c
c==========================================================
subroutine fcn(n,x,y,yprime)
c==========================================================
implicit double precision (a-h,l,o-z)
include "cuser.i"
c
dimension y(n),yprime(n)
c
a=0.d0
if (y(1)/clambda.lt.1.d0)
& a=cfact*dsqrt(1.d0-y(1)/clambda)
P=Pj*(1.d0+a)
V=Vj*(1.d0-a/gamma)
DU=D*V
yprime(1) = 0.d0
if (y(1).lt.1.d0.and.y(1).ge.0.d0)
& yprime(1)=2.d0/treac*dsqrt(1.d0-y(1))/(DU*U0)
c
return
end