program stubeic

c        Computes a uniformly-spaced, planar grid and an initial  
c        solution for the shock-tube problem.  Assumptions:
c        1) Planar domain, 2) Perfect gas.
c
c-----------------------------------------------------------------------

      Implicit none

      Integer nim
      Integer njm

      Parameter ( nim = 200 )
      Parameter ( njm =  50 )

      Integer iprob
      Integer i, j
      Integer m
      Integer ni, nj

      Real q (nim,njm,4)  
      Real x (nim,njm)
      Real y (nim,njm)

      Real a1, p1, r1, u1
      Real a4, p4, r4, u4
      Real gam
      Real xl, xr, xd, yh
      Real rmach, alpha, reyn, time
      Real gascon
      Real aref, eref, pref, rref, tref  

c-----------------------------------------------------------------------

c.....Set gas variables.

      gam    = 1.4
      gascon = 1716.0

c.....Set reference conditions (pressure = 10 psi, temperature = 520 K).

      pref = 10.0 * 144.0
      tref = 520.0
      rref = pref / ( tref * gascon )
      aref = sqrt( gam * gascon * tref )
      eref = rref * aref * aref

c.....Set problem to 10/1 pressure ratio.

      p4 = 1.0 * pref
      r4 = 1.0 * rref
      u4 = 0.0

      p1 = 0.1 * pref
      r1 = 0.125 * rref
      u1 = 0.0

      xl = 0.0
      xr = 1.0
      xd = 0.5
      yh = 0.1

c.....Compute acoustic velocities.

      a1 = sqrt( gam * p1 / r1 )
      a4 = sqrt( gam * p4 / r4 )

c.....Generate an evenly spaced, planar grid.

      write(*,*) ' '
      write(*,*) 'Enter ni and nj'
      read (*,*) ni, nj

      do  i = 1, ni
        do  j = 1, nj
          x(i,j) = xl + ( xr - xl ) * ( i - 1.0 ) / ( ni - 1.0 )
          y(i,j) = yh * ( j - 1.0 ) / ( nj - 1.0 )
        enddo
      enddo

      open ( unit=7, file='stube.x.fmt' )
      write ( 7, * ) ni, nj
      write ( 7, * ) (( x(i,j), i=1,ni), j=1,nj ),
     &               (( y(i,j), i=1,ni), j=1,nj )

c.....Generate initial solution.

      do  i = 1, ni
        do  j = 1, nj
          q(i,j,2) = 0.0
          q(i,j,3) = 0.0
          if ( x(i,j) .lt. xd ) then
            q(i,j,1) = r4 / rref
            q(i,j,4) = ( p4 / ( gam - 1.0 ) ) / eref
          else
            q(i,j,1) = r1 / rref
            q(i,j,4) = ( p1 / ( gam - 1.0 ) ) / eref
          endif
        enddo
      enddo 

c.....Output the solution into file st0q.fmt (2d plot3d format)

      open ( unit=8, file='stubeic.q.fmt' )
      rmach = 0.0
      alpha = 0.0
      reyn  = 1.0
      time  = 0.0 
      write ( 8, * ) ni, nj
      write ( 8, * ) rmach, alpha, reyn, time
      write ( 8, * ) ((( q(i,j,m), i=1,ni), j=1,nj), m=1,4 )

c-----------------------------------------------------------------------

      stop
      end