Fortran Gauss Siedel iterative solver using OpenMP does not converge

You could use ordered depend something like the below in order to create a parallel Gauss-Seidel solver. the example below solves a Poisson problem in 3D.

SUBROUTINE gauss_seidel_par(u,f,N,itermax)
   REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) :: u
   REAL(dp), DIMENSION(:,:,:), INTENT(IN) :: f
   INTEGER, INTENT(IN) :: N, itermax
   INTEGER :: i, j, k, l
   REAL(dp) :: factor1, factor2
   factor1 = 1.0_dp/6.0_dp              ! Coefficient of u(i,j,k)
   factor2 = (2.0_dp/REAL(N+1.0_dp))**2 ! Delta^2
   !$omp parallel
   DO l = 1,itermax                     ! run until maximum number of iterations is reached
   !$omp do schedule(static,1) ordered(2)
   DO k = 2,N-1
     DO j = 2,N-1
       !$omp ordered depend(sink: k-1,j) depend(sink: k,j-1)
       DO i = 2,N-1
         u(i,j,k) = factor1*(u(i+1,j,k)+u(i,j+1,k)+u(i,j,k+1)+u(i-1,j,k)+u(i,j-1,k)+u(i,j,k-1) + factor2*f(i,j,k))
       ENDDO
       !$omp ordered depend(source)
     ENDDO
   ENDDO
   ENDDO
   !$omp end parallel
   END SUBROUTINE gauss_seidel_par