      program euler
       
      include 'common.h'       

      integer  i, nt, nv
      integer  ibeg, iend
      real*8   u(nvar, nx),v(nvar, nx), flux(nvar, nx)  
      real*8   x(nx)
      real*8   t, dt, cmax, cfl, tstop 
      real*8   tfreq, df, dx

c ** generate grid **

      call grid (x, dx)
      ibeg = nghost + 1
      iend = nx - nghost

      call init (v, x)
      call primtocon (v, u, ibeg, iend)
      
      dt    = 1.d-4
      cfl   = 0.8d0
      tstop = 0.2
      t     = 0.d0

c ** begin computation **

      do nt = 1, 9999
       
        cmax = 0.d0
        
        call boundary  (u)
        call contoprim (u, v, 1, nx)
        call riemann   (v, flux, cmax, ibeg-1, iend)

        do i  = ibeg, iend 
        do nv = 1, nvar 
          df = flux(nv,i) - flux(nv,i-1)
          u(nv,i) = u(nv,i) - dt/dx*df
        end do
        end do

        t = t + dt
          
        write (*,'(A,I4,A,F8.6,A,F8.6)') 
     &         'out = ', nt,'  time = ',t,'  dt = ',dt

        dt = min(dt*1.1d0, cfl*dx/cmax)

        if (t .ge. tstop) then 
          call write_to_file (v, x, ibeg, iend)
          stop
        end if
      end do
       
      end

c **************************************************       
      subroutine grid (x, dx)
c
c
c
c
c **************************************************       
       
      include  'common.h'
      integer i
      real*8  x(nx), dx, xbeg, xend       

      xbeg =  0.d0
      xend =  1.d0

      dx = (xend - xbeg)/(nx - 2.d0*nghost)
      
      do i = 1, nx
        x(i) = xbeg + (i - nghost - 0.5d0)*dx
      end do
             
      return
      end
c **************************************************       
      subroutine init (v,x)
c
c
c
c
c **************************************************       
       
      include  'common.h'
      real*8   v(nvar, nx), x(nx)
      integer  i       
       
      do i = 1,nx        
        if (x(i) .le. 0.5) then 
         v(irho,i) = 1.0
         v(ivlx,i) = 0.0
         v(iprs,i) = 1.0
        else
         v(irho,i) = 0.125
         v(ivlx,i) = 0.0
         v(iprs,i) = 0.1
        end if  
      end do
      
      return
      end 
c **************************************************       
      subroutine boundary (v)
c
c
c
c
c **************************************************       
       
      include 'common.h'
      real*8   v(nvar, nx)  
      integer  i, nv

      do i  = 1, nghost
      do nv = 1, nvar
        v(nv, i) = v(nv, nghost + 1)
      end do  
      end do
      
      do i  = nx - nghost, nx
      do nv = 1, nvar
        v(nv, i) = v(nv, nx - nghost - 1)
      end do  
      end do
      
      return
      end       

c **************************************************       
      subroutine write_to_file (v, x, ibeg, iend)
c
c
c
c
c **************************************************       

      include 'common.h'
      real*8   v(nvar, nx), x(nx)  
      integer  i, ibeg, iend
      

      open (unit=20, file="data.out", status='unknown')

      do i = ibeg, iend
        write (unit=20, fmt="(4f10.4)") x(i), v(irho, i), v(ivlx, i),
     &                                  v(iprs, i)
      end do

      close(20)
      return
      end
            
c **************************************************       
      subroutine contoprim (u, v, ibeg, iend)
c
c
c
c
c **************************************************       
       
      include 'common.h'
      
      integer  ibeg, iend, i
      real*8   u(nvar, nx), v(nvar, nx)  
      real*8   rho, vel2
      
      do i = ibeg, iend
        rho           = u(irho,i)

        v(irho,i) = rho
        v(ivlx,i) = u(immx,i)/rho
        vel2      = v(ivlx,i)*v(ivlx,i) 
        v(iprs,i) = (gamma - 1.d0)*(u(ienr,i) - 0.5d0*rho*vel2)
      end do    

      return
      end       

c **************************************************       
      subroutine primtocon (v, u, ibeg, iend)
c
c
c
c
c **************************************************       
             
      include 'common.h'
      
      integer  ibeg, iend, i
      real*8   u(nvar, nx), v(nvar, nx)  
      real*8   rho, vx
      
      do i = ibeg, iend
      
        rho = v(irho,i)
        vx  = v(ivlx,i)

        u(irho,i) = rho
        u(immx,i) = rho*vx
        u(ienr,i) = v(iprs,i)/(gamma - 1.d0) + 0.5d0*rho*vx*vx
        
      end do    

      return
      end       

c **************************************************       
      subroutine riemann (v, flux, cmax, ibeg, iend)
c
c
c
c
c **************************************************       
       
      include 'common.h'
      
      integer  ibeg, iend, i, nv
      real*8   v(nvar, nx)  
      real*8   flux(nvar, nx)
      real*8   fr(nvar), fl(nvar), cmax, cmax_loc       
      real*8   rho_l, vx_l, p_l, kin_l
      real*8   rho_r, vx_r, p_r, kin_r
      real*8   mx_l, e_l, cs_l
      real*8   mx_r, e_r, cs_r

     
      do i  = ibeg, iend 

        rho_l = v(irho, i)
        vx_l  = v(ivlx, i)
        p_l   = v(iprs, i)

        rho_r = v(irho, i + 1)
        vx_r  = v(ivlx, i + 1)
        p_r   = v(iprs, i + 1)
      
        cs_r = sqrt(gamma*p_r/rho_r)
        cs_l = sqrt(gamma*p_l/rho_l)

        cmax_loc = max(abs(vx_r) + cs_r, abs(vx_l) + cs_l)
        cmax     = max(cmax, cmax_loc)

        mx_l = rho_l*vx_l
        mx_r = rho_r*vx_r
        
        kin_l = 0.5d0*rho_l*vx_l*vx_l
        kin_r = 0.5d0*rho_r*vx_r*vx_r

        e_l   = p_l/(gamma - 1.d0) + kin_l
        e_r   = p_r/(gamma - 1.d0) + kin_r
      
        flux(irho,i) = 0.5d0*(rho_l*vx_l + rho_r*vx_r  
     &                          + cmax_loc*(rho_l - rho_r))
        flux(immx,i) = 0.5d0*(mx_l*vx_l + mx_r*vx_r  
     &                          + cmax_loc*(mx_l - mx_r))
        flux(ienr,i) = 0.5d0*((e_l + p_l)*vx_l + (e_r + p_r)*vx_r 
     &                          + cmax_loc*(e_l - e_r))
        
        flux(immx,i) = flux(immx,i) + 0.5d0*(p_r + p_l)
        
      end do 
      
      return
      end       
      
