! Explicit Method for solving 2D Parabolic PDE
! Mustafa Khawari
! 10.11.99

module var2
  real, parameter :: dt = 2, a = 0.01
  integer, parameter :: Imax = 21, Jmax = 21, Nmax = 4000
  real :: T (Imax, Jmax, Nmax)
  real :: tim, dx, ddx, dy, ddy, Xtot, Ytot, x0, y0, x, y, rt
  integer :: i, j, n, npt
end
!---------------- The main of pogram ------------------------
use var2
open (6, file = "Para_2D_EM_.plt")

  npt = 5
  Xtot = 20.
  Ytot = 10.
  dx = Xtot / (Imax-1)
  dy = Ytot / (Jmax-1)
  x0 = Xtot/2
  y0 = Ytot/2
  ddx = (a * dt) / (dx**2)
  ddy = (a * dt) / (dy**2)

T = 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)
  end do
end do

do n = 1, Nmax - 1
  if (mod((n-1), npt) == 0) call plot2 ()
  
  do i = 1, Imax
    do j = 1, Jmax
      if ((i == 1 .or. i == Imax) .or. (j == 1 .or. j == Jmax)) then
       T (i, j, n+1) = T (i, j, n)
      else
        T (i, j, n+1) = T (i, j, n) + ddx * (T (i+1, j, n) - 2 * T (i, j, n) + T (i-1, j, n)) &
        & + ddy * (T (i, j+1, n) - 2 * T (i, j, n) + T (i, j-1, n))
      end if
    end do
  end do

  print *, n
end do
print *, "     d =", ddx + ddy
end
!======================================================================
subroutine plot2 ()

use var2
write (6, 10)

tim = (n-1) * dt
write (6, 11) tim, Imax, Jmax
do i = 1, Imax
  x = (i-1) * dx
  do j = 1, Jmax 
    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= ' , i4 , ' J= ' , i4 , ' F = POINT ' )
12    format ( 2f8.3, f20.6 )

return
end