```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                           :: i, j, 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
```