get_indices_in_int_array Subroutine

    subroutine get_indices_in_int_array(array, val, iloc, jloc)
        !! return the indices of all the `val` elements in 2d `array`.
        integer,dimension(:,:),intent(in) :: array
        integer,intent(in) :: val
        integer,dimension(:),allocatable,intent(out) :: iloc, jloc
        integer :: k !! counter
        ! do this by creating index matrices, and using mask to get the ones we want
        iloc = pack(spread([(k, k = 1, size(array,1))], dim=2, ncopies=size(array,2)), mask=array==val)
        jloc = pack(spread([(k, k = 1, size(array,2))], dim=1, ncopies=size(array,1)), mask=array==val)
    end subroutine get_indices_in_int_array