! The Alternating Direction Implicit (ADI) Method for solving 2D Elliptic PDE
! Mustafa Khawari
! 15.11.99

module var2
  integer, parameter :: Imax = 11, Jmax = 11, Kmax = 10000, Par = 500
  real, parameter :: Xtot = 10, Ytot = 10
  real :: T (Imax, Jmax, Kmax), TP (Imax, Jmax), D (Par), XX (Par)
  real ::  a, b, c, beta, err
  integer :: i, j, k
end
!---------------- The main of pogram ------------------------
use var2
open (8, file = "Ellipt_2D_ADI.plt")
  
dx = Xtot / (Imax-1)
dy = Ytot / (Jmax-1)
beta = dx/dy

T = 0.
TP = 0.
do j = 1, Jmax
  T (Imax, j, 1) = 100.
  TP (Imax, j) = 100.
end do

do k = 1, Kmax - 1 

  do i = 1, Imax
    do j = 1, Jmax
      if ((i == 1 .or. i == Imax) .or. (j == 1 .or. j == Jmax)) then
        T (i, j, k+1) = T (i, j, k)
        TP (i, j) = T (i, j, k)
      end if
    end do
  end do

  ! The first step
  D = 0. 
  a = 1
  b = -2 * (1 + beta**2)
  c = 1
 
  do j = 2, Jmax-1
    
    D (1) = - beta**2 * (TP (2, j-1) + T (2, j+1, k)) - TP (1, j)
    D (Imax-2) = - beta**2 * (TP (Imax-1, j-1) + T (Imax-1, j+1, k)) - TP (Imax, j)
    do i = 2, Imax-3
      D (i) = - beta**2 * (TP (i+1, j-1) + T (i+1, j+1, k))
    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 = beta**2
  b = -2 * (1 + beta**2)
  c = beta**2
  
  do i = 2, Imax-1
    
    D (1) = - (T (i-1, 2, k+1) + TP (i+1, 2)) - beta**2 * T (i, 1, k+1)
    D (Jmax-2) = - (T (i-1, Jmax-1, k+1) + TP (i+1, Jmax-1)) - beta**2 * T (i, Jmax, k+1)
    do j = 2, Jmax-3
      D (j) = - (T (i-1, j+1, k+1) + TP (i+1, j+1)) 
    end do
    
    call TDMA (Jmax, a, b, c, D, XX)
    do j = 2, Jmax - 1
      T (i, j, k+1) = XX (j-1)
    end do
  end do
  

  err = abs (T (Imax-1, Jmax-1, k+1) - T (Imax-1, Jmax-1, k))
  if (err <= 0.000001) exit
print *, k
end do 
call plot2 (Imax, Jmax, dx, dy)
end
!======================================================================
subroutine plot2 (IM, JM, dx, dy)
use var2
write (8, 10)    

write (8, 11) k, IM, JM
do i = 1, IM
  x = (i-1) * dx
  do j = 1, JM 
    y = (j-1) * dy
    write (8, 12) x, y, T (i, j, k)
  end do
end do

10    format ( 'VARIABLES=X  Y  C' )
11    format ( 'ZONE T="', i6 , '" I= ' , i6 , ' J= ' , i6 , ' F = POINT ' )
12    format ( 2f10.3, f20.6 )
return
end
!====================================================================
subroutine TDMA (IM, a, b, c, D, X)
real :: g (IM), f (IM), X (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

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