program omp_red_black

    use omp_lib

    implicit none

    integer, parameter                  :: DSIZE_X = 10, DSIZE_Y= 5000, ITERNUM = 10000
    double precision, dimension(0:DSIZE_X-1, 0:DSIZE_Y-1)   :: grid
    double precision                    :: sum = 0.0
    double precision                    :: start_time, final_time
    integer                             :: i, j, n

    ! starting the clock
    start_time = omp_get_wtime()

    ! Parallelized this loop.
    ! Parallelization is done in the following way to benefit
    ! from Fortran coloumn major order array idexing
    !$omp parallel do private(i) schedule(static)
    do j = 0, DSIZE_Y - 1
        do i = 0, DSIZE_X - 1
            grid(i, j) = mod((i * j), 10)
        end do
    end do
    !$omp end parallel do


    ! Parallelizing the outer loop is not possible.
    ! Each iteration depends on the value of previous iteration
    ! so we will parallelize one level below

    do n = 1, ITERNUM
    !$omp parallel private(i, j)
        ! Update red points
        !$omp do schedule(static)
        do j = 1, DSIZE_Y - 2, 2
            do i = 1, DSIZE_X - 2, 2
                grid(i, j) = 0.25 * (grid(i-1, j) + grid(i+1, j) +&
                            grid(i, j-1) + grid(i, j+1))
            end do
        end do
        !$omp end do nowait

        !$omp do schedule(static)
        do j = 2, DSIZE_Y - 2, 2
            do i = 2, DSIZE_X - 2, 2
                grid(i, j) = 0.25 * (grid(i-1, j) + grid(i+1, j) +&
                            grid(i, j-1) + grid(i, j+1))
            end do
        end do
        !$omp end do

        ! Update black points

        !$omp do schedule(static)
        do j = 2, DSIZE_Y - 2, 2
            do i = 1, DSIZE_X - 2, 2
                grid(i, j) = 0.25 * (grid(i-1, j) + grid(i+1, j) +&
                            grid(i, j-1) + grid(i, j+1))
            end do
        end do
        !$omp end do nowait

        !$omp do schedule(static)
        do j = 1, DSIZE_Y - 2, 2
            do i = 2, DSIZE_X - 2, 2
                grid(i, j) = 0.25 * (grid(i-1, j) + grid(i+1, j) +&
                            grid(i, j-1) + grid(i, j+1))
            end do
        end do
        !$omp end do
    !$omp end parallel
    end do

    !$omp parallel do private(i) schedule(static) reduction(+:sum)
    do j = 0, DSIZE_Y - 1
        do i = 0, DSIZE_X - 1
            sum = sum + grid(i, j)
        end do
    end do
    !$omp end parallel do

    ! Calculating total execution time
    final_time = omp_get_wtime() - start_time;

    write (*, '(a,f14.2)') 'Sum: ', sum
    write (*, '(a,f14.6)') 'Total time:  ', final_time
    stop
end program omp_red_black
