c
c
c =====================================================
subroutine ic(maxmx,meqn,mbc,mx,x,dx,q)
c =====================================================
implicit none
c
include 'heairacesf6.i'
include 'cuser.i'
c
integer maxmx, meqn, mbc, mx
double precision dx
double precision q(meqn,1-mbc:maxmx+mbc)
double precision x(1-mbc:maxmx+mbc)
integer i, m
double precision eta, wght, f, p, yair
do i=1-mbc, maxmx+mbc
do m=1,meqn
q(m,i) = 0.0d0
enddo
enddo
c
c set up the Riemann problem at the diaphragm
do i=1,mx
eta = (1.0d0+tanh((x(i)-dloc)/zicki))*0.5d0
yair = fraction_air + (1.0d0-fraction_air)*eta
c set the temperature consistently
q(8,i) = t_hemx*(1.0d0-eta)+t_air*eta
c set density and velocity
wght = 1.0d0/((1.0d0-yair)/eqs_Wo_he+yair/eqs_Wo_air)
p = p_hemx*(1.0d0-eta)+p_air*eta
q(1,i) = p*wght/eqs_Ro/q(8,i)
q(2,i) = rhou_hemx*(1.0d0-eta) + rhou_air*eta
c consistently determine total energy
f = (cp_air*yair+cp_he*(1.0d0-yair))*wght/eqs_Ro-1.0d0
q(5,i) = f*p + 0.5d0*(q(2,i)**2)/q(1,i)
c -no AIR, all sf6
q(6,i) = q(1,i)*yair
q(7,i) = 0.0d0
c - 1d simulation, v,w=0
q(3,i) = 0.d0
q(4,i) = 0.d0
end do
IF(trunc.eq.1) THEN
c # overwrite the diaphragm getting rid of all the He
c # and then paint the shock in at x=sloc
do i = 1,mx
q(1,i) = rho_air
q(2,i) = rhou_air
q(3,i) = 0.d0
q(5,i) = e_air
c # all AIR
q(6,i) = 1.0*rho_air
q(7,i) = 0.d0
q(8,i) = t_air
if (x(i).lt.sloc) then
q(1,i) = rho_Sair
q(2,i) = rhou_Sair
q(3,i) = 0.d0
q(5,i) = e_Sair
c # all AIR
q(6,i) = 1.0*rho_Sair
q(7,i) = 0.d0
q(8,i) = t_Sair
end if
end do
endif
c
c make the sf6 curtain
do i=1,mx
if((x(i).gt.crtnL-10.0d0*zicki)
$ .and.(x(i).lt.crtnR+10.0d0*zicki)) then
eta = (1.0d0+tanh((x(i)-crtnL)/zicki))
$ *(1.0d0-tanh((x(i)-crtnR)/zicki))*0.25d0
c set the temperature consistently
q(8,i) = t_air*(1.0d0-eta)+t_sf6*eta
c set density and velocity
wght = 1.0d0/((1.0d0-eta)/eqs_Wo_air+eta/eqs_Wo_sf6)
q(1,i) = p_air*wght/eqs_Ro/q(8,i)
q(2,i) = rhou_sf6*eta + rhou_air*(1.0d0-eta)
c consistently determine total energy
f = (cp_air*(1.0d0-eta)+cp_sf6*eta)*wght/eqs_Ro-1.0d0
q(5,i) = f*p_air + 0.5d0*(q(2,i)**2)/q(1,i)
c -no AIR, all sf6
q(6,i) = q(1,i)*(1.0d0-eta)
q(7,i) = q(1,i)*eta
end if
end do
return
end