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 David Hill
c
IMPLICIT NONE
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)
! ---- shocked (behind shock) first gass
DOUBLE PRECISION u1, P1, rho1, gamma1, c1,T1,rgas1
! ---- unshocked first gass
DOUBLE PRECISION u0, P0, rho0, gamma0, c0,T0,rgas0
! ---- unshocked second gass
DOUBLE PRECISION ub, Pb, rhob, gammab, cb,Tb,rgasb
! ---- molecular weights
DOUBLE PRECISION wmol1, wmolb
! ---- the speed of the shock
DOUBLE PRECISION amach
! ---- used in creating the contact discont
DOUBLE PRECISION rhoavg, rhodiff, rhoval
DOUBLE PRECISION x0,x1,y0,y1,z0, z1, amp, zpert
DOUBLE PRECISION xs,pert_amp, pi, xinterface
DOUBLE PRECISION gas1,amp2
real rand
INTEGER i,j,k
! ---- make pi
pi = 4.0d0*ATAN(1.0d0)
! ---- A Mach 1.5 shock
amach = 10.d0
! ---- one gamma gass
gamma0 = 1.4
gamma1 = 1.4
gammab = 1.4
! ---- different desities
wmol1 = 1.0D0
wmolb = 5.0D0
rgas0 = 8317.D0/wmol1
rgas1 = 8317.D0/wmol1
rgasb = 8317.D0/wmolb
! ----- Unshocked first gass is the reference
u0 = 0.0d0
P0 = 1.0d0
rho0 = 1.0d0
c0 = dsqrt(gamma0*P0/rho0)
T0 = P0/rho0/rgas0
! ---- Unshocked second gass
ub = 0.0d0
Pb = P0
Tb = T0
rhob = Pb/Tb/rgasb
! ---- Behind the shock in the first gass
! ---- From the R-H conditions
gas1=(gamma0-1.D0)/(gamma0+1.D0)
P1=P0*(1.0D0+2.D0*gamma0/(gamma0+1.D0)*(amach*amach-1.D0))
rho1=rho0*(P1/P0+gas1)/(gas1*P1/P0+1.D0)
T1=P1/P0*rho0/rho1*t0
u1=((rho1-rho0)*c0*amach+rho0*u0)/rho1
! ---- the average and jump between the densities
! ---- of the two still gasses
rhoavg = (rho0+rhob)/2.0d0
rhodiff = (rho0-rhob)/2.0d0
amp = 20.0d0
amp2 = .5
pert_amp = 0.1d0
! ---- place the contact discont at the physical location 1/2
x0 = 0.5d0
! ---- place the shock at the physical location 1/4
xs = 0.25
! ---- across the domain, establish the two still gasses
! ---- with an interface, but with out the shock
DO i = 1, mx
DO j = 1, my
! ---- compute the intface location
xinterface = x0 + pert_amp*abs(sin(2.0*pi*y(j)*2.5d0))
!*abs(sin(2.0*pi*z(k)*2.5d0))
! ---- compute the local value of density
rhoval = rhoavg - rhodiff * dtanh(amp*( x(i) - xinterface))
! ---- initialize the conserved variables
q(1,i,j) = rhoval
q(2,i,j) = 0.d0
q(3,i,j) = 0.d0
q(4,i,j) = 0.d0
q(5,i,j) = P0/(gamma0-1.0d0)
q(6,i,j) = rhoval * dtanh(amp* (x(i) - xinterface))
q(7,i,j) = -rhoval * dtanh(amp2* (x(i) - xinterface))
END DO
END DO
! ---- Initialize the shock
DO i = 1, mx
IF ( x(i).le.xs ) THEN
DO j = 1,my,1
q(6,i,j) = q(6,i,j)*rho1/q(1,i,j)
q(7,i,j) = q(7,i,j)*rho1/q(1,i,j)
q(1,i,j) = rho1
q(2,i,j) = rho1 * u1
q(3,i,j) = 0.0d0
q(4,i,j) = 0.0d0
q(5,i,j) = 0.5d0*rho1*u1*u1 + P1/(gamma0-1.0d0)
END DO
END IF
END DO
RETURN
END