c
c
c =====================================================
subroutine ic(maxmx,maxmy,maxmz,meqn,mbc,mx,my,mz,
& x,y,z,deltax,deltay,deltaz,q)
c =====================================================
c
implicit double precision (a-h,o-z)
c
include "cuser.i"
c
dimension q(meqn, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
dimension x(1-mbc:maxmx+mbc),y(1-mbc:maxmy+mbc),
& z(1-mbc:maxmz+mbc)
c
do 60 k = 1, mz
do 60 j = 1, my
do 60 i = 1, mx
rho = 0.d0
p = 0.d0
velo = 0.d0
dz = z(k)-deltaz*3.d0/8.d0 - z0
do 70 kk = 1, 4
dy = y(j)-deltay*3.d0/8.d0 - y0
do 80 jj = 1, 4
dx = x(i)-deltax*3.d0/8.d0 - x0
do 90 ii = 1, 4
dist = dsqrt(dx**2.+dy**2.+dz**2.)
rhog = rhoamb
pg = pamb
velog = 0.d0
distr = dist/R0
if(distr.le.0.3d0)then
dist=0.3d0*R0
endif
call von_neumann(E0,R0,rhoamb,pamb,dist,3.d0,
& gamma,1.d-3,rhog,pg,velog)
rho = rho + rhog
p = p + pg
velo = velo + velog
dx = dx + deltax/4.d0
90 continue
dy = dy + deltay/4.d0
80 continue
dz = dz + deltaz/4.d0
70 continue
rho = rho / 64.d0
p = p / 64.d0
velo = velo / 64.d0
dz = z(k) - z0
dy = y(j) - y0
dx = x(i) - x0
dist = dsqrt(dx**2.+dy**2.+dz**2.)
distr = dist/R0
if(distr.le.0.3d0)then
dist=0.3d0*R0
endif
u = velo*dx/dist
v = velo*dy/dist
w = velo*dz/dist
q(1,i,j,k) = rho
q(2,i,j,k) = rho*u
q(3,i,j,k) = rho*v
q(4,i,j,k) = rho*w
q(5,i,j,k) = p/gamma1 + 0.5d0*rho*(u*u+v*v+w*w)
60 continue
return
end