The answer is 'yes, a smooth, flat solution $f$ exists when $g$ is smooth and flat'.

Here is one way to show this:  I'll first do the case in which $g$ is *even*, i.e., $g(-p,-q)=g(p,q)$ and, for convenience, I'll assume that $g$ is defined on the entire $pq$-plane. (See the remark at the end about the local case.)  

Let $(u,v) = (p^2{-}q^2,\,2pq)$, and note that there is a (unique) function $\bar g$ on the $uv$-plane such that $g(p,q) = \bar g(u,v)$ and that $\bar g$ is smooth and flat at $(u,v)=(0,0)$.  I will look for a solution of the form $f(p,q) = \bar f(u,v)$.  By the Chain Rule, the equation $(\star)$ then becomes
$$
2\sqrt{u^2{+}v^2} \, \frac{\partial \bar f}{\partial u} = \bar g,
$$
so 
$$
\frac{\partial \bar f}{\partial u}(u,v) = \frac{\bar g(u,v)}{2\sqrt{u^2{+}v^2}}.
$$
The right hand side is smooth and flat at $(u,v) = (0,0)$, so one has a solution in the form
$$
\bar f(u,v) = \int_0^u \frac{\bar g(t,v)}{2\sqrt{t^2{+}v^2}}\,\mathrm{d}t\,.
$$
This $\bar f$ is flat and gives a solution to the problem, in fact, the unique solution that satisfies $\bar f(0,v) = 0$.

In the general case, one can write $g = g_0 + g_1$, where $g_0$ is even and $g_1$ is *odd*, i.e., $g_1(-p,-q) = -g_1(p,q)$.  So, to finish, by the linearity of the equation, it only remains to solve the equation when $g$ is odd.  This can be done by using the above solution on the half-planes $p>0$ and $p<0$ and being a little careful about the matching.  However, the right way to think about it in the odd case is that $\bar g$ is actually a section of a nontrivial flat line bundle over the punctured $uv$-plane, and the above integral is then taken using parallel translation in the flat line bundle along segements of the form $\sigma(t) = (tu,v)$ for $0\le t\le 1$, thereby yielding a section $\bar f$ of this nontrivial line bundle over the punctured $uv$-plane.  This defines a solution $f(p,q)$ of the equation on the punctured $pq$-plane that is odd and that vanishes to infinite order at $(p,q)=(0,0)$, i.e., it is a flat smooth solution, as desired.

For the local problem, one just assumes that $g$ is defined on an open neighborhood defined by $|p^2+q^2|<\epsilon^2$ for some $\epsilon>0$ and the line integrals in the formula will then still work to give the desired solution $f$ on the same neighborhood.

**Added remark:** It's probably worth pointing out that the above result can be applied to prove that the equation $(\star)$ can be solved for a smooth $f$ for a given $g$ if and only if $g$ is smooth and all of its 'balanced partials' (e.g., $g$, $g_{pq}$, $g_{ppqq}$, $g_{pppqqq}$, etc.) vanish at $(p,q) = (0,0)$.  This can be done using the above result plus the theorem that any formal Taylor series is the actual Taylor series of some smooth function.