program sodic

c        Establishes initial plot3d grid, solution files, and a restart 
c        file for nparc2d from the given flow conditions for Sod's 
c        shocktube problem. 
c        
c-----------------------------------------------------------------------

      parameter ( nbmax =    1 )

      parameter ( nimax =  110 )
      parameter ( njmax =   25 )

      dimension ni (nbmax)
      dimension nj (nbmax)

      dimension q  (nimax,njmax,nbmax,4)
      dimension x  (nimax,njmax,nbmax)
      dimension y  (nimax,njmax,nbmax)

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

c...Set some constants

      gamma  = 1.4

      gp1    = gamma + 1.0
      gm1    = gamma - 1.0

c...Generate an evenly-spaced grid (units are nondimensional)

      nblock = 1
      ni(1)  = 100
      nj(1)  = 21

      xa = 0.0
      xb = 1.0
      ya = 0.0
      yb = 0.2

      do  i = 1, ni(1)
        do  j = 1, nj(1)
          x(i,j,1) = xa + ( xb - xa ) * ( i - 1.0 ) / ( ni(1) - 1.0 )
          y(i,j,1) = xa + ( yb - ya ) * ( j - 1.0 ) / ( nj(1) - 1.0 )
        enddo
      enddo

c...Generate left side solution. (units are nondimensional).

      pres1 = 1.0
      rho1  = 1.0
      u1    = 0.0
      v1    = 0.0
      et1   = pres1 / gm1

c...Generate right side solution. (units are nondimensional).

      pres2 = 0.1
      rho2  = 0.125
      u2    = 0.0
      v2    = 0.0
      et2   = pres2 / gm1

c...Load the solution with the diaphragm at x = 0.5.

      do  i = 1, ni(1)
        do  j = 1, nj(1)
          if ( x(i,j,1) .lt. 0.5 ) then
            q(i,j,1,1) = rho1
            q(i,j,1,2) = rho1 * u1
            q(i,j,1,3) = rho1 * v1
            q(i,j,1,4) = et1
          else
            q(i,j,1,1) = rho2
            q(i,j,1,2) = rho2 * u2
            q(i,j,1,3) = rho2 * v2
            q(i,j,1,4) = et2
          endif
        enddo
      enddo

c...Output the grid file

      open ( 21, file='nsxic.fmt', form='formatted' )
      write ( 21, * ) nblock
      write ( 21, * ) ( ni(mk), nj(mk), mk = 1, nblock ) 
      do  mk = 1, nblock 
        write ( 21, * ) (( x(i,j,mk), i=1,ni(mk)), j=1,nj(mk)),
     &                  (( y(i,j,mk), i=1,ni(mk)), j=1,nj(mk))
      enddo

c...Output the solution file 

      t1 = 0.0
      t2 = 0.0
      t3 = 1.0
      t4 = 0.0

      open ( 22, file='nsqic.fmt', form='formatted' )
      write ( 22, * ) nblock
      write ( 22, * ) ( ni(mk), nj(mk), mk = 1, nblock ) 
      do  mk = 1, nblock 
        write ( 22, * ) t1, t2, t3, t4
        write ( 22, * ) (((q(i,j,mk,m),i=1,ni(mk)),j=1,nj(mk)),m=1,4)
      enddo

c...Nondimensionalize q for output to nparc

      rref   = 1.0
      aref   = sqrt( gamma * pres1 / rho1 )
      rlen   = 1.0

      do  mk = 1, nblock
        do  i = 1, ni(mk)
          do  j = 1, nj(mk)
            x(i,j,mk)   = x(i,j,mk) / rlen
            y(i,j,mk)   = y(i,j,mk) / rlen
            q(i,j,mk,1) = q(i,j,mk,1) / ( rref )
            q(i,j,mk,2) = q(i,j,mk,2) / ( rref * aref )
            q(i,j,mk,3) = q(i,j,mk,3) / ( rref * aref )
            q(i,j,mk,4) = q(i,j,mk,4) / ( rref * aref * aref )
          enddo
        enddo
      enddo

c...Output the nparc2d restart file (unformatted, fort.2)

      nc1 = 0

      write ( 2 ) nc1

      do  mk = 1, nblock
	write ( 2 ) ni(mk), nj(mk)
	write ( 2 ) (( x(i,j,mk), i=1,ni(mk)), j=1,nj(mk)),
     &              (( y(i,j,mk), i=1,ni(mk)), j=1,nj(mk))
        write ( 2 ) ((( q(i,j,mk,m), i=1,ni(mk)), j=1,nj(mk)), m=1,4 )
      enddo

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

      stop
      end