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 NOWAIT
     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 NOWAIT
    return
  end subroutine smooth_y

end module mysubs


program orphan3
  !! flawed !!

  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)
  enddo
  !$OMP END PARALLEL

  call check(a)  ! check result

end program orphan3