PROGRAM normic

!   Computes a uniformly-spaced, 3d grid for a 2.0 ft long straight 
!   duct with a diameter of 1.0 ft and initializes a flow with a normal 
!   shock midway through the duct.
!
!-----------------------------------------------------------------------

      Implicit none

      Integer, parameter :: nim = 202
      Integer, parameter :: njm =  11
      Integer, parameter :: nkm =   5

      Integer i, j, k, m
      Integer nc1
      Integer ni, nj, nk
      Integer nblock

      Real q (nim,njm,nkm,5)  
      Real x (nim,njm,nkm)
      Real y (nim,njm,nkm)
      Real z (nim,njm,nkm)

      Real a1, et1, m1, p1, r1, t1, u1, v1, w1
      Real a2, et2, m2, p2, r2, t2, u2, v2, w2
      Real gam, gm1, gp1
      Real gascon
      Real pi
      Real alpha, reyn, time
      Real xl, xr, xd, rad
      Real thet
      Real p2p1, r2r1, u2u1

!-----------------------------------------------------------------------

   pi     = 4.0 * atan(1.0)

!..Set gas variables.

   gam    = 1.4
   gp1    = gam + 1.0
   gm1    = gam - 1.0
   gascon = 1716.0

!..Set conditions on supersonic side of shock.

   m1  = 1.3
   p1  = 10.0 * 144.0
   t1  = 520.0
   r1  = p1 / ( t1 * gascon )
   a1  = sqrt( gam * gascon * t1 )
   u1  = m1 * a1
   v1  = 0.0
   w1  = 0.0
   et1 = p1 / gm1 + 0.5 * r1 * u1 * u1

!..Create conditions on subsonic side of shock.

   r2r1 = gp1 * m1*m1 / ( 2.0 + gm1 * m1*m1 )
   u2u1 = 1.0 / r2r1
   p2p1 = 1.0 + 2.0*gam*(m1*m1-1.0)/gp1
   m2   = sqrt( ( 1.0 + 0.5*gm1*m1*m1 ) / ( gam*m1*m1 - 0.5*gm1 ) )
      
   p2   = p2p1 * p1
   r2   = r2r1 * r1
   u2   = u2u1 * u1
   v2   = 0.0
   w2   = 0.0
   et2  = p2 / gm1 + 0.5 * r2 * u2 * u2 

!..Set duct size and shock location.
   
   xl  = 0.0
   xr  = 2.0
   xd  = 1.0
   rad = 0.5

!..Generate an evenly spaced, 3D grid.

   ni = 202
   nj =  11
   nk =   5

   do  i = 1, ni
     do  j = 1, nj
       do  k = 1, nk
         thet     = 10.0 * ( k - 1.0 ) / ( nk - 1.0 ) * pi / 180.0
         x(i,j,k) = xl + ( xr - xl ) * ( i - 1.0 ) / ( ni - 1.0 )
         y(i,j,k) = rad * ( j - 1.0 ) / ( nj - 1.0 ) * cos(thet)
         z(i,j,k) = rad * ( j - 1.0 ) / ( nj - 1.0 ) * sin(thet)
       enddo
     enddo
   enddo

   open ( unit=7, file='norm.x', form='unformatted', status='unknown' )
   nblock = 1
   write(7) nblock
   write(7) ni, nj, nk
   write(7) ((( x(i,j,k), i=1,ni), j=1,nj), k=1,nk),  &
            ((( y(i,j,k), i=1,ni), j=1,nj), k=1,nk),  &
            ((( z(i,j,k), i=1,ni), j=1,nj), k=1,nk)

!..Load the initial solution using state 1 as reference.

   do  i = 1, ni
     do  j = 1, nj
       do  k = 1, nk
         if ( x(i,j,k) .lt. xd ) then
           q(i,j,k,1) = r1      / r1
           q(i,j,k,2) = r1 * u1 / ( r1 * a1 )
           q(i,j,k,3) = r1 * v1 / ( r1 * a1 )
           q(i,j,k,4) = r1 * w1 / ( r1 * a1 )
           q(i,j,k,5) = et1     / ( r1 * a1 * a1 )
         else
           q(i,j,k,1) = r2      / r1
           q(i,j,k,2) = r2 * u2 / ( r1 * a1 )
           q(i,j,k,3) = r2 * v2 / ( r1 * a1 )
           q(i,j,k,4) = r2 * w2 / ( r1 * a1 )
           q(i,j,k,5) = et2     / ( r1 * a1 * a1 )
         endif
       enddo
     enddo
   enddo

!..Output the solution to file norm.q (3D plot3d unformated)

   open ( unit=8, file='norm.q', form='unformatted', status='unknown' )

   alpha = 0.0
   reyn  = 1.0
   time  = 0.0 

   write (8) nblock
   write (8) ni, nj, nk
   write (8) m1, alpha, reyn, time
   write (8) (((( q(i,j,k,m), i=1,ni), j=1,nj), k=1,nk), m=1,5 )

!..Output the NPARC restart file.

   nc1 = 0

   write (2) nc1
   write (2) ni, nj, nk
   write (2) ((( x(i,j,k), i=1,ni), j=1,nj), k=1,nk),  &
             ((( y(i,j,k), i=1,ni), j=1,nj), k=1,nk),  &
             ((( z(i,j,k), i=1,ni), j=1,nj), k=1,nk)
   write (2) (((( q(i,j,k,m), i=1,ni), j=1,nj), k=1,nk), m=1,5 )

!-----------------------------------------------------------------------

   stop
   end