!********************module*******************************************
module para
  implicit  none
  real(8),parameter :: x0 = 0.d0,   &
                     & xm = 1.d0,   &
                     & y0 = 0.d0,   &
                     & ym = 1.d0,   &
                     & omg = 1.4d0, &
                     & eps = 1.d-8
  integer,parameter :: m = 41, &
                     & itrmax = 100000
  integer :: i,j,k,LL
  real(8) :: dx,dy,pi
end module para
!****************main*************************************************
program poisson_equation
  use para
  implicit none
  real(8),dimension(0:m) :: x,y
  real(8),dimension(0:m,0:m) :: u,f
  
  call calc_param(x,y)
  call bcon(u)
  call icon(u)
  call calc_Rhand(x,y,f)
  call sor(u,f)
  
  stop
  
end program poisson_equation
!**********calculate parameter****************************************
subroutine calc_param(x,y)
  use para
  implicit none
  real(8),dimension(0:m) :: x,y
  
  x(0) = x0
  x(m) = xm
  y(0) = y0
  y(m) = ym
  
  dx = (xm-x0)/dble(m-1)
  dy = (ym-y0)/dble(m-1)
  
  do i = 1,m-1
      x(i) = dx*dble(m-1)
      y(i) = dy*dble(m-1)
  end do
  
  pi = 4.d0*datan(1.d0)
  
  return
  
end subroutine calc_param
!***************boundary condition************************************
subroutine bcon(u)
  use para
  implicit none
  real(8),dimension(0:m,0:m) :: u
  
  do i = 0,m
    u(i,0) = 0.d0
    u(i,m) = 0.d0
    u(0,i) = 0.d0
    u(m,i) = 0.d0
  end do
  
  return
  
end subroutine bcon
!***************initial condition*************************************
subroutine icon(u)
  use para
  implicit none
  real(8),dimension(0:m,0:m) :: u
  
  do j = 1,m-1
    do i = 1,m-1
      u(i,j) = 0.d0
    end do
  end do
  
  return
  
end subroutine icon
!*************calculate right side************************************
subroutine calc_Rhand(x,y,f)
  use para
  implicit none
  real(8),dimension(0:m)     :: x,y
  real(8),dimension(0:m,0:m) :: f
  
  do j = 0,m
    do i = 0,m
      f(i,j) = -5.d0*pi*pi*dsin(pi*x(i))*dsin( 2.d0*pi*y(j))
    end do
  end do
  
  return
  
end subroutine calc_Rhand
!*************solver SOR method***************************************
subroutine sor(u,f)
  use para
  implicit none
  integer :: n
  real(8),dimension(0:m,0:m) :: u,f,uk,d
  real(8) :: max
  
  do j = 0,m
    do i = 0,m
      uk(i,j) = 0.d0
    end do
  end do
  
  n = 0
  
  do LL = 1,itrmax
     
     n = n+1
     
     do j = 1,m-1
       do k = 1,m-1
         uk(k,j) = u(k,j)
       end do
     end do
     
     do j = 1, m-1
       do i = 1,m-1
         u(i,j) = uk(i,j)+omg/4.d0*(u(i-1,j)+u(i,j-1) &
                & +uk(i+1,j)+uk(i,j+1)-dx*dx*f(i,j)-4.0*uk(i,j))
                 
         d(i,j) = dabs(u(i,j)-uk(i,j))**2
       end do
     end do
     
     max = 0.0
     
     do j = 1,m-1
       do k = 1,m-1
         if(max < d(k,j)) then
           max = d(k,j)
         end if
       end do
     end do
        
     if(max < eps) exit
     
  end do
  
  return
  
end subroutine sor
!*********************************************************************