program poisson
  real (kind(1d0)),pointer,dimension(:,:) :: g1,g2,p
  integer :: i,j,k,img, ht, n, niter
  integer, parameter :: SIZE=400

  ht=SIZE+2
  niter=40000

  allocate (g1(SIZE,ht),g2(SIZE,ht))

  g1=0.0
  g1(SIZE/8+1:7*SIZE/8,1)=1.0
  
! Top
  g2(:,1)=g1(:,1)
! Bottom
  g2(:,ht)=g1(:,ht)

  do n=1,niter

    do j=2,ht-1
! Left
      g2(1,j)=(g1(1,j-1)+g1(1,j+1)+g1(2,j))/3
! Body
      do i=2,SIZE-1
        g2(i,j)=0.25*(g1(i,j-1)+g1(i,j+1)+g1(i-1,j)+g1(i+1,j))
      enddo
! Right
      g2(SIZE,j)=(g1(SIZE,j-1)+g1(SIZE,j+1)+g1(SIZE-1,j))/3
    enddo

    p=>g1
    g1=>g2
    g2=>p

  enddo

  img=10
  open(unit=img,file='poisson.pgm')
  write(img,'(a2)')"P2"
  write(img,'(I0,x,I0,x,I0)')SIZE,SIZE,255
  do j=2,ht-1
    do i=1,SIZE
      write(img,'(I0)')int(255*g1(i,j))
    enddo
  enddo
end
