! ADI Method for solving 2D Parabolic PDE
! Mustafa Khawari
! 14.11.99

module var2
  real, parameter :: dt = 1, aa = 0.01
  integer, parameter :: Imax = 21, Jmax = 21, Nmax = 1000, Par = 1000
  real :: T (Imax, Jmax, Nmax), TP (Imax, Jmax), D (Par), XX (Par)
  real :: tim, dx, ddx, dy, ddy, Xtot, Ytot, x0, y0, x, y, rt, a, b, c
  integer :: i, j, n, npt
end
!---------------- The main of pogram ------------------------
use var2
open (6, file = "Para_2D_ADI.plt")
  
  npt = 5
  Xtot = 20.
  Ytot = 10.
  dx = Xtot / (Imax-1)
  dy = Ytot / (Jmax-1)
  x0 = Xtot/2
  y0 = Ytot/2
  ddx = (aa * dt) / (dx**2)
  ddy = (aa * dt) / (dy**2)
  
T = 0.
TP = 0.
do i = 1, Imax
  x = (i-1) * dx
  do j = 1, Jmax
    y = (j-1) * dy
    rt = (x-x0) ** 2 + (y-y0) ** 2
    T (i, j, 1) =  8 * exp(-rt/2)
    TP (i, j) = T (i, j, 1)
  end do
end do


do n = 1, Nmax - 1 
! The first step
D = 0. 
  a = -ddx
  b = 1 + 2*ddx
  c = -ddx
  
  do j = 2, Jmax-1
    
    D (1) = ddy * T (2, j-1, n) + (1-2*ddy) * T (2, j, n) + ddy * T (2, j+1, n) + ddx * TP (1, j)
    D (Imax-2) = ddy * T (Imax-1, j-1, n) + (1-2*ddy) * T (Imax-1, j, n) + ddy * T (Imax-1, j+1, n) + ddx * TP (Imax, j)
    do i = 2, Imax-3
      D (i) = ddy * T (i+1, j-1, n) + (1-2*ddy) * T (i+1, j, n) + ddy * T (i+1, j+1, n)
    end do
    
    call TDMA (Imax, a, b, c, D, XX)
    do i = 2, Imax - 1
      TP (i, j) = XX (i-1)
    end do  
  end do
    
! The second step 
D = 0.  
  a = -ddy
  b = 1 + 2*ddy
  c = -ddy
  
  do i = 2, Imax-1
    
    D (1) = ddx * TP (i-1, 2) + (1-2*ddx) * TP (i, 2) + ddx * TP (i+1, 2) + ddy * T (i, 1, n+1)
    D (Jmax-2) = ddx * TP (i-1, Jmax-1) + (1-2*ddx) * TP (i, Jmax-1) + ddx * TP (i+1, Jmax-1) + ddy * T (i, Jmax, n+1)
    do j = 2, Jmax-3
      D (j) = ddx * TP (i-1, j+1) + (1-2*ddx) * TP (i, j+1) + ddx * TP (i+1, j+1)
    end do
    
    call TDMA (Jmax, a, b, c, D, XX)
    do j = 2, Jmax - 1
      T (i, j, n+1) = XX (j-1)
    end do
  end do
  
  if (mod (n-1, npt) == 0)  call plot2 (Imax, Jmax)
print *, n
end do 
  
print *, "     d =", ddx + ddy
end
!======================================================================
subroutine plot2 (IM, JM)
use var2
write (6, 10)
      
tim = (n-1) * dt
write (6, 11) tim, IM, JM
do i = 1, IM
  x = (i-1) * dx
  do j = 1, JM 
    y = (j-1) * dy
    write (6, 12) x, y, T (i, j, n)
  end do
end do

10    format ( 'VARIABLES=X  Y  C' )
11    format ( 'ZONE T="', f6.1 , '" I= ' , i6 , ' J= ' , i6 , ' F = POINT ' )
12    format ( 2f10.3, f20.6 )
return
end
!====================================================================
subroutine TDMA (IM, a, b, c, D, XX)
real :: g (IM), f (IM), XX (IM), D (IM)

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

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