To make things easier, first change the $x,y,z$ coordinates to a new system of Cartesian coordinates $x',y',z'$ parallel to the original one but with the origin at $(0,0,-1)$, so $x'=x$, $y'=y$, $z'=z+1$. In the new coordinate system the vector $\vec g$ is:
$$\vec g = \frac{(x',y',z')}{(x^{\prime 2}+y^{\prime 2}+z^{\prime 2})^{3/2}} = \frac{\vec r}{r^3}\,,$$
where $\vec r=(x',y',z')$ is the radial vector field (from the new origin of coordinates) and $r$ its magnitude. This brings up the symmetry of the problem (and its physical connection; $\vec g$ is proportional to the electric field of a static electric charge at the new origin of coordinates, or to the the magnetic field produced by a Dirac's magnetic monopole at that point). Although you can work in Cartesian coordinates, to take advantage of the symmetry it is probably easier to work using standard spherical coordinates $r,\theta,\phi$, and that is what I am going to do, later converting back to Cartesian coordinates. Before we start it is a good idea to check that $\nabla\cdot\vec g=0$, otherwise the problem has no solution. But the divergence of this vector field is well known ($\vec g$ is essentially the gradient of the Green's function of the Poisson equation):
$$ \nabla\cdot\frac{\vec r}{r^3} = 4\pi\delta^3(\vec r)\,.$$
Although outside the origin the divergence is zero, because it is not zero everywhere in $R^3$, it is not possible to find a vector $\vec f$ that is non-singular everywhere $R^3$ minus the origin and such that $\nabla\times\vec f=\vec g$. However it is going to be possible to find $\vec f$ for the required region $z>0$. Note that the solution to this equation is far from unique, because adding to a solution any gradient also gives a solution. We need only to find one solution, all other solution are obtained form this one by adding an arbitrary gradient. To make the calculations easier we can impose some reasonable condition that will restrict the number of available solutions (physicsits call this to "fix the gauge"). For example, if you want to work in Cartesian coordinates you could require your solution to have $f_z=0$ (physicists call this choice the "axial gauge"). Because $\vec g$ is radial, and the radial component of $\vec f$ does not contribute to it, I am going to make the radial component of $\vec f$, $f_r$ equal to zero. With this choice, using the expression of the curl in spherical coordinates we have for $f_\theta$ and $f_\phi$:
$$
\begin{align}
\frac{1}{r\sin{\theta}}\left[\frac{\partial}{\partial\theta}\left(\sin\theta f_\phi\right)-\frac{\partial f_\theta}{\partial\phi}\right] &= \frac{1}{r^2} \\
\frac{\partial}{\partial r}\left(rf_\phi\right) &= 0 \\
\frac{\partial}{\partial r}\left(rf_\theta\right) &= 0
\end{align}
$$
The last two equations show that in this gauge ($f_r=0$), $f_\theta$ and $f_\phi$ have to be proportional to $1/r$, which makes also the radial part of the first equation correct. In fact, with this in mind, it is easily seen by inspection of the first equation that the following choice is a solution:
$$f_\phi = \frac{C-\cos\theta}{r\sin\theta},~~f_\theta = 0\,,$$
where $C$ is some constant. Because we want the solution to be valid if $z>0$ we
must choose $C=1$. In other words, if $\hat u_\phi$ is the unit vector along the $\phi$ direction, then:
$$\vec f = \frac{1-\cos\theta}{r\sin\theta}\hat u_\phi$$
is a solution to the problem in spherical coordinates. The only thing that remains is to convert back to Cartesian coordinates. If $\hat u_x, \hat u_y, \hat u_z$ are the unit vectors along the Cartesian axes (old and new!), and using:
$$\hat u_\phi = -\sin\phi\,\hat u_x+\cos\phi\,\hat u_y$$
and
$$\cos\phi = \frac{x'}{r\sin\theta},~~\sin\phi = \frac{y'}{r\sin\theta},
~~\cos\theta = \frac{z'}{r},~~r^2\sin^2\theta = x'^2+y'^2$$
we get:
$$\vec f = \frac{y'(z'-r)}{r(x'^2+y'^2) }\hat u_x -\frac{x'(z'-r)}{r(x'^2+y'^2) }\hat u_y\,.$$
Reverting to the original $x,y,z$ coordinate system (just replace $x'~y'$ by $x~y$, and $z'$ by $z+1$) we get the final answer:
$$\vec f = \frac{y(z+1-r)\,\hat u_x-x(z+1-r)\,\hat u_y}{r(x^2+y^2)}$$
with
$$ r=\sqrt{x^2+y^2+(z+1)^2}\,.$$
Note that this solution satisfy the physicists "axial gauge" $f_z=0$, although I was not looking for it. Also note that this solution is regular for $z>-1$, as required, and becomes singular in the half line $x=y=0, z<=-1$.
$\nabla^2 =\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial}{\partial r}\right)+\frac{1}{r^2}\frac{\partial^2}{\partial\phi^2}+\frac{\partial^2}{\partial z^2}$ is supposed to be used for a scalar field only.
For the vector field, however, it is easy to check that Laplacian is given by:
$$\nabla^2 \mathbf{u}= \hat{\boldsymbol\rho}\left(\nabla^2 u_r-\frac{u_r}{r^2}-\frac{2}{r^2}\frac{\partial u_\phi}{\partial \phi}\right)+\hat{\boldsymbol\phi}\left(\nabla^2 u_\phi-\frac{u_\phi}{r^2}+\frac{2}{r^2}\frac{\partial u_r}{\partial \phi}\right)+\hat{\mathbf{z}}\nabla^2u_z.$$
For the vector field given, we have $u_r=u_z =0, u_\phi = 1$, and Laplacian will be equal to
$$\nabla^2 \hat{\boldsymbol{\phi}}=-\frac{\hat{\boldsymbol{\phi}}}{r^2}.$$
So, it means that $\boldsymbol\nabla\times(\boldsymbol\nabla\times\mathbf{u})=\frac{\hat{\boldsymbol{\phi}}}{r^2}=\boldsymbol\nabla(\boldsymbol\nabla\cdot\mathbf{u})-\nabla^2\mathbf{u}=-(-\frac{\hat{\boldsymbol{\phi}}}{r^2})$, so the identity holds.
Best Answer
Let's formulate the definition of curl slightly more precisely in the form of a definition/theorem. I'll also not use boldface objects, simply for ease of typing
In the above paragraph, "nice oriented surface" means nice enough so that Stoke's theorem can be applied; for example a smooth two-dimensional, oriented manifold-with-boundary (or however much you wanna weaken the hypothesis... because each book presents it with varying levels of generality... just as long as Stoke's theorem can be applied).
Note that for the above definition of curl to make sense, we have to first show the existence and uniqueness of such a vector field $H$. We shall start by showing the uniqueness of $H$. So, we assume such a $H$ exists, and then prove its components are determined according to the formula $(2)$; this will complete the proof of uniqueness.
Proof of Uniqueness
I'll carry out the computation in detail for proving $H_x = \dfrac{\partial F_z}{\partial y} - \dfrac{\partial F_y}{\partial z}$, and leave the other two to you (it's simply a matter of renaming $x,y,z$). We prove this equality pointwise of course. So, fix a point $p \in A$; then for any $\delta > 0$ such that the closed cube $C_{p,\delta} = p + [-\delta,\delta]^3$ (which is the closed cube centered at $p$ of sidelength $2\delta$) lies entirely inside $A$ (note that since $A$ is open, there are infinitely many such $\delta>0$), we define $M^{\delta} := \{p_1\}\times [p_2-\delta, p_2 + \delta]\times [p_3-\delta, p_3 + \delta]$. This is a piece of a plane which we shall orient so that it has a constant outward normal vector field $n = e_1 \equiv \boldsymbol{i}$. Now, we calculate: note that $\partial M^{\delta}$ has $4$-pieces, and the unit tangent vector along these boundary paths is constant, so (if you're very careful with signs... which I hope I didn't make any sign mistakes), we get \begin{align} \dfrac{1}{|M^{\delta}|} \int_{\partial M^{\delta}}\langle F, \tau \rangle\, dl &= \dfrac{1}{4\delta^2} \bigg[\int_{-\delta}^{\delta} F_2(p_1, p_2 + y, p_3-\delta) - F_2(p_1, p_2 + y, p_3+\delta) \, dy\bigg]\\ &+\dfrac{1}{4\delta^2}\bigg[ \int_{-\delta}^{\delta} F_3(p_1, p_2+\delta, p_3+z) - F_3(p_1, p_2-\delta, p_3+z)\, dz\bigg] \end{align} For each term, we apply the mean-value theorem for integrals (which we can use since everything is continuous), to obtain some $\eta \in [p_2-\delta, p_2+\delta]$ and some $\zeta\in [p_3-\delta, p_3+\delta]$ such that
\begin{align} \dfrac{1}{|M^{\delta}|} \int_{\partial M^{\delta}}\langle F, \tau \rangle\, dl &= \dfrac{1}{4\delta^2}\bigg[2\delta \cdot F_2(p_1, \eta, p_3-\delta) - 2\delta \cdot F_2(p_1, \eta, p_3+\delta)\bigg] \\ &+ \dfrac{1}{4\delta^2}\bigg[2\delta \cdot F_3(p_1, p_2+\delta, \zeta) - 2\delta \cdot F_3(p_1, p_2-\delta, \zeta)\bigg] \\\\ &=\dfrac{F_3(p_1, p_2+\delta, \zeta) - F_3(p_1, p_2-\delta, \zeta)}{2\delta} -\dfrac{F_2(p_1, \eta, p_3+\delta) - F_2(p_1, \eta, p_3-\delta)}{2\delta} \\ &= \dfrac{\partial F_3}{\partial y}(p_1, \alpha, \zeta) - \dfrac{\partial F_2}{\partial z}(p_1, \eta, \beta), \end{align} for some $\alpha\in [p_2-\delta, p_2+\delta], \beta\in [p_3-\delta, p_3+\delta]$, using the mean-value theorem for derivatives (which can certainly be applied since we assumed $F$ is $C^1$).
Quick summary: what we showed so far is that for every $p\in A$ and every $\delta>0$ such that the cube $C_{p,\delta}$ lies inside $A$, if we define $M^{\delta}$ as above to be the plane centered at $p$ with normal pointing in $e_1$ direction, then, there exist points $a_{p,\delta},b_{p,\delta} \in M^{\delta} \subset C_{p,\delta}$ inside the surface (in particular inside the cube) such that
\begin{align} \dfrac{1}{|M^{\delta}|}\int_{\partial M^{\delta}}\langle F, \tau\rangle \, dl &= \dfrac{\partial F_3}{\partial y}(a_{p,\delta}) - \dfrac{\partial F_2}{\partial z}(b_{p,\delta}) \end{align} From here, it is a simple matter of using the continuity of partial derivatives. Here's the full $\epsilon,\delta$ argument to finish it off: let $p\in A$ and $\epsilon> 0$ be arbitrary. By our hypothesis of $(1)$, there is an open $U$ such that blablabla. Now, for this given $\epsilon > 0$, let's choose $\delta > 0$ small enough so that
(so really we have to take a minimum of several $\delta$'s). Then, we choose the oriented plane $M^{\delta}$ as defined above (this plane lies inside $U$ by construction, because of how small $\delta$ is). \begin{align} \left|\left(\dfrac{\partial F_3}{\partial y}(p) - \dfrac{\partial F_2}{\partial z}(p)\right) - H_1(p)\right| &= \left|\left(\dfrac{\partial F_3}{\partial y}(p) - \dfrac{\partial F_2}{\partial z}(p)\right) - \langle H(p), n(p)\rangle\right| \\\\ &\leq \left|\dfrac{\partial F_3}{\partial y}(p) - \dfrac{\partial F_3}{\partial y}(a_{p,\delta})\right| + \left|-\dfrac{\partial F_2}{\partial z}(p) + \dfrac{\partial F_2}{\partial z}(b_{p,\delta})\right|\\ &+ \left|\left(\dfrac{\partial F_3}{\partial y}(a_{p,\delta}) - \dfrac{\partial F_2}{\partial z}(b_{p,\delta})\right) - \dfrac{1}{|M^{\delta}|}\int_{\partial M^{\delta}}\langle F, \tau\rangle\, dl \right| \\ &+ \left|\dfrac{1}{|M^{\delta}|}\int_{\partial M^{\delta}}\langle F, \tau\rangle\, dl - \langle H(p), n(p) \rangle \right| \\\\ &\leq 4\epsilon \end{align} (each absolute value is $\leq \epsilon$ based on everything I've said above, and by the choice of $\delta$). Since the point $p$ and $\epsilon > 0$ are arbitrary, the inequality above shows that \begin{align} H_1 &= \dfrac{\partial F_3}{\partial y} - \dfrac{\partial F_2}{\partial z} \end{align}
As a recap of the proof idea: choose a small plane $M^{\delta}$ with outward normal pointing along $e_1$; it is the flatness of the plane (which is inherently adapted to cartesian coordinates), along with the ease of boundary parametrization which makes the resulting line integral easy to calculate. Then, simply calculate everything, and use the mean-value theorems for derivatives and integrals (this is one way to fill in the gaps for the arugments you typically see in physics texts, which say "let's keep things up to first order only" and where they use $\approx$ everywhere); finally we complete it off with a standard $\epsilon,\delta$ continuity argument.
Some remarks is that for this argument to work I've had to assume $F$ is $C^1$, so that I can apply the mean-value theorems twice, and finally finish it off with a continuity argument. I'm not sure if this proof can be strengthened so that we only have to assume $F$ is differentiable (rather than $C^1$).
Proof of Converse
Now we show the existence of such a vector field $H$; for this we'll show that $\text{curl}F$, defined by equation $(2)$ satisfies the conditions of $(1)$. Like I mentioned in the comments, I'm not sure how to do this without already appealing to Stokes theorem. With Stokes' theorem, this becomes quite simple.
Let $p\in A$, $\epsilon > 0$ and let $S\subset A$ be any "nice surface". Since $\langle\text{curl}(F), n\rangle$ is a continuous function on $S$, there is an open neighbourhood $U$ around $p$ in $S$ such that for all $q\in U$, \begin{align} \left|\langle \text{curl}F(q), n(q)\rangle - \langle \text{curl}F(p), n(p)\rangle\right| & \leq \epsilon \end{align} Now, for any "nice surface" $M\subset U$ (with unit normal being the restriction of the one already on $S$), we have by Stokes theorem: \begin{align} \left|\dfrac{1}{|M|}\int_{\partial M}\langle F,\tau\rangle \, dl - \langle \text{curl} F(p), n(p) \rangle\right| &= \dfrac{1}{|M|}\left|\int_M \langle \text{curl}F, n\rangle \, dA - \int_M \langle \text{curl} F(p), n(p) \rangle\, dA\right| \\ &\leq \dfrac{1}{|M|}\int_M \left|\langle\text{curl }F, n\rangle - \langle\text{curl }F(p), n(p)\rangle \right| \, dA \\ & \leq \dfrac{1}{|M|} \epsilon |M| \\ &= \epsilon. \end{align} This completes the proof of existence.