! Crank Nicolson Method for solving Parabolic PDE
! Mustafa Khawari
! 9.11.99

implicit none
  real, parameter :: dt = 0.1, aa = 0.01
  integer, parameter :: Imax = 21, Nmax = 500
  real :: X (Imax), T (Imax, Nmax), D (Imax-2), XX (Imax-2)
  real :: tim, dd, dx, XTot, a, b, c
  integer :: i, n, h

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

XTot = 1.
dx = XTot / (Imax-1)
h = int (Imax/2) + 1
dd = (aa * dt) / (dx**2)

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
    T (i, 1) = 200 - 200 * X (i)
  end if
end do

a =  dd/2
b = -(1 + dd)
c =  dd/2

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

  D (1) = -dd/2 * T (1, n) + (dd-1) * T (2, n) - dd/2 * T (3, n) - dd/2 * T (1, n+1)
  D (Imax-2) = -dd/2 * T (Imax-2, n) + (dd-1) * T (Imax-1, n) - dd/2 * T (Imax, n) - dd/2 * T (Imax, n+1)
  do i = 2, Imax-3
    D (i) = -dd/2 * T (i, n) + (dd-1) * T (i+1, n) - dd/2 * 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.4 )
end
!=====================================================
subroutine TDMA (IM, a, b, c, D, X)
real :: g (IM-2),f (IM-2), D (IM-2), X (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