print_matrix.f90

  subroutine print_matrix(A, var_name)
    real, intent(in) :: A(:,:)
    character(len = *), intent(in) :: var_name
    integer :: i, M, N
    M = size(A, 1)
    N = size(A, 2)
    print *, var_name, '='
    do i = 1, M
      print *, (A(i, j), j = 1,N)
    end do
  end subroutine print_matrix