program jacobi
   implicit none
!hpf$ nosequence
   integer, parameter :: n = 20, iterations = 200
   integer count, i, j
   integer, dimension(n, n) :: field, work_space
!hpf$ align work_space(:,:) with field
!hpf$ distribute (*, block) :: field

!initialisation of the field
   field = 3
   field(ubound(field, dim=1),:) = 70
   field(lbound(field, dim=1),:) = -1
   field(:, ubound(field, dim=2)) = 40
   field(:, lbound(field, dim=2)) = 40

!print the matrix before iterations
   write (6, '(1x)')
   do i = 1, size(field, dim=1)
      write (6, '(1x, 20i3)') field(i, :)
   end do

!Jacobi iteration
   work_space = field
   do count = 1, iterations/2
      forall (i = 2:size(field, dim=1)-1, j=2:size(field, dim=2)-1)
         work_space(i, j) = (field(i, j-1) + field(i-1, j)             &
                             + field(i+1, j) + field(i, j+1))/4
         field(i,j) = (work_space(i, j-1) + work_space(i-1, j)         &
                       + work_space(i+1, j) + work_space(i, j+1))/4
      end forall
   end do

!print the matrix after iterations
   write (6, '(1x)')
   do i = 1, size(field, dim=1)
      write (6, '(1x, 20i3)') field(i, :)
   end do

end program jacobi