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

module var2
  integer, Parameter :: Kmax = 10000, Imax = 11, Jmax = 11 
  real, parameter ::  Xtot = 10, Ytot = 10
  real :: T (Imax, Jmax, Kmax)
  integer :: i, j, k
  real :: beta, err, dx, dy
end 

use var2
open (8, file = "Elliptic_2D_PGaussSeidel.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

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)
      else
        T (i, j, k+1) = 1/(2 * (1 + beta**2)) * (T (i+1, j, k) + T (i-1, j ,k+1) + beta**2 * (T (i, j+1, k) + T (i, j-1, k+1)))
      end if
    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