$\def\RR{\mathbb{R}}\def\CC{\mathbb{C}}\def\NN{\mathbb{N}}$I'm going to try once more. As my previous flawed answers should have made obvious, the real analytic category is not my home, so read with caution.
Conveniently, the two big results I need are stated as Corollary 5.43 and 5.44 here:
Cor. 5.43 $M$ has a proper real analytic embedding into some $\RR^n$. I'll write $(z_1, \ldots, z_n)$ for the coordinates on $\RR^n$.
Cor. 5.44 (specialized and rephrased) Let $p$ be a point of $M$ and let $M$ be embedded in $\RR^n$ as above. Let $U$ be a neighborhood of $p$ in $\RR^n$. Let $f: M \cup U \to \RR$ be a function which is real analytic on $U$ and on $M$. Let $d$ be a nonnegative integer. Then there is a real analytic function $g: \RR^n \to \RR$ such that $g|_M=f$ and $(g-f)|_U$ vanishes to order $d$ at $p$.
We also need the following:
Easy Lemma: Let $f: \RR^n \to \RR$ be a real analytic function for which $f(z_1, \ldots, z_{n-1},0)$ is identically zero. Then $f(z_1, z_2, \ldots, z_{n-1}, z_n)/z_n$ extends to a real analytic function $\RR^n \to \RR$.
Proof This is a local statement, and can be easily checked locally using convergent power series. $\square$
Key Lemma: Let $g: \RR^n \to \RR$ be a real analytic function vanishing to order $d$ at $0$. Then we can write $g(z)$ as a finite sum $\sum z^a p_a(z)$ where $z^a$ are monomials of degree $d$ and $p_a(z)$ are real analytic functions.
Proof
By induction on $d$ and $n$. The base case $d=0$ is trivial: Write $g(z) = 1 \cdot g(z)$; the base case $n=0$ is vacuous. We assume $d$ and $n>0$. Define $g(z_1, \ldots, z_n)$ by
$$f(z_1, \ldots, z_n) = f(z_1, \ldots, z_{n-1}, 0) + z_n g(z_1, \ldots, z_{n-1}).$$
By the Easy Lemma, $g$ is analytic and, by a basic computation, it vanishes to order $d-1$. By induction on $d$, we can write $g = \sum z^c r_c(z)$ where the monomials $z^c$ have degree $d-1$. By induction on the number of variables, we can write $f(z_1, \ldots, z_{n-1},0)$ as $\sum z^b p_b(z_1, \ldots, z_{n-1})$ where the monomial $z^b$ have degree $d$. So
$$f(z_1, \ldots, z_n) = \sum z^b p_b(z_1, \ldots, z_{n-1}) + \sum (z^c z_n) q_c(z_1,\ldots, z_n),$$
an expression of the desired form. $\square$
We now prove the result.
Lemma Let $f$ be a real analytic function on $M$ that vanishes to second order at $p$. Then $(D f)(p)=0$.
Proof Embed $M$ in $\RR^n$ by Cor 5.43. Extend $f$ to a function on $\RR^n$ also vanishing to order $2$ at $p$ and, without loss of generality, translate $p$ to $0$. By the Key Lemma, we can write $f(z) = \sum z_i z_j c_{ij}(z)$. Then
$$(D f)(0) = \sum (D z_i)(0) \cdot 0 \cdot c_{ij}(0) + \sum 0 \cdot (D z_j)(0) \cdot c_{ij}(0) + \sum 0 \cdot 0 \cdot (D c_{ij})(0) =0. \quad \square$$
Lemma $(D f)(p)$ depends only on $(df)(p)$.
Proof Suppose that $(d f_1)(p) = (d f_2)(p)$. Define $q(z)$ by the equation $f_2(z) = f_1(z) + (f_2(p)-f_1(p)) + q(z)$, then $q$ vanishes to order $2$ at $p$ so $(D q)(p) =0$. Also, $D$ of a constant function is $0$. So $(D f_1)(p) = (D f_2)(p)$. $\square$
The previous lemma shows that there is a well defined map $v(p) : T^{\ast}_p M \to \RR$
so that $(D f)(p) = v(p)(d f)$. (Well, we also need to see that $C^{\omega}(M) \to T^{\ast}_p M$ is surjective, but that is obvious because the functions $z_i$ span the cotangent space.) It is easy to see that $v(p)$ is linear. So $v(p)$ is an element of $(T_{\ast})_p M$. We have now created a section $v: M \to T_{\ast} M$ so that $(D f)(p) = \langle v(p), df \rangle$. It remains to see that $v$ is real analytic.
This question is local. Let $(z_1, \ldots, z_d)$ give local coordinates on $M$ near some point $p$. Then we can write
$$v(p) = \sum_{i=1}^d v_i(p) \frac{\partial}{\partial z_i}.$$
We need to show that the functions $v_i$ are real analytic. Since $v_i = D z_i$, this is clear. QED.
Best Answer
For real-analytic manifolds, coherence of the structure sheaf always holds. The 1-sentence reason is that one can pass to real and imaginary parts on Oka's coherence theorem in several complex variables.
Before discussing a proper proof of that 1-sentence executive summary, I should address that in the real-analytic setting however an analogue of one of Oka's "coherence" results over $\mathbf{C}$ can fail (leading one to see phrases in the literature such as "non-coherent real-analytic spaces"): there is no real counterpart to the "analytic Nullstellensatz" of Oka. To explain this, recall Oka's result that for a complex-analytic set $X$ in a complex manifold $V$ (i.e., a closed subset that is locally on $V$ given by the vanishing of finitely many holomorphic functions) the ideal sheaf of sections of $O_V$ vanishing on $X$ is always locally finitely generated (or equivalently a coherent $O_V$-module, since $O_V$ is coherent by Oka's big theorem). But this fails in the real-analytic case. In contrast, there exist real-analytic sets $Z$ in real-analytic manifolds $U$ (i.e., a closed subset that is locally on $U$ given by the vanishing of finitely many real-analytic functions) such that the ideal sheaf $I_Z$ of sections of $O_U$ vanishing on $Z$ is not locally finitely generated; such a $Z$ is called "non-coherent" in $U$ (because in such cases the subsheaf $A = O_U/I_Z$ of the sheaf of $\mathbf{R}$-valued continuous functions on $Z$ can fail to be a coherent sheaf of rings, though I have never cared enough to check up on a proof of this consequence in some cases). An explicit example of $Z$ and $U$ with $I_Z$ not coherent inside $O_U$ is given near the start of https://arxiv.org/pdf/math/0612829.pdf.
To be more detailed about the proof of coherence of the structure sheaf of a real-analytic manifold, the assertion is of local nature and so only depends on the (local) dimension; i.e., for manifolds of (pure) dimension $n$ it is equivalent to showing that the sheaf $O_U$ of real-analytic functions on every small open ball $U$ around the origin in $\mathbf{R}^n$ is coherent. By definition, this coherence amounts to local finite generation for the kernel of any $O_U$-linear map $\varphi:O_U^{\oplus N} \rightarrow O_U$ for any small $U$. This map locally "extends" over an open in $\mathbf{C}^n$; i.e., by working locally on $U$ we can arrange that that exists an open $V \subset \mathbf{C}^n$ satisfying $V \cap \mathbf{R}^n = U$ and holomorphic $F_1, \dots, F_N$ on $V$ whose restriction to $U$ recovers the $N$ components $f_1, \dots, f_N$ of $\varphi$.
By Oka's big coherence theorem (not the analytic Nullstellensatz above), the resulting map $\Phi: O_V^{\oplus N} \rightarrow O_V$ extending $\varphi$ has kernel that is locally finitely generated. Working locally on $V$ around points of $U$, we can thereby arrange that there exist $s_1, \dots, s_r \in O(V)^{\oplus N}$ generating $\ker \Phi$. We claim that the real and imaginary parts of the restrictions to $U$ of these $N$-tuples belong to $\ker \varphi$ and generate it.
Our problem is local on $V$ near $U$, so for generation we can focus on $(g_1,\dots,g_N) \in (\ker \varphi)(U)$. By working locally on $V$ around a point in $U$ we can arrange that $V$ is connected and that each $g_j$ extends (necessarily uniquely) to a holomorphic $G_j$ on $V$. Then the global section $\sum G_j F_j \in O(V)$ has restriction to $U$ equal to $\sum g_j f_j = \varphi(g_1,\dots,g_N)=0$, so $\Phi(G_1,\dots,G_N)=\sum G_j F_j = 0$ by connectedness of $V$ and analyticity considerations. Thus, working locally on $V$ some more we can arrange that $(G_1,\dots,G_N) = \sum_{k=1}^r a_k s_k$ for some $a_{1}, \dots, a_{r} \in O(V)$, so restricting to $U$ gives $$(g_1,\dots,g_N) = \sum_{k=1}^r a_{k}|_U \cdot s_k|_U.$$
Each function $g_j$ is $\mathbf{R}$-valued, whereas $a_{k}|_U$ and the $N$ components of $s_k|_U$ are $\mathbf{C}$-valued, and the real and imaginary parts of these various $\mathbf{C}$-valued functions are real-analytic on $U$. Thus, $$(g_1,\dots,g_N) = \sum_{k=1}^r ({\rm{Re}}(a_{k}|_U){\rm{Re}}(s_k|_U) - {\rm{Im}}(a_{jk}|_U) {\rm{Im}}(s_k|_U)).$$ Finally, the condition $s_k \in (\ker \Phi)(V) \subset O(V)^{\oplus N}$ says that the $N$ components $s_{k1},\dots,s_{kN} \in O(V)$ satisfy $\sum F_j s_{kj} = 0$, so $\sum f_j \cdot s_{kj}|_U = 0$ as $\mathbf{C}$-valued functions on $U$. But each $f_j$ is $\mathbf{R}$-valued, so $\sum f_j {\rm{Re}}(s_{kj}|_U)$ and $\sum f_j {\rm{Im}}(s_{kj}|_U)$ both vanish on $U$; i.e., ${\rm{Re}}(s_k|_U)$ and ${\rm{Im}}(s_k|_U)$ belong to $(\ker \varphi)(U)$. Thus, ${\rm{Re}}(s_1|_U), \dots, {\rm{Re}}(s_r|_U), {\rm{Im}}(s_1|_U), \dots, {\rm{Im}}(s_r|_U)$ generate $\ker \varphi$ over the (now shrunken) $U$. This shows that in the initial setup (before we began shrinking $U$) the kernel of $\varphi$ is locally finitely generated.