module mod_laplace
use mod_util
contains

subroutine laplace_fortran(M,N,LDIM,x,r,y,omega,niter)
!*** arguments
integer, intent(in)::M,N,LDIM,niter      ! sizes
real, intent(in)::omega              ! relaxation parameter
real,dimension(LDIM,N),intent(inout)::x ! solution
real,dimension(LDIM,N),intent(in)::r    ! right-hand side
real,dimension(LDIM,N),intent(inout)::y ! work
!*** local
integer:: k
do k=1,niter
   call laplace_onepass(M,N,LDIM,y,x,r,omega)
   call laplace_onepass(M,N,LDIM,x,y,r,omega)
enddo
end subroutine laplace_fortran

subroutine laplace_onepass(M,N,LDIM,y,x,r,omega)
implicit none
!*** purpose
! one step for scaled Laplace equation on M by N grid
!*** arguments
integer, intent(in)::M,N,LDIM      ! sizes
real, intent(in)::omega              ! relaxation parameter
real,dimension(LDIM,N),intent(out)::y ! solution
real,dimension(LDIM,N),intent(in)::x ! solution
real,dimension(LDIM,N),intent(in)::r    ! right-hand side
!** local
integer i,j,it
real:: res
!***
! print *,'laplace start'
    do j=2,N-1
       do i=2,M-1
          res=r(i,j)-(x(i,j)-0.25*(x(i-1,j)+x(i+1,j)+x(i,j-1)+x(i,j+1)))
          y(i,j)=x(i,j)+omega*res
          !print *,"onerow:",i,j,Y(i,j),R(i,j),X(i,j),res,omega
          !print *,X(i-1,j-1),X(i-1,j),X(i-1,j+1)
          !print *,X(i,j-1),X(i,j),X(i,j+1)
          !print *,X(i+1,j-1),X(i+1,j),X(i+1,j+1)
       enddo
    enddo
! print *,'laplace end'
end subroutine laplace_onepass
end module mod_laplace