c
c =====================================================
subroutine ic(maxmx,maxmy,meqn,mbc,mx,my,x,y,dx,dy,q)
c =====================================================
c
c Copyright (C) 2002 Ralf Deiterding
c Brandenburgische Universitaet Cottbus
c
implicit double precision (a-h,o-z)
c
include "cuser.i"
include "ck.i"
c
dimension q(meqn,1-mbc:maxmx+mbc,1-mbc:maxmy+mbc)
dimension x(1-mbc:maxmx+mbc),y(1-mbc:maxmy+mbc)
dimension Xk(leNsp), Yk(leNsp)
c
c SI-units
c
rhol = 0.072d0
ul = 0.d0
pl = 7123.d0
c
rhor = 0.18075d0
ur = -487.34d0
pr = 35594.d0
c
c Convert kg/m**3 into g/cm**3
rhol = rhol * 1.d-3
rhor = rhor * 1.d-3
c Convert Pa = J/m**2 into dynes/cm**2
pl = pl * 1.d1
pr = pr * 1.d1
c Convert m/sec into cm/sec
ul = ul * 1.d2
ur = ur * 1.d2
c
do k=1,Nsp
Xk(k) = 0.d0
enddo
Xk(1) = 1.d0
Xk(6) = 2.d0
Xk(9) = 7.d0
c
sum = 0.d0
do k = 1, Nsp
sum = sum + Xk(k)
enddo
do k=1,Nsp
Xk(k) = Xk(k)/sum
enddo
c
call ckxty (Xk,iwork,rwork,Yk)
call ckmmwy (Yk,iwork,rwork,W)
c
Tl = (pl*W)/(rhol*RU)
Tr = (pr*W)/(rhor*RU)
c
write(6,*) 'left state'
write(6,*) rhol, ul, pl, Tl
write(6,*) 'right state'
write(6,*) rhor, ur, pr, Tr
c
do 150 i=1,mx
do 150 j=1,my
alpha = datan(y(j)/x(i))
call cellave(x(i)-dx/2.d0,y(j)-dy/2.d0,dx,dy,wl)
rho = (1.d0-wl)*rhor + wl*rhol
u = (1.d0-wl)*ur + wl*ul
p = (1.d0-wl)*pr + wl*pl
T = (1.d0-wl)*Tr + wl*Tl
do k=1, Nsp
q(k,i,j) = rho*Yk(k)
enddo
q(Nsp+1,i,j) = rho*u*dcos(alpha)
q(Nsp+2,i,j) = rho*u*dsin(alpha)
q(Nsp+3,i,j) = rho*(0.5d0*u**2+
& avgtabip(T,Yk,hms,Nsp)) - p
q(Nsp+4,i,j) = T
150 continue
return
end
c
c