Can we specify the value of harmonic forms at a point? Let $M$ be a smooth $d$-dimensional oriented Riemannian manifold, and let $1 < k < d$ be fixed. 
Let $p \in M$, and let $\alpha_p \in \bigwedge^k(T_pM)^*$.

Does there exist an open neighbourhood $U$ of $p$, and a closed and co-closed $k$-form $\omega \in \Omega^k(U)$ satisfying $\omega_p=\alpha_p$?


This question is equivalent to the following question:

Do closed and co-closed frames for $\bigwedge^k(T^*M)$ always exist locally?

Indeed, if we can specify the value of a form in a point, we can take a basis for  $\bigwedge^k(T_pM)^*$, and so obtain forms which form a frame at $p$. Since "being a frame" is an open condition, we have a local frame. On the other hand, suppose that local closed and co-closed frames exist. Then, by choosing a linear combination with constant coefficients of that frame, we can realize any given value at $p$. 
Comment: In general we cannot expect such a frame to be induced from coordinates. Indeed, when we specialize to even dimension $d$, and  $k=\frac{d}{2}$, then, for a generic metric $g$, there are no coordinate systems where even one wedge $\mathrm{d}x^{i_1}\wedge\cdots\wedge\mathrm{d}x^{i_n}
$ is co-closed. 
 A: The answer to your question is 'yes', but it may take me a little while to look up the appropriate references, since I'm traveling now.  
The main point is that the (usually overdetermined) system of PDE defined by $\mathrm{d}\alpha = \mathrm{d}^*\alpha = 0$ for $\alpha\in\Omega^k(M)$ is involutive and hence, if the metric is real-analytic, it's formally integrable, so one can specify the $0$-jet of a local solution arbitrarily.  Meanwhile, because the elliptic complex $\mathrm{d}+\mathrm{d}^*:\Omega^*(M)\to\Omega^*(M)$ satisfies Spencer's $\delta$-estimate, this local formal integrability carries over to the smooth case, which is what you are asking.
I'll get back home later in the week, and I'll be able to post appropritate references then.
A: Take a torus $\mathbb{T}^d$ of dimension $d$ and introduce on it such a metric that some part of it would be isometric to a (sufficiently small) neighbourhood  $U$ of $p\in M$. By Hodge theory, harmonic $k$ - forms on it are in one-to-one correspondence to cohomology classes in $H^k(\mathbb{T}^d,\mathbb{R})\cong \bigwedge^k H^1(\mathbb{T}^d,\mathbb{R})$. Then each element of $\bigwedge^k T_p^* M$ is represented by some harmonic form simply by dimensional consideration.
(To avoid problems, we may assume that the metric is flat everywhere except a small patch. Then these forms won't vanish at $p$.)
[EDIT] All right, I will add some details. 


*

*Choose a chart $\phi: U\to U'$ where $U$ is a neighbourhood of $p$ and 
$U'\subset \mathbb{R}^d$ is a domain in the euclidean space with the usual metric 
$g_0$,
$$ds^2=\sum_i (dx^i)^2.$$
We assume that $\phi(p)=0$ and that the induced metric $\phi^*g_0$ coincides with 
the original metric at $p$.  

*Introduce a new metric $g$ on $\mathbb{R}^d$ with the following properties. 
(a) $g=g_0$ outside a neighbourhood of the origin. (b) $|g-g_0|<\epsilon$ everywhere. (c) The chart $\phi$ is an isometry near $p$. (d) Second derivatives  $\frac{\partial^2 g}{\partial x_i\partial x_j}$ are uniformly bounded. 
(It is not difficult to achieve all this using a partition of unity.)

*Changing the metric $g$ to a periodic one in the obvious way, we may interpret 
$U'$ as a domain in a Riemannian manifold $M'$ with this metric which is homeomorphic to a torus $\mathbb{T}^d$. 
It is well known that the space of harmonic k-forms on $M'$  is isomorphic to $H^k(M',\mathbb{R})\cong \bigwedge^k H^1(\mathbb{T}^d,\mathbb{R})$, so its dimension is 
$\binom{d}{k}$. Moreover, in the outer domain where $g=g_0$ these forms are linear combinations of $\omega_{i_1,\dots,i_k}=dx^{i_1}\wedge\dots\wedge dx^{i_k}$. Note that if $\omega$ is a harmonic form on $M'$ then $\phi^*\omega$ is a harmonic form on $M$ near $p$. 
What is left is to make sure  that no harmonic forms on $M'$ vanish at the origin. (While the forms $\omega_{i_1,\dots,i_k}$ may be extended to  the whole $M'$,  they do not have the above nice form near the origin.)  This part is a little technical, but it certainly can be done. It is possible to expand the 
"good" domain (where $g=g_0$) as much as we like while keeping $\frac{\partial^2 g}{\partial x_i\partial x_j}$ bounded. The coefficients at $dx^{i_1}\wedge\dots\wedge dx^{i_k}$ of a harmonic form cannot behave too wildly (due to the Weitzenböck identity), in particular, they cannot suddenly vanish.
P.S. I think the idea is simple and transparent, but the proof becomes somewhat messy once you fill all the details in.  A more conceptual approach suggested by Bryant may be better after all.
A: $\newcommand{\M}{\mathcal{M}}$
$\newcommand{\ep}{\epsilon}$
I will sketch here a different approach than the one given by Robert's answer. (It is loosely based on Alex's answer). We want to prove that around every point $p \in \M$ there exist a local frame for $\bigwedge^k(T^*\M)$ whose elements are closed and co-closed. 

For the Euclidean metric this is immediate: We have the standard (constant) frame $dx^I=dx^{i_1} \wedge \ldots dx^{i_k} $. Since every metric is locally close to being Euclidean on small neighbourhoods, the idea is to use an approximation argument:

Given a Riemannian metric $g$, we denote the space of $g$-harmonic forms of degree $k$ by $H^k_{g}$.
We view $H^k_{g}$, as a subspace of $\Omega^k(\M)$ which is "changing continuously" with the metric $g$. Suppose $g_{\ep} \to g_0$ in the $C^1$-norm where $g_0$ is the Euclidean metric; Then $H^k_{g_{\ep}} \to H^k_{g_0}$ in the following sense: there exist a family of bases of $H^k_{g_{\ep}}$, which converges to a basis of $H^k_{g_{0}}  $ in $C^1$; this basis of $H^k_{g_{0}}$ forms a local frame for $\bigwedge^k(T^*\M)$. Since being a frame is an open condition, those bases for $H^k_{g_{\ep}}$ are local frames for sufficiently small $\ep$. 
For the full details, see Appendix A in my paper here.
Some more details:
Even though the claim is local, and the approximation scheme is also inspired by a local phenomena, the implementation of the proof is based on a combination of local and global arguments. The reason is that on a closed manifold, being closed and co-closed is equivalent to being harmonic, and the dimension of the space of harmonic forms is a finite number which is a topological invariant of the manifold; it does not depend on the chosen metric. 
Thus, given a family of metrics $g_{\ep} \to g_0$ on a closed manifold $\M$, we consider the behaviour of the finite-dimensional subspaces $H^k_{g_{\ep}}$ (all of the same dimension) as $\ep \to 0$.
That is, we look at the map $g \to H^k_{g}=\ker \Delta_g$. It turns out that this map is continuous in some appropriate sense; this relies on a certain "stability property of kernels of linear operators". It turns out that a crucial factor in the existence of such a stability phenomenon is the assumption that all the kernels have the same finite dimension. The convergence of kernels does not always hold when the dimensions are not equal or infinite.
