Let $\Omega\subseteq\mathbb{R}^n$ be an open set, and let $u:\bar{\Omega}\rightarrow\mathbb{R}$ be a function of class $C^1(\bar{\Omega})$. Given $\lambda\in\mathbb{R}$, let $\Gamma_\lambda=\{x\in\Omega:\,u(x)=\lambda\}$.
Coarea formula: Suppose $|\nabla u|>0$ on $\bar{\Omega}$, and let $f\in L^1(\Omega)$. Then $$\int_\Omega f\,dx=\int_\mathbb{R}\int_{\Gamma_\lambda} \frac{f}{|\nabla u|}\,d\sigma\,d\lambda.$$
I present a proof that I think would work, and I want to ask about a (very important) step inside the proof. I know that coarea formula can be proved in a more general setting, but I am interested on solving this particular proof.
For each $p\in\Omega$, there is $i\in\{1,\ldots,n\}$ with $|u_{x_i}(p)|>0$. By continuity, there exists $r_p>0$ such that for all $x\in B(p,r_p)\subseteq\Omega$ we have $|u_{x_i}(x)|>0$. As $\Omega$ is Lindelöf, $\Omega=\cup_{j=1}^{\infty} B_j$, where $B_j$ is some of such balls. Partition of the unity: there exists $\psi_j\in C_c^{\infty}(B_j)$ with $0\leq\psi_j\leq 1$ and $\sum_{j=1}^{\infty}\psi_j(x)=1$ for all $x\in\Omega$. Define $f_j=f\,\psi_j$ on $\Omega$. From now on, fix one of the balls $B=B_j$ and for simplicity of notation assume that $|u_{x_n}|>0$ on $B$. By the implicit function theorem, $\Gamma_\lambda\cap B=\{x\in B:\,u(x)=\lambda\}=\{(x',\varphi(x',\lambda)):\,x'\in\tilde{B}\}$, where $\tilde{B}$ open in $\mathbb{R}^{n-1}$ and $\varphi(\cdot,\lambda)\in C^1(\tilde{B})$.
Question: can I assume that $\tilde{B}$ is chosen independently of $\lambda$? Do we have $\varphi(x',\cdot)\in C^1$?
On the one hand, $$\int_\mathbb{R}\int_{\Gamma_\lambda}f_j\,d\sigma\,d\lambda=\int_\mathbb{R}\int_{\Gamma_\lambda\cap B}f_j\,d\sigma\,d\lambda=\int_\mathbb{R}\int_{\tilde{B}}f_j(x',\varphi(x',\lambda))\sqrt{1+|\nabla_{x'}\varphi(x',\lambda)|^2}\,dx'\,d\lambda.$$ On the other hand, if the question has a positive answer,
\begin{equation} \begin{split} \int_\Omega f_j\,|\nabla u|\,dx &=\int_B f_j\,|\nabla u|\,dx \\ &\stackrel{\substack{\text{change}\\\text{variable}}}{=}\int_\mathbb{R}\int_{\tilde{B}}f_j(x',\varphi(x',\lambda))|\nabla u(x',\varphi(x',\lambda))|\,|\varphi_\lambda (x',\lambda)|\,dx'\,d\lambda. \end{split} \end{equation}
Finally use $u(x',\varphi(x',\lambda))=\lambda$ for all $x'\in \tilde{B}$, and derivating with respect to $x'$ and $\lambda$ (if the question has a positive answer), one arrives at the equality between the last two big expressions.
If the answer to the question were negative, are there any alternatives to make this proof work? I don't know, something like changing the balls by rectangles...
No, you cannot assume $\tilde B$ to be independent of $\lambda$. Consider the example $u(x)=x_n$ and $B=B^n(0,1)\subset\mathbb R^n$. Now $\Gamma_\lambda=\{\lambda\}\times B^{n-1}(0,\sqrt{1-\lambda^2})$. Your set $\tilde B$ is now the ball $B^{n-1}(0,\sqrt{1-\lambda^2})$, and its radius shrinks to zero as you approach the boundary of $B$ (when $\lambda\to\pm1$).
Here is a possible alternative approach.
Take any point $x\in\Omega$. Since $\nabla u(x)\neq0$, there is $i\in\{1,\dots,n\}$ so that $u_i(x):=\partial u(x)/\partial x_i\neq0$. By continuity there is a ball $B^n(x,\epsilon)$ so that $u_i\neq0$ in $B^n(x,\epsilon)$ and $\bar B^n(x,\epsilon)\subset\Omega$.
For notational simplicity assume that $i=n$ and $u_n(x)>0$. Denote $x'=(x_1,\dots,x_{n-1})\in\mathbb R^{n-1}$. Suppose $\delta<\epsilon$; we will fix the value later. Consider the tube $T=B^n(x,\epsilon)\cap (B^{n-1}(x',\delta)\times\mathbb R)$. The tube has two "caps": $C_\pm=\{x\in\partial B^n(x,\epsilon)\cap\partial T;\pm x_n>0\}$.
Since $u_n>0$, the function is supposed to be bigger on $C_+$ than $C_-$, but this is only true if $\delta$ is small enough. We need an auxiliary claim: For sufficiently small $\delta>0$ we have $\inf_{C_+}u>u(x)>\sup_{C_-}u$. Proof: In the limit $\delta\to0$ the infimum on $C_+$ goes to $u(x+\epsilon e_n)$ and the supremum on $C_-$ to $u(x-\epsilon e_n)$ by continuity, where $e_n=(0,\dots,0,1)$. Since $u_n>0$, we have $u(x+\epsilon e_n)>u(x)>u(x-\epsilon e_n)$.
Take any such $\delta>0$. Denote $a_+:=\inf_{C_+}u$ and $a_-:=\sup_{C_-}u$. Let $U=\{x\in T;a_-<u(x)<a_+\}$. This $U$ is now a "convenient neighborhood of $x$". The functions given by the implicit function theorem are1 all defined over $\tilde B:=B^{n-1}(x',\delta)$. The mental image is that $U$ is foliated by the level sets of $u$, and they all have the same projection $\tilde B$.
In the definition of $U$ we cut the tube $T$ so that all level sets are "full" (their projections are the whole $\tilde B$). But if $\delta$ is too large, this cut might make $U$ empty. It would be unfortunate if $x\notin U$.
Every point $x\in\Omega$ has a neighborhood like this, and the rest of the proof works from here.
1 Elaboration: As you wrote in your question, it follows from the implicit function theorem that the level sets of $u$ can be written as graphs in the set $B=B^n(x,\epsilon)$. The problem was that the functions (whose graphs we have) were defined over different sets in $\mathbb R^{n-1}$. For any $x\in U$, consider the level set $L=U\cap u^{-1}(u(x))$. We want to show that the projection of $U$ is $\tilde B$. Now $a_-<u(x)<a_+$. Take any $y'\in\tilde B$ and consider the function $h(t)=u(y',t)$. The function $h$ is strictly increasing. The point $(y',t)$ is in $\bar T$ when $|y'-x'|^2+(t-x_n)^2\leq\epsilon^2$ or $t\in[x_n-\sqrt{\epsilon^2-|y'-x'|^2},x_n+\sqrt{\epsilon^2-|y'-x'|^2}]=:[b_-,b_+]$. By definition $h(b_+)\geq a_+>a_-\geq h(b_-)$. Shrinking $T$ to $U$ corresponds to shrinking $[b_-,b_+]$ to $[c_-,c_+]$, where $h(c_\pm)=a_\pm$. Since $h|_{[c_-,c_+]}$ is strictly increasing and continuous, there is exactly one $t\in(c_-,c_+)$ so that $h(t)=u(x)$ (recall that $a_-<u(x)<a_+$). This means that the point $(y',t)\in U$ projects down to $y'$ and $u$ has the value $u(x)$ on it.