Existence of solution via divergence and curl Let us consider a bounded domain $V$ of $\mathbb{R}^3$ with surface $S$ and let $C$ be a divergence-free vector field and $d$ be some real function on $V$. Let everything be sufficiently well-behaved (smoothness and all that). Now I want to find the solution $F$ to the following
$$
\nabla\cdot F=d, \quad \nabla\times F=C
$$
Where $F$ may possibly satisfy some boundary condition on $S$, e.g., $F\cdot \hat{n}=0$ on $S$. How would one show that such a solution exists (possibly through explicit construction)?
Background. Now from the Wiki page on Helmholtz decomposition, if we were working in $\mathbb{R}^3$ and $C$ decays sufficiently fast at $\infty$, then this $F$ can be easily constructed with the Green functions/Newtonian potential, i.e., if
$$
G(r)=\frac{1}{4\pi |r|}
$$
And $G*f$ denotes the convolution, then we have $F=-\nabla \phi+\nabla \times A$ where $\phi=G*d$ and $A=G*C$. Indeed, notice that
$$
\nabla\times(\nabla\times (G*C))=\nabla(\nabla\cdot(G*C))-\nabla^2(G*C)
$$
The second term is $=C$ since $-\nabla^2 G$ is the dirac-delta function. The first term goes to zero due to the decay of $C$, indeed, approximating $G*C$ by restricting to some bounded domain $V$ of $\mathbb{R}^3$ with surface $S$ and eventually taking $V \to \mathbb{R}^3$, we have
\begin{align}
\nabla\cdot(G*C)&= -\int (\nabla'G(r-r')) C(r') \\
&= -\int\nabla'\cdot(G(r-r')C(r'))\\
&= -\oint_S \hat{n}'\cdot G(r-r')C(r') \to 0
\end{align}
where I have used the fact that $\nabla \cdot C=0$ and $C$ decays sufficiently fast at $\infty$. However, since we're working on a bounded domain, this term doesn't quite disappear, so the method doesn't quite work.

Using the vector calculus identity for curl of a curl, you can rewrite your system as $\nabla^2 F = \nabla d - \nabla\times C$. As $d$ and $C$ are known then you can define a new vector field $a =  \nabla d - \nabla\times C$, which by the Helmholtz decomposition is unrestricted. Then your system is equivalent to the Poisson equation for each component of $F$, $$\nabla^2 F = a, \qquad F\cdot \hat{n}=0,$$ and you now have three equations for three unknowns, plus boundary conditions.
I want to then say that the solution is $F(x) = \int dy \,G(x-y) a(y)$ where $G(x-y)$ is the Green's function for the Laplace operator on your domain, plus a harmonic piece which satisfies $\nabla^2 F =0$. I've only worked with Green's functions on infinite domains though, so I don't want to make a statement about whether there exists a Green's function on an arbitrary bounded domain. But this is surely a standard problem which you could find an answer for online.
