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(v)
    ! check results
    ! here: print total sum
    
    use kinds
    implicit none
    
    real(kind=REAL8), dimension(:)  :: v
    real(kind=REAL8)                :: sum1, sum2
    integer                         :: n
    
    n = size(v, 1)
    sum1 = sum(v)
    sum2 = n**3 + n**2
    
    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

end module mysubs


program loop1

  use kinds
  use mysubs
  implicit none

  integer, parameter                :: n = 30
  integer, parameter                :: work = 10000

  real(kind=REAL8), dimension(n,n)  :: a
  real(kind=REAL8), dimension(n)    :: v
  integer                           :: i, j
      
  v = 0
  call init(a)   ! initial values for a

  do i=1, n
     do j=1, n
        v(i) = v(i) + f(a(i,j), work)
     enddo
  enddo

  call check(v)  ! check result

end program loop1