!!
!! Liebmann's method to compute heat transfer across a 2-D surface
!! J. Overbey 8/27/08, updated for OpenMP 5/3/12
!!
!! Use liebmann-viz.sh in the main project directory to visualize the resulting
!! table using gnuplot
!!
!! This version is based on the code in src-liebmann-4-openacc but uses OpenMP
!! rather than OpenACC.
!!
program liebmann_example
    implicit none

    integer, parameter :: SIZE = 4096
    integer, parameter :: INTERIOR_SIZE = SIZE - 2
    integer, parameter :: OUTPUT_SIZE = 128
    real, parameter    :: BOUNDARY_VALUE = 5.0
    real, parameter    :: EPSILON = 0.001

    call main()

contains

subroutine main()
    real :: surface(SIZE, SIZE)
    integer :: i, j
    integer :: num_threads
    integer :: start_time, end_time, count_rate

    call system_clock(start_time, count_rate)
    call liebmann(surface, num_threads)
    call system_clock(end_time)

    do i = 1, SIZE, SIZE/OUTPUT_SIZE
        do j = 1, SIZE-1, SIZE/OUTPUT_SIZE
            write (*,'(F4.2" ")',advance="no") surface(i, j)
        end do
        write (*,'(F4.2" ")') surface(i, SIZE)
    end do

    print *, "Threads: ", num_threads
    print *, "Elapsed Time (seconds):", real(end_time - start_time) / real(count_rate)
end subroutine

subroutine liebmann(surface, num_threads)
    real, dimension(SIZE, SIZE), intent(out) :: surface
    integer, intent(out) :: num_threads
    real, dimension(SIZE, SIZE) :: prev, next
    real :: delta
    integer :: i, j
    logical :: done

    done = .false.

    !$omp parallel
    !$omp master
    num_threads = omp_get_num_threads()
    !$omp end master
    !$omp end parallel

    call init_with_boundaries(prev)
    call init_with_boundaries(next)

    done = .false.
    !$omp parallel
    do while (.not. done)
        !$omp do schedule(static) private(i)
        do j = 2, SIZE-1
            do i = 2, SIZE-1
                next(i,j) = &
                    (prev(i-1, j) + &
                     prev(i+1, j) + &
                     prev(i, j-1) + &
                     prev(i, j+1)) / 4.0
            end do
        end do
        !$omp end do
        delta = 0.0
        !$omp do schedule(static) private(i) reduction(max:delta)
        do j = 2, SIZE-1
            do i = 2, SIZE-1
                prev(i,j) = &
                    (next(i-1, j) + &
                     next(i+1, j) + &
                     next(i, j-1) + &
                     next(i, j+1)) / 4.0
                 delta = max(delta, abs(prev(i,j)-next(i,j)))
            end do
        end do
        !$omp end do
        !$omp single
        if (delta < EPSILON) then
            done = .true.
        end if
        !$omp end single
    end do
    !$omp end parallel
    surface = prev
end subroutine

subroutine init_with_boundaries(surface)
    real, dimension(SIZE, SIZE), intent(out) :: surface

    surface          = 0.0
    surface(1, :)    = BOUNDARY_VALUE
    surface(SIZE, :) = BOUNDARY_VALUE
    surface(:, 1)    = BOUNDARY_VALUE
    surface(:, SIZE) = BOUNDARY_VALUE
end subroutine

end program