program mac_2d

c   *** This program solves the 2D Euler equations using MacCormack's
c   *** method.  The equations are solved in non-dimensional form with the
c   *** inlet density, temperature, and velocity as well as the duct length
c   *** used as scaling parameters.  Thus the code solves for
c   *** rho*=rho/rho1, u*=u/u1, p*=p/(rho1*u1^2), Et*=Et/(rho1*u1^2)
c   *** as a function of x*=x/L.

      include '2d_mac.h'

      real cons(Lx,Ly,4), R(Lx,Ly,4), R_old(Lx,Ly,4), prim(Lx,Ly,4),
     &     quasi(Lx,4), work(Lx,Ly,8), sum1(5), error(4)

c   *** Open files and write headers.

      open(unit=1,file='2d_mac.inp')
      open(unit=3,file='grid.dat')
      open(unit=4,file='avg_solution.dat')
      open(unit=8,file='vel.out',form='unformatted')
      open(unit=12,file='rms_error.dat')
      open(unit=13,file='final_error.dat')
      open(unit=14,file='min_max.dat')
      open(unit=15,file='rhs_rms.dat')
      open(unit=19,file='quasi.dat')

       write( 4, 4)
4     format(
     & '#    x/L       rho/rho1      u/u1        v/u1        p/p1')
       write(12,12)
12     format(
     & '#   time      error rho    error u     error v     error p')
       write(13,13)
13     format(
     & '#    x/L      error rho    error u     error v     error p')
       write(14,14)
14     format('#  time   rho_max    u_max    v_max    p_max  rho_min',
     &      '    u_min    v_min    p_min')
       write(15,15)
15     format('#   time        R1_rms      R2_rms      R3_rms',
     &     '      R4_rms       M_min')
19    format(
     & '#    x/L       rho/rho1      u/u1        v/u1        p/p1')

c   *** Read inputs from the file 2d_mac.inp.

      call input( i_restart, nt_max, n_skip, i_plot3d,
     &            i_cfl, dt0, cfl0 )

c   *** Initialize constants etc.

      call init( )

c   *** Read the grid file.

      call read_grid( )
 
c   *** Compute the grid metrics.

      call metrics( )

c   *** Either read or generate a new velocity field.

      if( i_restart .eq. 1 ) then
         open(unit=2,file='vel.in',form='unformatted')
         call read_field( cons, prim )
      else
         call initial_field( cons, prim, quasi )
      end if

c   *** Find the quasi-1-D exact solution.

      call quasi_exact( M_inlet, p_ratio, grid, cent, Nx, Ny, 
     &     quasi(1,1), quasi(1,2), quasi(1,3), quasi(1,4), x_s )

c   *** Compute the initial cfl number and initial statistics.

      if( i_cfl .ne. 0 ) call get_cfl( cons, prim, cfl_by_dt )

      call stats( prim, R, work(1,1,1), work(5,1,1), work(9,1,1) )
      write(14,21) time, (work(i,1,1),i=1,8)

c   *** Enter the time stepping loop, which consists of an outer loop over
c   *** the number of time steps and an inner loop over the number of
c   *** Runge-Kutta substeps.

      R_old = 0.0

      do nt=1, nt_max

c      *** Compute the time step.

         if( i_cfl .eq. 0 ) then
            dt = dt0
         else
            dt = cfl0 / cfl_by_dt
         end if

c      *** Enter the Runge-Kutta loop.

         do nrk=1, nrk_max

c         *** Compute the current right hand side.  Choose between MacCormack
c         *** and first order upwind.

            if( i_scheme .eq. 0 ) then
              call rhs_mac(cons, prim, R, nrk, work(1,1,1), work(1,1,5))
            else if( i_scheme .eq. 1 ) then
              call rhs_up1(cons, prim, R, nrk, work(1,1,1), work(1,1,5))
            else 
              print *, 'Error: no scheme for i_scheme = ', i_scheme
            end if

c         *** Perform the solution update.

            do n=1, 4
               do j=2, Ny+1
                  do i=2, Nx+1
                     cons(i,j,n) = ?
                     R_old(i,j,n) = R(i,j,n)
                  end do
               end do
            end do

c      *** Apply the boundary conditions.

            call bc_in_ex( cons, prim, beta(nrk)*dt )
            call bc_uw_lw( cons, prim )

c      *** Compute the primitive variables from the conservative solution.

            call get_prim( cons, prim )

         end do
 
         time = time + dt

c      *** Compute the cfl number.

         if( i_cfl .ne. 0 ) call get_cfl( cons, prim, cfl_by_dt )

c      *** Compute and write averaged variables and deviation from 
c      *** the quasi-1-D solution.

         if( mod(nt,n_skip) .eq. 0 .or. nt .eq. 1 ) then

            call stats( prim, R, work(1,1,1), work(5,1,1), work(9,1,1) )
            write(14,21) time, (work(i,1,1),i=1,8)
            write(15,20) time, (work(i,1,1),i=9,13)
20          format(1p,6e12.4)
21          format(f8.4,8f9.4)

            error = 0.0
            do i=2, Nx+1
               sum1 = 0.0
               do j=2, Ny+1
                  do n=1, 4
                     sum1(n) = sum1(n) + prim(i,j,n)*vol(i,j)
                  end do
                  sum1(5) = sum1(5) + vol(i,j)
               end do
               do n=1, 4
                  sum1(n) = sum1(n) / sum1(5)
                  error(n) = error(n) + ( sum1(n) - quasi(i,n) )**2
               end do
            end do
            do n=1, 4
               error(n) = sqrt( error(n)*Nx_inv )
            end do
            write(12,20) time, (error(n), n=1, 4)

c      *** Check to make sure the solution is not blowing up.

            if( work(8,1,1) .lt. 0.01 ) then
            print *, 'Solution blowing up, terminating on time step ',nt
               goto 99
            end if

         end if

      end do

c   *** Write plot3d file if desired.

99    if( i_plot3d .eq. 1 ) then

         open(unit=10,file='sol.xyz')
         open(unit=11,file='sol.q')

         write(10,*) Nx+2, Ny+2
         do n=1, 2
            do j=1, Ny+2
               do i=1, Nx+2
                  write(10,*) grid(i,j,n)
               end do
            end do
         end do

         write(11,*) Nx+2, Ny+2
         write(11,*) 0.0, 0.0, 0.0, 0.0
         do n=1, 4
            do j=1, ny+2
               do i=1, Nx+2
                  write(11,*) prim(i,j,n)
               end do
            end do
         end do

      end if

c   *** Compute and write the solution averaged in y.

      do n=1, 4
         error(n) = 0.0
      end do

      do i=1, Nx+2
         x = cent(i,1,1)
         sum1 = 0.0
         do j=2, Ny+1
            do n=1, 4
               sum1(n) = sum1(n) + prim(i,j,n)*vol(i,j)
            end do
            sum1(5) = sum1(5) + vol(i,j)
         end do
         do n=1, 4
            sum1(n) = sum1(n) / sum1(5)
            error(n) = sum1(n) - quasi(i,n)
         end do
         sum1(4) = sum1(4)*p_fac
         write(13,20) x, (error(n), n=1, 4)
         write( 4,20) x, (  sum1(n), n=1, 4)
         write(19,20) x, (quasi(i,n), n=1, 3), quasi(i,4)*p_fac
      end do

c   *** Write the final solution for use in a restart.

      call write_field( cons, 8 )

      stop
      end


      function Area( x )

      pi = acos( -1.0 )

      Area = ?

      return
      end

      
      function dAreadx( x )

      pi = acos( -1.0 )

      dAreadx =  ?

      return
      end


      subroutine bc_in_ex( cons, prim, dt_bc )

      include '2d_mac.h'

      real cons(Lx,Ly,4), prim(Lx,Ly,4)

      i=1
      do j=1, Ny+2
         prim(i,j,1) = ?
         prim(i,j,2) = ?
         prim(i,j,3) = ?
         prim(i,j,4) = ?
         cons(i,j,1) = ?
         cons(i,j,2) = ?
         cons(i,j,3) = ?
         cons(i,j,4) = ?
      end do

      i=Nx+2
      do j=1, Ny+2
         dx_inv = 1.0/( cent(i,j,1) - cent(i-1,j,1) )
         c_exit = ?
         drhodx = ?
         dudx   = ?
         dvdx   = ?
         dpdx   = ?
         prim(i,j,1) = ?
         prim(i,j,2) = ?
         prim(i,j,3) = ?
         prim(i,j,4) = ?
         cons(i,j,1) = ?
         cons(i,j,2) = ?
         cons(i,j,3) = ?
         cons(i,j,4) = ?

      end do
            
      return
      end


      subroutine bc_uw_lw( cons, prim )

      include '2d_mac.h'

      real cons(Lx,Ly,4), prim(Lx,Ly,4)

c   *** Lower boundary.  Assumes streamline

      j=1
      do i=1, Nx+2
         U_wall = ( prim(i,j+1,2)*sy(i,j+1,2) - 
     &              prim(i,j+1,3)*sy(i,j+1,1)   )*ds_x_inv(i,j+1)
         ub =  2*U_wall*sy(i,j+1,2)*ds_x_inv(i,j+1) - prim(i,j+1,2)      
         vb = -2*U_wall*sy(i,j+1,1)*ds_x_inv(i,j+1) - prim(i,j+1,3)      
         pb = 2*prim(i,j+1,4) - prim(i,j+2,4) 
         rhob = prim(i,j+1,1)*(pb/prim(i,j+1,4))
         cons(i,j,1) = rhob
         cons(i,j,2) = rhob*ub
         cons(i,j,3) = rhob*vb
         cons(i,j,4) = c0_inv*pb + 0.5*rhob*( ub**2 + vb**2 )
         prim(i,j,1) = rhob
         prim(i,j,2) = ub
         prim(i,j,3) = vb
         prim(i,j,4) = pb
      end do

c   *** Upper boundary.  Assumes streamline

      j=Ny+2
      do i=1, Nx+2
         U_wall = ( prim(i,j-1,2)*sy(i,j,2) - 
     &              prim(i,j-1,3)*sy(i,j,1)   )*ds_x_inv(i,j)
         ub =  2*U_wall*sy(i,j,2)*ds_x_inv(i,j) - prim(i,j-1,2)      
         vb = -2*U_wall*sy(i,j,1)*ds_x_inv(i,j) - prim(i,j-1,3)      
         pb = 2*prim(i,j-1,4) - prim(i,j-2,4) 
         rhob = prim(i,j-1,1)*(pb/prim(i,j-1,4))
         cons(i,j,1) = rhob
         cons(i,j,2) = rhob*ub
         cons(i,j,3) = rhob*vb
         cons(i,j,4) = c0_inv*pb + 0.5*rhob*( ub**2 + vb**2 )
         prim(i,j,1) = rhob
         prim(i,j,2) = ub
         prim(i,j,3) = vb
         prim(i,j,4) = pb
      end do

      return
      end


      subroutine get_cfl( cons, prim, cfl_by_dt )

      include '2d_mac.h'

      real cons(Lx,Ly,4), prim(Lx,Ly,4)

      cflx_by_dt = 0.0
      cfly_by_dt = 0.0

c   *** Base the CFL on the velocity components normal to the cell faces

      do j=2, Ny+1
         do i=2, Nx+1
           u_x = ( 0.5*( prim(i,j,2) + prim(i-1,j,2) )* sx(i,j,1) +
     &             0.5*( prim(i,j,3) + prim(i-1,j,3) )* sx(i,j,2)
     &           )*ds_y_inv(i,j)
           u_y = ( 0.5*( prim(i,j,2) + prim(i,j-1,2) )* sy(i,j,1) +
     &             0.5*( prim(i,j,3) + prim(i,j-1,3) )* sy(i,j,2)
     &           )*ds_x_inv(i,j)
           c = sqrt( gamma_gas*prim(i,j,4)/prim(i,j,1) )
           cflx_by_dt = max( cflx_by_dt, (abs(u_x)+c)*ds_x_inv(i,j) )
           cfly_by_dt = max( cfly_by_dt, (abs(u_y)+c)*ds_y_inv(i,j) )
         end do
      end do

      cfl_by_dt = cflx_by_dt + cfly_by_dt

      return
      end


      subroutine get_prim( cons, prim )

      include '2d_mac.h'

      real cons(Lx,Ly,4), prim(Lx,Ly,4)

      do j=1, Ny+2
         do i=1, Nx+2
            prim(i,j,1) = ?
            prim(i,j,2) = ?
            prim(i,j,3) = ?
            prim(i,j,4) = ?
         end do
      end do
            
      return
      end


      subroutine init( )

      include '2d_mac.h'

      pi = acos(-1.0)

      c0 = gamma_gas - 1.0
      c0_inv = 1.0/c0

      p_fac = gamma_gas*M_inlet**2
      p_fac_inv = 1.0/p_fac

      p_exit = p_ratio*p_fac_inv

      Nx_inv = 1.0/Nx
      NxNy_inv = 1.0/(Nx*Ny)

c   *** Define weights for Runge-Kutta time stepping.

      if( nrk_max .eq. 1 ) then
         gamma(1) = 1.0
         zeta( 1)  = 0.0
      else if( nrk_max .eq. 2 ) then
         gamma(1) = 1.0
         gamma(2) = 0.5
         zeta(1)  = 0.0
         zeta(2)  =-0.5
      else if( nrk_max .eq. 3 ) then
         gamma(1) = 8.0/15.0
         gamma(2) = 5.0/12.0
         gamma(3) = 3.0/4.0
         zeta(1)  = 0.0
         zeta(2)  =-17.0/60.0
         zeta(3)  =-5.0/12.0
      else
         print *, 'ERROR: chosen Runge-Kutta scheme not yet implemented'
         print *, 'you specified nrk = ', nrk
         stop
      end if

      do n=1, nrk_max
         beta(n) = gamma(n) + zeta(n)
      end do

      time = 0.0

      return
      end


      subroutine initial_field( cons, prim, quasi )

      include '2d_mac.h'

      real cons(Lx,Ly,4), prim(Lx,Ly,4), quasi(Lx,4)

c   *** Sovle the quasi-1-D problem without a shock.

      call quasi_exact( M_inlet, 0.0, grid, cent, Nx, Ny, 
     &     quasi(1,1), quasi(1,2), quasi(1,3), quasi(1,4), x_s )

c   *** Initialize the conservative variables on the interior.  Estimate
c   *** the vertical velocity from the slope of the upper boundary.

      do i=2, Nx+1
         slope = ( grid(i+1,Ny+2,2)-grid(i,Ny+2,2) )/
     &           ( grid(i+1,Ny+2,1)-grid(i,Ny+2,1) )
         h_inv = 1.0/grid(i,Ny+2,2)
         do j=2, Ny+1
            prim(i,j,1) = ?
            prim(i,j,2) = ?
            prim(i,j,3) = quasi(i,2)*slope*cent(i,j,2)*h_inv
            prim(i,j,4) = ?
            cons(i,j,1) = ?
            cons(i,j,2) = ?
            cons(i,j,3) = ?
            cons(i,j,4) = ?
         end do
      end do

c   *** Initialize the variables on the inlet and exit boundaries

      i=1
      do j=1, Ny+2
         cons(i,j,1) = cons(i+1,j,1)
         cons(i,j,2) = cons(i+1,j,2)
         cons(i,j,3) = 0.0
         cons(i,j,4) = cons(i+1,j,4)
         prim(i,j,1) = prim(i+1,j,1)
         prim(i,j,2) = prim(i+1,j,2)
         prim(i,j,3) = 0.0
         prim(i,j,4) = prim(i+1,j,4)
      end do

      i=Nx+2
      do j=1, Ny+2
         cons(i,j,1) = cons(i-1,j,1)
         cons(i,j,2) = cons(i-1,j,2)
         cons(i,j,3) = 0.0
         cons(i,j,4) = cons(i-1,j,4)
         prim(i,j,1) = prim(i-1,j,1)
         prim(i,j,2) = prim(i-1,j,2)
         prim(i,j,3) = 0.0
         prim(i,j,4) = prim(i-1,j,4)
      end do

c   *** Initialize the variables on the upper and lower boundaries

      call bc_uw_lw( cons, prim )
 
      return
      end


      subroutine input( i_restart, nt_max, n_skip, i_plot3d,
     &                  i_cfl, dt0, cfl0 )

      include '2d_mac.h'

      read(1,10) i_restart, nt_max, n_skip, i_plot3d,
     &           i_cfl, dt0, cfl0, i_scheme, nrk_max, eps_diss, 
     &           M_inlet, p_ratio
10    format(i6,/,i6,/,i6,/,i6,//,
     &       i6,/,e15.4,/,f12.5,/,i6,/,i6,/,f12.5,//,
     &       f12.5,/,f12.5)

      if( i_scheme .eq. 0 ) nrk_max = 2

      print *, 'i_restart, nt_max, n_skip, i_plot3d = '
      print *,  i_restart, nt_max, n_skip, i_plot3d
      print *, 'i_cfl, dt0, cfl0, i_scheme, nrk_max, eps_diss = '
      print *,  i_cfl, dt0, cfl0, i_scheme, nrk_max, eps_diss
      print *, 'M_inlet, p_ratio = '
      print *,  M_inlet, p_ratio

      return
      end


      subroutine metrics( )

      include '2d_mac.h'

c   *** Compute cell center coordinates.

      do n=1, 2
         do j=1, Ny+2
            do i=1, Nx+2
               cent(i,j,n) = 0.25*( grid(i,j  ,n) + grid(i+1,j  ,n) +
     &                              grid(i,j+1,n) + grid(i+1,j+1,n)   )
            end do
         end do
      end do

c   *** Compute the metrics for an east face.

      do i=1, Nx+3
         do j=1, Ny+2
            dxe = grid(i,j+1,1) - grid(i,j,1)
            dye = grid(i,j+1,2) - grid(i,j,2)
            sx(i,j,1) = ?
            sx(i,j,2) = ?
            ds_y_inv(i,j) = 1.0 / sqrt( dxe**2 + dye**2 )
         end do
      end do
            
c   *** Compute the metrics for a north face.

      do i=1, Nx+2
         do j=1, Ny+3
            dxn = grid(i+1,j,1) - grid(i,j,1)
            dyn = grid(i+1,j,2) - grid(i,j,2)
            sy(i,j,1) = ?
            sy(i,j,2) = ?
            ds_x_inv(i,j) = 1.0 / sqrt( dxn**2 + dyn**2 )
         end do
      end do
            
c   *** Compute the cell volume.

      do j=1, Ny+2
         do i=1, Nx+2
            a = sqrt( sx(i  ,j  ,1)**2 + sx(i  ,j  ,2)**2 )
            b = sqrt( sy(i  ,j  ,1)**2 + sy(i  ,j  ,2)**2 )
            c = sqrt( sx(i+1,j  ,1)**2 + sx(i+1,j  ,2)**2 )
            d = sqrt( sy(i  ,j+1,1)**2 + sy(i  ,j+1,2)**2 )
            p = sqrt( (grid(i+1,j+1,1)-grid(i  ,j,1))**2 +
     &                (grid(i+1,j+1,2)-grid(i  ,j,2))**2   )
            q = sqrt( (grid(i  ,j+1,1)-grid(i+1,j,1))**2 +
     &                (grid(i  ,j+1,2)-grid(i+1,j,2))**2   )
            arg = (2*p*q)**2 - ( b**2 + d**2 - a**2 -c**2 )**2
            vol(i,j) = 
     &      0.25*sqrt( (2*p*q)**2 - ( b**2 + d**2 - a**2 -c**2 )**2 )
            vol_inv(i,j) = 1.0 / vol(i,j)
         end do
      end do
            
      return
      end


      subroutine quasi_exact( M_inlet, p_ratio, grid, cent, Nx, Ny, 
     &                        rho, u, v, p, x_s )

c   *** Solve the quasi-1D equations exactly.
c   *** Arguments:
c   ***   M_inlet - inlet Mach number (input)
c   ***   p_ratio - exit to inlet pressure ratio (input)
c   ***   grid(Nx+3,Ny+3) - grid
c   ***   cent(Nx+3,Ny+3) - cell center locations
c   ***   Nx      - number of points in x (input)
c   ***   Ny      - number of points in x (input)
c   *** rho(Nx+2) - density scaled by inlet density (output)
c   ***   u(Nx+2) - velocity scaled by inlet velocity (output)
c   ***   v(Nx+2) - velocity scaled by inlet velocity (output)
c   ***   p(Nx+2) - pressure scaled by inlet density times inlet velocity^2 
c   ***             (output)
c   ***   x_s     - shock location scaled by duct length (output)

      parameter( gamma_gas=1.4 )

      real M_inlet, M_save, M, M1, M1_sq, M2_sq
      real grid(Nx+3,Ny+3,2), cent(Nx+3,Ny+3,2)
      real rho(Nx+2), u(Nx+2), v(Nx+2), p(Nx+2)

c   *** Define functions.  f(M) = (A/A*)^2, df(m) = df/dM.

      f( M) = ( ( c2*( 1.0 + c1*M**2 ) )**c3 )/M**2
      df(M) = 2.0/M*( ( c2*( 1.0 + c1*M**2 ) )**c4 - f(M) )

c   *** Define frequently used constants

      c1 = 0.5*( gamma_gas - 1.0 )
      c2 = 2.0/( gamma_gas + 1.0 )
      c3 = ( gamma_gas + 1.0 ) / ( gamma_gas - 1.0 )
      c4 = c3 - 1.0
      c5 = gamma_gas / ( gamma_gas - 1.0 )
      c6 = 1.0 / ( gamma_gas - 1.0 )
      c7 = 1.0/(gamma_gas*M_inlet**2)

      p_exit = c7*p_ratio

c   *** Enter the iteration loop for the shock position.  
c   *** The shock position, x_s, is set to an initial guess.

      x_0 = grid(2   ,1,1)
      x_L = grid(Nx+2,1,1)
      x_s = 0.7654321*( x_L - x_0 ) + x_0
      if( p_ratio .eq. 0.0 ) x_s = 2.0*x_L

      do k=1, 10

c      *** Set reference quantities to values at the inlet.

         x_ref = cent(2,1,1)
         A_ratio_sq_ref = f(M_inlet)
         T_factor = 1.0 + c1*M_inlet**2
         p_factor = T_factor**c5
         M_save = M_inlet

c      *** Enter the loop to march down the inlet.

         do i=2, NX+1

            x  = cent(i  ,1,1)
            xp = cent(i+1,1,1)

c      *** Compute the local area ratio

            A_ratio_sq = A_ratio_sq_ref*(Area(x)/Area(x_ref))**2

c      *** Use Newton iteration to determin the local Mach number.
c      *** As an initial guess, use the Mach number from the prior x station.

            M = M_save
            do j=1, 10
               error = f(M) - A_ratio_sq
               M = M - error/df(M)
            end do

c      *** Compute the remainder of the flow variables.

            T = T_factor / ( 1.0 + c1*M**2 )
            pp = p_factor / ( 1.0 + c1*M**2 )**c5
            rho(i) = pp / T
            p(i) = c7*pp
            u(i) = M/M_inlet*sqrt( T )
            v(i) = 0.0

            M_save = M

c      *** Compute the jumps across the shock wave.
c      *** Note that X_s need not lie on a mesh point.

            if( x .le. x_s .and. x_s .lt. xp .and. M .gt. 1.0 ) then

               x = x_s

               A_ratio_sq = A_ratio_sq_ref*(Area(x)/Area(x_ref))**2

               M1 = M
               do j=1, 10
                  error = f(M1) - A_ratio_sq
                  M1 = M1 - error/df(M1)
               end do

               M1_sq = M1**2
               M2_sq= ( 1.0 + c1*M1_sq )/( gamma_gas*M1_sq - c1 )
               p_jump = 1.0 + c2*gamma_gas*( M1_sq - 1.0 )
               rho_jump = ( M1_sq/c2 ) / ( 1.0 + c1*M1_sq )
               p0_jump= ( rho_jump**gamma_gas/p_jump )**c6

c         *** Adjust the pressure and area ratio reference values
c         *** to account for the loss in total pressure.

               P_factor= p_factor*p0_jump
               A_ratio_sq_ref = A_ratio_sq_ref*p0_jump**2

               M_save = sqrt( M2_sq )

            end if

         end do

c      *** Use secant bisection to update the shock position.

         if( p_ratio .eq. 0.0 ) p_exit = p(Nx+1)

         error = p(Nx+1) - p_exit

         if( k .eq. 1 ) then
            x_s_old = x_s
            x_s = x_s + 0.1
            if( abs(error) .lt. 1.e-7 ) return
         else
            derror = ( error - error_old )/( x_s - x_s_old )
            x_s_old = x_s
            x_s = x_s - error/derror
            if( abs(error) .lt. 1.e-7 ) goto 90
         end if

         error_old = error

      end do

c   *** If the program gets here the shock position iteration failed.

      print *, 'WARNING shock position iteration did not converge'
      print *, 'final error = ', error

90    continue

      rho(   1) = rho(   2)
      u(     1) = u(     2)
      v(     1) = v(     2)
      p(     1) = p(     2)
      rho(Nx+2) = rho(Nx+1)
      u(  Nx+2) = u(  Nx+1)
      v(  Nx+2) = v(  Nx+1)
      p(  Nx+2) = p(  Nx+1)

      return
      end


      subroutine read_field( cons, prim )

      include '2d_mac.h'

      real cons(Lx,Ly,4), prim(Lx,Ly,4)
      real M_inlet0

      read(2) Lx0, Ly0
      read(2) time0
      read(2) M_inlet0, p_ratio0

      if( Lx0 .ne. Lx .or. Ly0 .ne. Ly ) then
         print *, 'Restart dimensions do not match'
         print *, 'Lx0, Ly0 = ', Lx0, Ly0
         stop
      end if

      if( abs(M_inlet0-M_inlet) .gt. 1.0e-10 .or.
     &    abs(p_ratio0-p_ratio) .gt. 1.0e-10      ) then
         print *,
     &   'WARNING gas parameters from restart and input do not match'
         print *, 'M_inlet, p_ratio ='
         write(6,10) 'restart file', M_inlet0, p_ratio0
         write(6,10) 'input   file', M_inlet , p_ratio
10       format(a12,1p5e12.4)
      end if

      read(2) cons

      call get_prim( cons, prim )

      time = time0

      return
      end


      subroutine read_grid( )

      include '2d_mac.h'

      read(3,*) Nxg, Nyg

      if( Nxg .ne. Nx+1 .or. Nyg .ne. Ny+1 ) then
         print *, 'Error grid dimensions wrong, Nxg, Nx+1, Nyg, Ny+1 = '
         print *,  Nxg, Nx+1, Nyg, Ny+1
         stop
      end if

      do j=2, Ny+2
         do i=2, Nx+2
            read(3,*) grid(i,j,1), grid(i,j,2)
         end do
      end do

      do j=2, Ny+2
         grid(1   ,j,1) = 2*grid(2   ,j,1) - grid(3   ,j,1)
         grid(1   ,j,2) = 2*grid(2   ,j,2) - grid(3   ,j,2)
         grid(Nx+3,j,1) = 2*grid(Nx+2,j,1) - grid(Nx+1,j,1)
         grid(Nx+3,j,2) = 2*grid(Nx+2,j,2) - grid(Nx+1,j,2)
      end do

      do i=1, Nx+3
         grid(i,1   ,1) = 2*grid(i,2   ,1) - grid(i,3   ,1)
         grid(i,1   ,2) = 2*grid(i,2   ,2) - grid(i,3   ,2)
         grid(i,Ny+3,1) = 2*grid(i,Ny+2,1) - grid(i,Ny+1,1)
         grid(i,Ny+3,2) = 2*grid(i,Ny+2,2) - grid(i,Ny+1,2)
      end do

      return
      end


      subroutine rhs_mac( cons, prim, R, nrk, flux_x, flux_y )

      include '2d_mac.h'

      real cons(Lx,Ly,4), prim(Lx,Ly,4), R(Lx,Ly,4), 
     &     flux_x(Lx,Ly,4), flux_y(Lx,Ly,4)
      real f(4), g(4)

      if( nrk .eq. 1 ) then
         i_shift = 0
         j_shift = 0
      else
         i_shift = -1
         j_shift = -1
      end if

c      *** Compute fluxes across west faces for the interior.

      do j=2, Ny+1
         do i=2, Nx+2
            is = i + i_shift
            f(1) = ?
            f(2) = ?
            f(3) = ?
            f(4) = ?
            g(1) = ?
            g(2) = ?
            g(3) = ?
            g(4) = ?
            do n=1, 4
               flux_x(i,j,n) = ?
            end do
         end do
      end do

c      *** Compute fluxes across south faces for the interior.

      do i=2, Nx+1
         do j=2, Ny+2
            js = j + j_shift
            f(1) = ?
            f(2) = ?
            f(3) = ?
            f(4) = ?
            g(1) = ?
            g(2) = ?
            g(3) = ?
            g(4) = ?
            do n=1, 4
               flux_y(i,j,n) = ?
            end do
         end do
      end do

c      *** Compute fluxes for the j=2, j=Ny+2 boundaries.

      do j=2, Ny+2, Ny
         do i=2, Nx+1
            if( j .eq. 2 ) then
               jj  = j
               jj1 = j+1
            else
               jj  = j-1
               jj1 = j-2
            end if
            p_wall = 0.5*( 3*prim(i,jj,4) - prim(i,jj1,4) )
            flux_y(i,j,1) = ?
            flux_y(i,j,2) = ?
            flux_y(i,j,3) = ?
            flux_y(i,j,4) = ?
         end do
      end do

c      *** Sum fluxes and divide by the volume.  Add dissipation also.

      fac = eps_diss/dt
      do n=1, 4
         do j=2, Ny+1
            do i=2, Nx+1
               diss = fac*(cons(i-1,j,n)-2.0*cons(i,j,n)+cons(i+1,j,n) +
     &                     cons(i,j-1,n)-2.0*cons(i,j,n)+cons(i,j+1,n) )
               R(i,j,n) = ? + diss
            end do
         end do
      end do

      return
      end


      subroutine rhs_up1( cons, prim, R, nrk, flux_x, flux_y )

c   *** Compute the right hand side using a first order upwind scheme.
c   *** You do not need to complete this routine if you just want to
c   *** stick with MacCormack's scheme.

      include '2d_mac.h'

      real cons(Lx,Ly,4), prim(Lx,Ly,4), R(Lx,Ly,4), 
     &     flux_x(Lx,Ly,4), flux_y(Lx,Ly,4)
      real f(4), g(4)

c   *** Lots of good stuff goes in here.

      return
      end


      subroutine stats( prim, rhs, prim_max, prim_min, rhs_rms )

      include '2d_mac.h'

      real prim(Lx,Ly,4), rhs(Lx,Ly,4), prim_max(4), prim_min(4),
     &     rhs_rms(5)
      real mach

      do n=1, 4
         prim_max(n) = 0.0
         prim_min(n) = 1.0e+10
         rhs_rms( n) = 0.0
      end do
      rhs_rms(5) = m_inlet

      do j=2, Ny+1
         do i=2, Nx+1
            if( prim(i,j,4) .lt. prim_min(4) ) then
               ip_min = i
               jp_min = j
            end if
            do n=1, 4
               prim_max(n) = max(prim_max(n),prim(i,j,n))
               prim_min(n) = min(prim_min(n),prim(i,j,n))
               rhs_rms(n) = rhs_rms(n) + rhs(i,j,n)**2
            end do
            mach = 
     &      sqrt( prim(i,j,1)*( prim(i,j,2)**2 + prim(i,j,3)**2 ) /
     &            (gamma_gas*prim(i,j,4)) )
            rhs_rms(5) = min( rhs_rms(5), mach )
         end do
      end do

      do n=1, 4
         rhs_rms(n) = sqrt( rhs_rms(n)*NxNy_inv )
      end do

      return
      end


      subroutine write_field( cons, iunit )

      include '2d_mac.h'

      real cons(Lx,Ly,4)

      write(iunit) Lx, Ly
      write(iunit) time
      write(iunit) M_inlet, p_ratio
      write(iunit) cons

      return
      end