! The Line Gauss-Seidel Iteration Method for solving 2D Elliptic PDE
! Mustafa Khawari
! 14.11.99

module var2
  integer, Parameter :: Kmax = 10000, Imax = 11, Jmax = 11, Par = 500
  real, parameter ::  Xtot = 10, Ytot = 10
  real :: T (Imax, Jmax, Kmax), D (Par), XX (par)
  integer :: i, j, k
  real :: beta, err, dx, dy, a, b, c
end 

use var2
open (8, file = "Elliptic_2D_LGaussSeidel.plt")

dx = Xtot / (Imax-1)
dy = Ytot / (Jmax-1)
beta = dx/dy

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

D = 0.
  a = -1
  b = 2 * (1 + beta**2)
  c = -1
do k = 1, Kmax-1

  do j = 1, Jmax
    T (Imax, j, k+1) = 100.
  end do
  
  do j = 2, Jmax-1
    D (1) = beta ** 2 * (T (2, j+1, k) + T (2, j-1, k+1)) + T (1, j, k+1)
    D (Imax-2) = beta ** 2 * (T (Imax-1, j+1, k) + T (Imax-1, j-1, k+1)) + T (Imax, j, k+1)
    do i = 2, Imax-3
      D (i) = beta ** 2 * (T (i+1, j+1, k) + T (i+1, j-1, k+1))
    end do
    
    call TDMA (Imax, a, b, c, D, XX)
    do i = 2, Imax - 1
      T (i, j, k+1) = XX (i-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 ()
end
!======================================================================
subroutine plot2 ()
use var2
write (8, 10)
      
write (8, 11) k, Imax, Jmax
do i = 1, Imax
  x = (i-1) * dx
  do j = 1, Jmax
    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