If you call the standard $n$-simplex $\Delta^n$, i.e.
$$\Delta^n = \{ (x_0, x_1, \cdots, x_n) : x_i \geq 0 \forall i, \sum_i x_i = 1\}$$
then the function
$\phi : \Delta^n \to \mathbb R$
given by $\phi(x_0, \cdots, x_n) = \sum_i x_i^2$
almost has the properties we want. It has critical values $1, 1/2, 1/3, \cdots, 1/(n+1)$ located on the barycentres of the various facets of the simplex. $1$ is on the vertices, and $1/(n+1)$ is at the barycentre of $\Delta^n$, i.e. at $x_i = 1/(n+1) \forall i$.
The problem with this function if you use the characteristic maps of the simplicies of a smooth triangulation of $M$ to lift this function to a function $f : M \to \mathbb R$, it will only be piecewise smooth, with lack of smoothness on the $(m-1)$-skeleton -- specifically the derivatives normal to the $(m-1)$-dimensional simplices generally won't exist. I see a few ways to address this, my preference being to apply a bump-function construction in a regular neighbourhood of the $(m-1)$-skeleton to dampen the normal derivatives of $f$, so that it becomes $C^2$-smooth on $M$
One way to accomplish this (although its not the prettiest formula) is to take the homeomorphism of $\Delta^n$ given by a formula where $\vec x = (x_0, \cdots, x_n)$
$$\vec x \longmapsto \alpha \vec x + \beta \frac{(x_0, \cdots, \hat x_j, \cdots x_n)}{1-x_j}$$
Where $\alpha + \beta = 1$, and we are assuming $x_j = \min\{x_0, \cdots, x_n\}$, and $\hat x_j = 0$. $\alpha$ is only a function of $x_j$, say it is an approximation to $\alpha(x_j) = 2x_j$, with domain $[0,1/2]$ but we will demand the 1st derivative to $\alpha$ at $0$ is zero, and $\alpha(1/2)=1$, and $\alpha'(1/2)=0$, i.e. otherwise $\alpha$ is increasing.
If you think carefully about this it won't be quite $C^1$ as there are issues on the co-dimension two subspaces where $x_i=x_j$ for $i \neq j$. But that can be fixed by essentially by doing this construction for all the boundary facets, i.e. not just the co-dimension one facets. You have to be careful not to destroy the Morse properties of the function, which requires careful choice of $\alpha$.
To be more precise, if both $x_j$ and $x_k < 1/2$ then
$$\vec x \longmapsto \alpha \vec x + \beta \frac{(x_0, \cdots, \hat x_j, \cdots, \hat x_k, \cdots x_n)}{1-x_j-x_k}.$$
i.e. I'm suggesting doing a simultaneous damping on the nearness to the various strata -- this avoids ruining the Morse nature of the function, higher-order critical points, etc.
This isn't a complete answer to your question but I'll put it up here for now and will revise it when/if I have a more transparent response. I suspect there's a really clean analytic formula for the function, but the inspiration hasn't hit. In the mean time, perhaps someone else will see a more efficient technique.