OP's actual question follows directly from Theorem 5.3 in Ref. 2, but that leaves the obvious question: How to prove Theorem 5.3? It seems the only really satisfying answer would be to outline a complete proof of the Liouville-Arnold Theorem. This is what we intend to do in this answer.
On one hand, from the commuting Hamiltonian vector fields $X_k:=-\{F_k, \cdot\}$, we generate an abelian flow $$\begin{align} g(t)~:=~& \exp(\sum_{k=1}^n t^k X_k), \cr
g(t)g(t^{\prime})~=~&g(t+t^{\prime}), \qquad
t,t^{\prime}~\in~\mathbb{R}^n .\end{align} \tag{1}$$
We focus on a connected component of a compact level-set ${\cal M}_f$. We conclude from the compactness assumption that the $t$-parameters generate a periodic torus action $t^k\sim t^k+\Pi^k_{\ell}(F)$ (for each connected component), where $\Pi^k_{\ell}(F)$ is a period matrix.
On the other hand, from Frobenius theorem for the involutive distribution $\Delta:={\rm span}_{\mathbb{R}}(X_1,\ldots,X_n)$, there exists an $n$-foliation. In other words, there exists an atlas of local coordinates of the form $(\varphi^1, \ldots,\varphi^1, G_1, \ldots, G_n)$, where firstly the $\varphi^k$-coordinates parametrize the leaves, and secondly the $G_k$-coordinates (like the $F_{\ell}$-coordinates) label the leaves. Hence the $G_k$-coordinates must be functions of $F_{\ell}$-coordinates. Two overlapping local coordinates systems (say, primed and unprimed) are related via
$$\begin{align} \varphi^{\prime k} ~=~& \varphi^{\prime k}(\varphi,G), \cr
G^{\prime}_k ~=~&G^{\prime}_k(G), \cr
k~\in~&\{1, \ldots, n\}.\end{align}\tag{2}$$
We may w.l.o.g. assume that the local $\varphi^k$-coordinates always stratify the Hamiltonian vector fields
$$ \frac{\partial}{\partial \varphi^k}~=~X_k, \qquad k~\in~\{1, \ldots, n\}.\tag{3}$$
Notice that the abelian flow then becomes
$$g(t)~\stackrel{(1)+(3)}{=}~ \exp(\sum_{k=1}^n t^k \frac{\partial}{\partial \varphi^k}) , \qquad t~\in~\mathbb{R}^n .\tag{4}$$
It becomes clear that we then should identify the parameter $t^k$ with the coordinate $\varphi^k$. Eqs. (2) & (3) narrow coordinate transformations down to the form $\varphi^{\prime k} = \varphi^k+ h^k(F)$. It becomes clear that we can take unions of such local coordinate patches to create a $\varphi$-complete coordinate system (for each connected component) with a periodicity condition $\varphi^k\sim \varphi^k+\Pi^k_{\ell}(F)$. In other words, the $\varphi^k$ are un-normalized angle variables.
Given a connected component of a compact level set, which must be an $n$-torus $\cong (\mathbb{S}^1)^n$, we can find a tubular neighborhood ${\cal U}\subseteq {\cal M}$ with coordinates $(\varphi^1,\ldots,\varphi^n,F_1, \ldots, F_n)$, where the $F$-coordinates live in some contractible space, say, an $n$-box. We now restrict the construction to this tubular neighborhood ${\cal U}$.
The fundamental Poisson brackets are not necessarily on Darboux form$^1$
$$\begin{align} \{ \varphi^k, \varphi^{\ell}\} ~=~& \beta^{k\ell}, \cr
\{ \varphi^k, F_{\ell}\}~=~& \delta^k_{\ell}, \cr
\{ F_k, F_{\ell}\}~=~& 0, \cr
k,\ell ~\in~&\{1, \ldots, n\}.\end{align}\tag{5}$$
The corresponding symplectic form is the inverse structure
$$\begin{align}\omega~=~&\sum_{k=1}^n\mathrm{d}F_k \wedge \mathrm{d}\varphi^k + \beta, \cr
\beta~:=~& \frac{1}{2} \sum_{k,\ell=1}^n \mathrm{d}F_k ~\beta^{k\ell} \wedge \mathrm{d}F_{\ell}.\end{align} \tag{6}$$
However, there is a bi-grading of two anti-commuting differentials
$$\begin{align}\mathrm{d}~=~&\mathrm{d}^{(\varphi)}+ \mathrm{d}^{(F)}, \cr
\mathrm{d}^{(\varphi)}~:=~&\sum_{k=1}^n\mathrm{d}\varphi^k \frac{\partial}{\partial \varphi^k},\cr
\mathrm{d}^{(F)}~:=~&\sum_{k=1}^n\mathrm{d}F_k \frac{\partial}{\partial F_k} .\end{align} \tag{7}$$
The closedness $\mathrm{d}\beta=0$ implies firstly $\mathrm{d}^{(\varphi)}\beta=0$
(i.e. $\beta^{k\ell}$ can not depend on $\varphi$.), and secondly $\mathrm{d}^{(F)}\beta=0$. By Poincare Lemma there exists a $(0,1)$-form
$$\alpha~=~\sum_{k=1}^n\alpha^k(F)~\mathrm{d}F_k\tag{8}$$
on ${\cal U}$ such that the symplectic $(0,2)$-form
$$\beta|_{\cal U}~=~\mathrm{d}^{(F)}\alpha~=~\mathrm{d}\alpha.\tag{9}$$
Define new coordinates
$$\begin{align} q^k~:=~\varphi^k -\alpha^k(F)~\sim~& q^k+\Pi^k_{\ell}(F),
\cr k,\ell ~\in~&\{1, \ldots, n\},\end{align}\tag{10} $$
which are periodic with the same periods. It is easy to check that
$(q^1,\ldots,q^n,F_1, \ldots, F_n)$ is a $q$-complete Darboux coordinate neighborhood ${\cal U}$ with a symplectic potential $(1,0)$-form
$$\begin{align}\vartheta~=~&\sum_{k=1}^nF_k ~\mathrm{d}q^k, \cr
\mathrm{d}\vartheta~=~&\sum_{k=1}^n\mathrm{d}F_k \wedge \mathrm{d}q^k
~=~\omega|_{\cal U}.\end{align}\tag{11}$$
At this point we take a step backwards. (This is not needed for the AA existence proof, but in practice it is useful to outline the construction for more general coordinates. See also remark 7 below.) Let us more generally imagine that we have a tubular neighborhood ${\cal U}\subseteq {\cal M}$ with some $q$-complete (not necessarily Darboux) coordinates $(q^1,\ldots,q^n,F_1, \ldots, F_n)$
with a symplectic potential $(1,0)$-form
$$\vartheta~=~\sum_{k=1}^n p_k(q,F) ~\mathrm{d}q^k,\tag{12}$$
such that $q^k\sim q^k+\Pi^k_{\ell}(F)$.
Consistency with the non-degenerate symplectic 2-form
$$\begin{align}\omega|_{\cal U}~=~& \mathrm{d}\vartheta\cr
~\stackrel{(12)}{=}~& \underbrace{\sum_{k,\ell=1}^n\mathrm{d}q^k ~\frac{\partial p_{\ell}}{\partial q^k} \wedge \mathrm{d}q^{\ell}}_{=0\text{ because }\{ F_k, F_{\ell}\}=0 }
+\sum_{k,\ell=1}^n\mathrm{d}F_k ~\frac{\partial p_{\ell}}{\partial F_k} \wedge \mathrm{d}q^{\ell},\end{align}\tag{13}$$
implies:
that $\{q^k,q^{\ell}\}=0$.
the Maxwell relations
$$ \frac{\partial p_{\ell}(q,F)}{\partial q^k}
~=~ (k\leftrightarrow \ell) , \qquad k,\ell~\in~\{1, \ldots, n\}.\tag{14}$$
Eq. (14) is an integrability condition for the local existence of the Hamilton's characteristic function $W(q,F)$, cf. Section 6 below.
that the matrix $\frac{\partial p_{\ell}}{\partial F_k}$ is invertible. After possibly restricting to a smaller tubular neighborhood ${\cal V}\subseteq {\cal U}$, this proves that $(q^1,\ldots,q^n,p_1, \ldots, p_n)$ are Darboux coordinates on ${\cal V}$. (Here we have used the inverse function theorem.)
Define action variables
$$J_k(q_{\ast},F)~:=~ \oint_{C_k(q_{\ast})}\! \left. \vartheta \right|_{\text{fixed } F}, \qquad k~\in~\{1, \ldots, n\}, \tag{15}$$
where $C_k(q_{\ast})$ are a $k$'th 1-cycle of the tori in the $q$-space $M$ (aka. configuration space) that starts and ends at the same point $q_{\ast}\in M$.
Use Stokes' theorem
$$\begin{align} J_k(q^{\prime}_{\ast},F)-J_k(q_{\ast},F)~:=~& \oint_{\partial A}\! \left.\vartheta \right|_{\text{fixed } F}\cr
~=~& \iint_A\! \left. \omega \right|_{\text{fixed } F}~=~0, \cr
k~\in~&\{1, \ldots, n\}, \end{align}\tag{16}$$
to show that $J_k(q_{\ast},F)$ do not depend on $q_{\ast}$.
At this point we assume that the matrix $\frac{\partial J_k}{\partial F_{\ell}}$ is non-degenerate. (This is satisfied for the important case where $p_k(q,F)=F_k$.) After possibly restricting to a smaller tubular neighborhood ${\cal V}\subseteq {\cal U}$, this proves that $(q^1,\ldots,q^n,J_1, \ldots, J_n)$ are coordinates on ${\cal V}$. (Here we have used the inverse function theorem.)
Define Hamilton's characteristic function
$$W(\gamma,F)~:=~ \int_{\gamma}\!\left. \vartheta \right|_{\text{fixed } F}, \tag{17}$$
where $\gamma \subseteq M$ is an oriented curve in the $q$-space $M$ from a fiducial point $q_{\ast}(F)$ to an endpoint, which we (with a slight abuse of notation) denote $q$. One may again show that the definition (17) only depend on the homotopy class (i.e. winding numbers) of $\gamma$. Hence we may rewrite definition (17) in a slight different notation as
$$W(q,F)~:=~ \int_{q_{\ast}}^q\!\left. \vartheta \right|_{\text{fixed } F}, \tag{18}$$
where we use the multi-valued nature of the $q$-coordinate system to indicate the winding number. Note that
$$ p_k(q,F)~=~\frac{\partial W(q,F)}{\partial q^k}, \qquad k~\in~\{1, \ldots, n\}. \tag{19}$$
Define $W(q,J)~:=~W(q,F(J))$.
Similarly,
$$ p_k(q,J)~=~\frac{\partial W(q,J)}{\partial q^k}, \qquad k~\in~\{1, \ldots, n\}. \tag{20}$$
Then if we shift with a period
$$\begin{align} \Delta_k W ~:=~&W(q+\Pi_k,J)-W(q,J)~=~J_k, \cr k~\in~&\{1, \ldots, n\}.\end{align} \tag{21}$$
Define new angle variables
$$ w^k(q,J)~:=~\frac{\partial W(q,J)}{\partial J_k}, \qquad k~\in~\{1, \ldots, n\}. \tag{22}$$
In other words, $W$ is a type-2 generating function for the canonical transformation $(q,p)\to (w,J)$. Then the periods have unit-length
$$\begin{align} \Delta_{\ell} w^k~:=~& w^k(q+\Pi_{\ell},J)-w^k(q,J)~=~ \delta^k_{\ell}, \cr k,\ell~\in~&\{1, \ldots, n\}. \end{align}\tag{23}$$
$\Box$
Remark. If a tubular neighborhood ${\cal U}\subseteq {\cal M}$ of a level set torus is covered by a single Darboux coordinate chart $(q^1,\ldots,q^n,p_1,\ldots,p_n)$ such that the projection of a level set torus to the $q$-space $M$ only contain isolated turning points, then it is still possible to define action variables (15), Hamilton's characteristic function (17), etc, via appropriate generalizations. E.g. the analogue of eq. (17) is
$$W(\gamma,F)~:=~ \int_{\sigma_F\circ \gamma}\!\left. \vartheta \right|_{\text{fixed } F}, \tag{24}$$
where $\sigma_F$ denotes the $F$-dependent lift from $q$-space $M$ to the level set torus.
Example: The simple harmonic oscillator. The Hamiltonian
$$H~=~\frac{p^2}{2} + \frac{q^2}{2}\tag{25} $$
is an integral of motion. The system is integrable in the 2D phase space (except for the origin, which is a singular orbit). We can define an angle variable
$$\begin{align} \varphi~:=~&{\rm atan2}(q,p), \cr p+iq~=~&\sqrt{2H}e^{i\varphi}, \cr
q~=~&\sqrt{2H}\sin\varphi, \cr
p~=~&\sqrt{2H}\cos\varphi\cr
~=~&\pm \sqrt{2H-q^2}.\end{align}\tag{26} $$
The action variable reads
$$\begin{align} J~\stackrel{(15)}{=}~&\oint \!\left. p ~\mathrm{d}q \right|_{\text{fixed } H}\cr
~=~&2H\int_0^{2\pi}\! \mathrm{d}\varphi~\cos^2\varphi\cr
~=~&2\pi H.\end{align}\tag{27}
$$
Hamilton's characteristic function reads
$$\begin{align}W(q,H)~\stackrel{(18)}{=}~&\int_{q_{\ast}}^q \!\left. p ~\mathrm{d}q \right|_{\text{fixed } H}\cr
~=~&2H\int_0^{\varphi(q,H)}\! \mathrm{d}\varphi~\cos^2\varphi\cr
~=~&H\left. \{\varphi + \frac{1}{2}\sin 2 \varphi \}\right|_{\varphi=\varphi(q,H)}\cr
~=~&\left. \{ H~{\rm atan2}(q,p) + \frac{1}{2} q p\}\right|_{p=p(q,H)}.\end{align}\tag{28} $$
The angle variable becomes
$$\begin{align}2\pi w~\stackrel{(22)}{=}~& 2\pi\frac{\partial W(q,J)}{\partial J}\cr
~\stackrel{(27)}{=}~&\frac{\partial W(q,H)}{\partial H}\cr
~\stackrel{(28)}{=}~& {\rm atan2}(q,p) + \left\{H\frac{1}{1+(q/p)^2} \frac{-q}{p^2}+\frac{q}{2}\right\}\frac{\partial p(q,H)}{\partial H}\cr
~=~&{\rm atan2}(q,p)\cr
~=~&\varphi,\end{align} \tag{29}$$
as expected.
