c c ===================================================== subroutine ic(maxmx,maxmy,meqn,mbc,mx,my,x,y,dx,dy,q) c ===================================================== c c Copyright (C) 2003-2007 California Institute of Technology c Ralf Deiterding, ralf@amroc.net c implicit double precision (a-h,o-z) c include "cuser.i" c dimension q(meqn,1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc) dimension x(1-mbc:maxmx+mbc),y(1-mbc:maxmy+mbc) c c # circle of high-density gas c idisc = 2 r0 = rc0 x0 = xc0 y0 = yc0 c Yw = 0.d0 Ys = 1.d0 c do 60 i = 1, mx do 60 j = 1, my if (y(j).gt.wlev) then rho = rhoai p = rhoai*grav*y(j) + pai Y1 = Ys else call cellave(x(i)-dx/2.d0,y(j)-dy/2.d0,dx,dy,wl) rho = (1.d0-wl)*rhow + wl*rhos pwh = rhow*grav*y(j) + pai p = (1.d0-wl)*pwh + wl*ps Y1 = (1.d0-wl)*Yw + wl*Ys endif c W = 1.d0/(Y1/Wk(1)+(1.d0-Y1)/Wk(2)) X1 = Y1*W/Wk(1) c q(1,i,j) = rho q(2,i,j) = 0.d0 q(3,i,j) = 0.d0 q(5,i,j) = X1/(g(1)-1.d0)+(1.d0-X1)/(g(2)-1.d0) q(6,i,j) = X1*g(1)*pinf(1)/(g(1)-1.d0)+ & (1.d0-X1)*g(2)*pinf(2)/(g(2)-1.d0) q(4,i,j) = p*q(5,i,j)+q(6,i,j)+ & 0.5d0*(q(2,i,j)**2+q(3,i,j)**2)/q(1,i,j) 60 continue return end