c =====================================================
subroutine ic(maxmx,maxmy,meqn,mbc,mx,my,x,y,dx,dy,q)
c =====================================================
implicit none
include 'heairacesf6.i'
include 'cuser.i'
integer maxmx,maxmy,meqn,mbc,mx,my
double precision dx,dy
double precision x(1-mbc:maxmx+mbc)
double precision y(1-mbc:maxmy+mbc)
double precision q(meqn, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
double precision xconR,xconL
double precision pi
double precision freq
double precision freqL, freqR
double precision width,eta
double precision wght,f,yair,p
real randm
integer i,j,k,l,m,n
pi = 4.0d0*atan(1.0d0)
c # width of the shock tube
width = GeomMax(2)-GeomMin(2)
c # used to convert modes to freq
freqL = 2.d0*pi*omegaL/width
freqR = 2.d0*pi*omegaR/width
c
do 50 j = 1-mbc,maxmy+mbc
do 50 i = 1-mbc,maxmx+mbc
! --- this insures the any extended bits
! --- in the vector of state will start as zero
do l = 1,meqn
q(l,i,j) = 0.0d0
end do
50 continue
c
c set up the Riemann problem at the diaphragm
do j = 1,my
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,j) = 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,j) = p*wght/eqs_Ro/q(8,i,j)
q(2,i,j) = 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,j) = f*p + 0.5d0*(q(2,i,j)**2)/q(1,i,j)
c -no AIR, all sf6
q(6,i,j) = q(1,i,j)*yair
q(7,i,j) = 0.0d0
c - 1d simulation, v,w=0
q(4,i,j) = 0.d0
end do
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 j = 1,my
do i = 1,mx
q(1,i,j) = rho_air
q(2,i,j) = rhou_air
q(3,i,j) = 0.d0
q(5,i,j) = e_air
c # all AIR
q(6,i,j) = 1.0*rho_air
q(7,i,j) = 0.d0
q(8,i,j) = t_air
if (x(i).lt.sloc) then
q(1,i,j) = rho_Sair
q(2,i,j) = rhou_Sair
q(3,i,j) = 0.d0
q(5,i,j) = e_Sair
c # all AIR
q(6,i,j) = 1.0*rho_Sair
q(7,i,j) = 0.d0
q(8,i,j) = t_Sair
end if
end do
end do
else if ( trunc .eq. 2 ) then
c interpolate solution to current mesh
l = 1
do while ( qin(1,l+1) .lt. x(1) )
l = l + 1
enddo
do i=1, mx
if ( x(i) .lt. sloc ) then
do while ( qin(1,l+1) .lt. x(i) )
l = l + 1
enddo
if ( l .lt. nic ) then
eta = (x(i)-qin(1,l))/(qin(1,l+1)-qin(1,l))
else
eta = 0.0d0
endif
do j=1, my
do m=1,8
q(m,i,j) = qin(m+1,l)+(qin(m+1,l+1)-qin(m+1,l))*eta
enddo
enddo
endif
enddo
end if
c Create the gas SF6 curtain
do j = 1,my
do i = 1,mx
c # compute the location of the curtain with perturbations
xconL = crtnL +ampL*cos(freqL*y(j))
xconR = crtnR +ampR*cos(freqR*y(j))
if((x(i).gt.xconL-10.0d0*zicki)
$ .and.(x(i).lt.xconR+10.0d0*zicki)) then
eta = (1.0d0+tanh((x(i)-xconL)/zicki))
$ *(1.0d0-tanh((x(i)-xconR)/zicki))*0.25d0
c set the temperature consistently
q(8,i,j) = 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,j) = p_air*wght/eqs_Ro/q(8,i,j)
q(2,i,j) = 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,j) = f*p_air + 0.5d0*(q(2,i,j)**2)/q(1,i,j)
c -no AIR, all sf6
q(6,i,j) = q(1,i,j)*(1.0d0-eta)
q(7,i,j) = q(1,i,j)*eta
end if
end do
end do
return
end