Module bh01
  type Mesh
    integer :: sizex, sizey
    real, pointer,dimension(:,:)::grid
  end type Mesh
end module bh01

program laplace
use bh01
type(Mesh), pointer :: mygrid

interface
  subroutine allocategrid(mymesh, x, y)
    use bh01
    type(mesh),pointer :: mymesh
    integer :: x , y
  end subroutine allocategrid
end interface

interface
  function update(z,mesh, x, y)
    integer  x, y
    real, dimension(x,y) :: z
    real, dimension(x,y) :: mesh
  end function update
end interface

real, parameter :: Epsilon = 0.001
integer, parameter :: R = 128, C = 64
real :: error = 100.0, bsum
integer :: x, y, niter, kount = 1
integer, pointer, dimension(:,:) :: I, J
real, pointer, dimension(:,:) :: Z

allocate(Z(R,C))
allocate(I(R,C))
allocate(J(R,C))

I = SPREAD( [1:R], 2, C )     !Coordinates
J = SPREAD( [1:C], 1, R )

!Z is a solution of the difference equation
Z = -1.4 + 0.12*I - 0.23*J + 0.001*(I**2 - J**2)

call allocategrid(mygrid, R, C)

mygrid%grid = Z    ! Install boundary values

print *, "Solution range: ", minval(Z), maxval(Z)
bsum = SUM(Z(1,:)) + SUM(Z(R,:)) + SUM(Z(2:R-1,1)) + SUM(Z(2:R-1,C))

!Intialize interior to average on boundary
!The "-4" prevents counting corners twice
mygrid%grid(2:R-1, 2:C-1) = bsum/(2*R + 2*C - 4)

do while (error > Epsilon)
  kount = kount + 1
  error = update(Z(:,:), mygrid%grid(:,:), mygrid%sizex, mygrid%sizey)
if (mod(kount,1000).eq. 0) print *, 'error = ', error, kount
end do

  print *, 'error = ', error, kount
end program

!-----------------subroutine allocate-------------------------
subroutine allocategrid(mygrid, x, y)
use bh01
type(mesh),pointer :: mygrid
integer :: x , y

allocate(mygrid)

mygrid%sizex = x
mygrid%sizey = y

allocate(mygrid%grid(x,y))

return
end



!-----------------subroutine update-------------------------
function update(Z,grid, i, j)
real, dimension(i,j) :: grid
real, dimension(i,j) :: Z
integer :: i, j

forall(k=2:i-1,l=2:j-1) grid(k,l) = (grid(k-1,l)+ grid(k+1,l) + grid(k,l-1) + grid(k,l+1))/4.0

update = MAXVAL(abs(Z - grid))

return
end
