c c ===================================================== subroutine ic(maxmx,maxmy,maxmz,meqn,mbc,mx,my,mz, & x,y,z,dx,dy,dz,q) c ===================================================== c c Copyright (C) 2003-2007 California Institute of Technology c Ralf Deiterding, ralf@cacr.caltech.edu 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 r = dsqrt(y(j)*y(j)+z(k)*z(k)) do 60 i = 1, mx if (x(i).le.as) then rho = rhoo u = uo p = po gam = g(2) pin = pinf(2) else rho = rhoi u = ui p = pi gam = g(1) pin = pinf(1) endif c q(1,i,j,k) = rho q(2,i,j,k) = u*rho q(3,i,j,k) = 0.d0 q(4,i,j,k) = 0.d0 q(6,i,j,k) = 1.d0/(gam-1.d0) q(7,i,j,k) = gam*pin/(gam-1.d0) q(5,i,j,k) = p*q(6,i,j,k)+q(7,i,j,k)+ & 0.5d0*q(2,i,j,k)**2/q(1,i,j,k) call tintp(x(i),q(1,i,j,k),meqn) 60 continue c return end c c ************************************************************** c subroutine tintp(x,q,meqn) implicit double precision(a-h,o-z) c include "cuser.i" dimension q(meqn) c do i = 1, Nread-1 if (x.ge.xtab(i).and.x.le.xtab(i+1)) then alpha = (x-xtab(i))/(xtab(i+1)-xtab(i)) rho = (1.d0-alpha)*qtab(1,i) + alpha*qtab(1,i+1) u = (1.d0-alpha)*qtab(2,i) + alpha*qtab(2,i+1) p = (1.d0-alpha)*qtab(3,i) + alpha*qtab(3,i+1) q(1) = rho q(2) = u*rho q(5) = p*q(6) + q(7) + 0.5d0*rho*u**2 endif enddo c return end