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 orphan2 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) enddo !$OMP END PARALLEL call check(a) ! check result end program orphan2