When passing a multidimensional array as an argument to a subprogram, is is even more critical that the dimensions are known.

With a one-dimensional array, the subprogram must know the size of the array in order to avoid out-of-bounds errors.

With a multidimensional array, we still have the possibility of out-of-bounds errors. In addition, if the correct dimensions are not known to the subprogram, it will not be able to compute the correct location of an element from the subscripts! Look again at the memory map of a small 2D array:

        double precision :: matrix(2, 3)
        
        Memory map:
        
        2000    matrix(1,1)
        2008    matrix(2,1)
        2016    matrix(1,2)
        2024    matrix(2,2)
        2032    matrix(1,3)
        2040    matrix(2,3)
        

The memory address of an element is computed as:

address(r, c) = base + (total-rows * (c-1) + (r-1)) * sizeof(type)

For an array of double precision values based at address 2000:

address(r,c) = 2000 + (2 * (c-1) + (r-1)) * 8

address(2,3) = 2000 + (2 * 2 + 1) * 8 = 2040

Since the calculation of an address involves the total number of rows in the matrix, we must at least get this right just to be able access the array elements correctly. Getting the number of columns right is also essential to prevent out-of-bounds problems.

        subroutine print_matrix(matrix, rows, cols)
            ! Disable implicit declarations (i-n rule)
            implicit none
            
            ! Dummy variables
            integer, intent(in) :: rows, cols
            double precision, intent(in) :: matrix(rows, cols)
        
            ! Local variables
            integer :: r
        
            do r = 1, rows
                print *, matrix(r, 1:cols)
            enddo
        end subroutine