! Upwind Method for solving Hyperbolic PDE
! Mustafa Khawari
! 16.10.1399

implicit none
real, parameter :: dx = 0.05, dt = 0.1, u = 0.1
integer, parameter :: Imax = 101, Nmax = 351
real*8 :: X (Imax), T (Imax, Nmax)
real :: cr, tim
integer :: i, n, h

open (6, file = "Hyper_1D_Upwind_EM.plt")
write (6, 10)

h = int (0.5 / dx) + 1
cr = u * dt / dx

do i = 1, Imax
  X (i) = (i-1) * dx
end do

T = 0.
do i = 1, Imax
  if (i <= h) then
    T (i, 1) = 200 * X (i)
  else if (i > h .and. i <= (2*h-1)) then
    T (i, 1) = 200 - 200 * X (i)
  else 
    T (i, 1) = 0.
  end if
end do

  

do n = 1, Nmax - 1
  
  tim = ((n-1) * dt)
  write (6, 11) tim
  
  do i = 1, Imax 
    if (i == 1 ) then
      T (i, n+1) = T (i, n)
    else
      T (i, n+1) = T (i, n) - cr * (T (i, n) - T (i-1, n))
    end if
    write (6, 12) X (i), T (i, n)
  end do
  
end do

10    format ( 'VARIABLES=X  C' )
11    format ( 'ZONE T="' , f6.2 , '"' )
12    format ( f8.2 , f16.3 )
end