module mysubs
contains
subroutine init(a)
! initial values for a: a(i,j) = i-j
use kinds
implicit none
real(kind=REAL8), dimension(:,:) :: a
integer :: n, i, j
n = size(a, 1)
do i=1,n
do j=1,n
a(i,j) = i - j
enddo
enddo
return
end subroutine init
subroutine check(a)
! check results
! here: print total sum
use kinds
implicit none
real(kind=REAL8), dimension(:,:) :: a
real(kind=REAL8) :: sum1, sum2
integer :: n
n = size(a, 1)
sum1 = sum(a)
sum2 = 0.0
print *, 'sum is : ', sum1
print *, 'and should be: ', sum2
return
end subroutine check
function f(x, work)
! time consuming way of setting f = x + eps
use kinds
implicit none
real(kind=REAL8) :: f, x
integer :: work
integer :: i
real(kind=REAL8) :: t1, t2
t1 = x
do i= 1, work
t1 = cos(t1)
end do
t2 = x
do i= 1, work
t2 = sin(t2)
end do
f = x + t1 + t2 - 0.7390851332151607
end function f
subroutine smooth_x(a, b, work)
! smoothes a in x direction
use kinds
implicit none
real(kind=REAL8), dimension(:,:) :: a, b
integer :: work
integer :: i, j, n
n = size(a, 1)
!$OMP DO
do i= 2, n-1
do j=1, n
b(j,i) = (a(i-1,j) + f(a(i,j), work) + a(i+1,j))/3.0
end do
end do
!$OMP END DO
return
end subroutine smooth_x
subroutine smooth_y(a, b, work)
! smoothes a in y direction
use kinds
implicit none
real(kind=REAL8), dimension(:,:) :: a, b
integer :: work
integer :: i, j, n
n = size(a, 1)
!$OMP DO
do i= 1, n
do j=2, n-1
b(j,i) = (a(i,j-1) + f(a(i,j), work) + a(i,j+1))/3.0
end do
end do
!$OMP END DO
return
end subroutine smooth_y
end module mysubs
program orphan4
use kinds
use mysubs
implicit none
integer, parameter :: n = 20
integer, parameter :: tmax = 10
integer, parameter :: work = 10000
real(kind=REAL8), dimension(n,n) :: a, b
integer :: t
call init(a) ! initial values for a
b = a ! necessary for boundary values
!$OMP PARALLEL
do t=1, tmax
call smooth_x(a, b, work)
call smooth_y(b, a, work)
!$OMP SINGLE
print*, 't = ', t
!$OMP END SINGLE
enddo
!$OMP END PARALLEL
call check(a) ! check result
end program orphan4