The Regular Value (Preimage) Theorem

Fix $p \in f^{-1}(q)$. Since $df_p: T_pM \to T_qN$ is surjective and $\dim T_pM = m \ge n = \dim T_qN$, its kernel $K = \ker df_p$ has dimension $m-n$. Choose a linear complement $W$ so $T_pM = K \oplus W$ with $\dim W = n$; then $df_p|_W : W \to T_qN$ is an isomorphism.

Work in local coordinates centered at $p$ and $q$ (via charts), so $f$ becomes a smooth map $\mathbb{R}^m \to \mathbb{R}^n$ with $f(0)=0$ and $df_0$ surjective. Reorder coordinates $(x,y) \in \mathbb{R}^{m-n}\times\mathbb{R}^n$ so that $\partial f/\partial y$ at $0$ is the invertible $n\times n$ block (possible since $df_0$ has rank $n$).

Define $\Phi(x,y) = (x, f(x,y))$. Then $d\Phi_0 = \begin{pmatrix} I & 0 \\ \partial f/\partial x & \partial f/\partial y \end{pmatrix}$ has $\det d\Phi_0 = \det(\partial f/\partial y) \ne 0$, so $\Phi$ is a local diffeomorphism by the Inverse Function Theorem. In the new coordinates $(x,y') = \Phi(x,y)$, the equation $f=q$ (i.e. $f=0$) becomes exactly $y'=0$.

So near $p$, $f^{-1}(q)$ is the set $\{y'=0\}$, which in these coordinates is literally an $(m-n)$-dimensional coordinate slice — a smooth chart for $f^{-1}(q)$. Since $p$ was arbitrary, every point of $f^{-1}(q)$ has such a chart, and the transition maps between these charts are restrictions of the (smooth) transition maps of $M$, hence smooth. Therefore $f^{-1}(q)$ is a smooth $(m-n)$-dimensional submanifold.