! Crank Nicolson Method for solving Hyperbolic PDE
! Mustafa Khawari
! 9.11.1399

implicit none
  real, parameter :: dx = 0.05, dt = 0.2, u = 0.1
  integer, parameter :: Imax = 201, Nmax = 200
  real*8 :: X (Imax), T (Imax, Nmax), D (Imax-2), XX (Imax-2)
  real :: cr, tim, a, b, c
  integer :: i, n, h

open (6, file = "Hyper_1D_CrankNicolson_IM.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

a = -cr / 4
b = 1
c = cr / 4

do n = 1, Nmax - 1 
  tim = ((n-1) * dt)
  write (6, 11) tim

  D (1) = cr/4 * T (1, n) + T (2, n) - cr/4 * T (3, n) + cr/4 * T (1, n+1)
  D (Imax-2) = cr/4 * T (Imax-2, n) + T (Imax-1, n) - cr/4 * T (Imax, n) - cr/4 * T (Imax, n+1)
  do i = 2, Imax-3
    D (i) = cr/4 * T (i, n) + T (i+1, n) - cr/4 * T (i+2, n)
  end do
  
  call TDMA (Imax, a, b, c, D, XX)
  
  do i = 1, Imax
    if (i == 1 .or. i == Imax) then
      T (i, n+1) = 0.
    else
      T (i, n+1) = XX (i-1)
    end if
    write (6, 12) X (i), T (i, n)
  end do 
  
  print *, n
end do

10    format ( 'VARIABLES=X  C' )
11    format ( 'ZONE T="' , f6.2 , '"' )
12    format ( f8.2 , f16.3 )
end
!==================================================
subroutine TDMA (IM, a, b, c, D, X)
real*8 :: g (IM-2), f (IM-2), X (IM-2), D (IM-2)

g (1) = D (1) / b
f (1) = c / b
do i = 2, IM-2
  rt = (b - a * f (i-1))
  g (i) = (d (i) - a * g (i-1)) / rt
  if (i <= IM-3) f (i) = c / rt
end do

X (IM-2) = g (IM-2)
do i = IM-3, 1, -1
  X (i) = g (i) - f (i) * X (i+1)
end do
return
end