You want to think of the Haar measure $d\mu(U)$ as a way of measuring uniformity in the group $U(N)$ of unitary $N\times N$ matrices.
To form your intuition, consider $N=1$. You then have $U=e^{i\phi}$, with $0<\phi\leq 2\pi$ and $d\mu(U)=d\phi$ measures the perimeter of the unit circle. This is a uniform measure, because $d(\phi+\phi_0)=d\phi$ for any fixed phase shift $\phi_0$. You could write the requirement of uniformity in the form $d\mu(UU_0)=d\mu(U)$, with $U_0=e^{i\phi_0}$ the unitary matrix corresponding to the phase shift $\phi_0$.
Once your intuition is formed for $N=1$, you simply generalize to $N>1$ using the same definition of uniformity, $d\mu(UU_0)=d\mu(U)$ for any fixed $U_0\in U(N)$. For orthogonal (or symplectic) matrices you use the same definition of uniformity, with $U_0$ now restricted to the orthogonal or symplectic subgroup of $U(N).$
To explicitly write down the Haar measure $d\mu(U)$ in terms of the matrix elements of $U$ is only easily done for a few small values of $N$. (In particular, there is no relationship to random directions of rows or columns, as Yemon Choi pointed out.) You typically do not need such explicit expressions, since integrals with the Haar measure can be evaluated by using only the definition of uniformity.
In response to the follow-up question: If you wish to evaluate Haar-measure integrals of polynomials of matrix elements of $U$, you can use the socalled Weingarten functions.
http://en.wikipedia.org/wiki/Weingarten_function
Here is a Mathematica program to generate these,
http://arxiv.org/abs/1109.4244
If you need an explicit expression for the Haar measure, the steps to take are the following:
1) parameterize your matrix $U$ in terms of a set of real parameters $\{x_i\}$.
2) calculate the metric tensor $m_{ij}$, defined by $\sum_{ij}|dU_{ij}|^2 = \sum_{ij}m_{ij}dx_i dx_j$
3) obtain the Haar measure by equating $d\mu(U) = ($Det $m)^{1/2}\prod_i dx_i$
This is the general recipe. In practice, for many parameterizations the answer is in the literature. In particular, for the Haar measure in Euler angle parameterizations see:
http://arxiv.org/abs/math-ph/0205016
http://www.cft.edu.pl/~karol/pdf/ZK94.pdf
Best Answer
The general formula in the circular ensembles follows from equation 2.10 of arXiv:cond-mat/9612179: $$P(H)\propto\text{det}\,(1-H)^{\tfrac{1}{2}\beta(N-2M+1-2/\beta)},\;\;N-2M\geq 0,\qquad\qquad(\ast)$$ where $H=AA^\dagger$ and $A$ is an $M\times M$ principal submatrix of an $N\times N$ unitary matrix $U$ (circular unitary ensemble, $\beta=2$) or unitary symmetric matrix (circular orthogonal ensemble, $\beta=1$) or unitary selfdual matrix (circular symplectic ensemble, $\beta=4$).
The case in the OP is $\beta=1$, in that case $H=AA^\dagger=AA^\ast$ (in the physics notation that $\ast$ is complex conjugation and $\dagger$ is Hermitian conjugate). The assumption that $N-2M\geq 0$ is needed, because if $M>N/2$ there will be an additional set of $M-N/2$ eigenvalues of $H$ that are pinned to unity.
The case that $U$ is real orthogonal (circular real ensemble, $H=AA^\dagger=AA^\top$) has an additional factor, $$P(H)\propto\text{det}\,H^{-1/2}\,\text{det}\,(1-H)^{\tfrac{1}{2}(N-2M-1)}.\qquad\qquad(\ast\ast)$$ That case is considered in an earlier MO posting. Note that $(\ast\ast)$ agrees with the corresponding result in the OP (the second equation), which gives $P(A)$ instead of $P(H)$.
To make contact with the physics literature: $U$ represents a scattering matrix, $A$ is the reflection matrix, and the eigenvalues of $1-H=1-AA^\dagger$ are the transmission eigenvalues $T_1,T_2,\ldots T_M$. Equation (2.10) in the cited paper gives the joint probability distribution of the $T_i$'s, which can equivalently be expressed as the equation $(\ast)$ above.