This question has an answer in the language of high level mathematics. Can somebody explain this in the language of vector calculus.
Part I: Let us consider Cartesian coordinate system with origin $O$ and axes $x,y,z$. Let:
$$r=\sqrt{x^2+y^2+z^2}$$
$$\text{and }\vec{V}=0\ (\hat{i}) +\dfrac{\partial}{\partial z} \left( \dfrac{1}{r} \right) (\hat{j}) -\dfrac{\partial}{\partial y} \left( \dfrac{1}{r} \right) (\hat{k})$$
It is obvious that $\dfrac{1}{r}$ is defined everywhere except at the origin.
Now let us take the divergence of $\vec{V}$:
$$\vec{\nabla} \cdot \vec{V}=0$$
Since $\dfrac{1}{r}$ is not defined at the origin, $\vec{\nabla} \cdot \vec{V}=0$ is true everywhere except at the origin.
Since $\vec{\nabla} \cdot \vec{V} \neq 0$ at a point $P (0,0,0)$, will this prevent us from concluding $\vec{V}=\vec{\nabla} \times \vec{U}$ at points other than $P$? Why? Why not?
Part II: If $\vec{\nabla} \cdot \vec{V} \neq 0$ at points on a one dimensional arbitrary curve in space, will this prevent us from concluding $\vec{V}=\vec{\nabla} \times \vec{U}$ at other points not on the curve? Why? Why not?
NOTE - For both Part I and Part II:
If (Why/Why not) is beyond the scope of vector (multivariable) calculus, just reply as yes/no. However please try your best to explain (Why/Why not) in the language of vector (multivariable) calculus.
SEMI ANSWER: Please point out the limitations
I have stumbled upon a derivation in the language of elementary vector calculus. Please point out if there are any limitations in my derivation. In the context of advanced mathematics (de Rham cohomology or Poincare lemma), it seems to me that there are limitations.
Derivation:
To prove: At all points where $\vec{B}$ is defined (whatever be the domain of $\vec{B}$), if $\vec{\nabla} \cdot \vec{B}=0$, then $\vec{B}=\vec{\nabla} \times \vec{A}$
Proof:
At all points where $\vec{B}$ is defined (whatever be the domain of $\vec{B}$): \begin{align} \vec{B} &= B_x\ (\hat{i}) + B_y\ (\hat{j}) + B_z\ (\hat{k})\\ &= B_x\ (\hat{i}) + B_y\ (\hat{j}) + \int^{(x,y,z)}_{(x,y,\infty)} \dfrac{\partial B_z}{\partial z}\ dz\ (\hat{k})\\ &= B_x\ (\hat{i}) + B_y\ (\hat{j}) - \int^{(x,y,z)}_{(x,y,\infty)} \left( \vec{\nabla} \cdot \vec{B} - \dfrac{\partial B_z}{\partial z}\ \right) dz\ (\hat{k})\\ &\text{{Since $\vec{\nabla} \cdot \vec{B}=0$}}\\ &= B_x\ (\hat{i}) + B_y\ (\hat{j}) - \int^{(x,y,z)}_{(x,y,\infty)} \left( \dfrac{\partial B_x}{\partial x} + \dfrac{\partial B_y}{\partial y}\ \right) dz\ (\hat{k})\\ &= B_x\ (\hat{i}) + B_y\ (\hat{j})\ + \left[ \dfrac{\partial}{\partial x} \left(- \int^{(x,y,z)}_{(x,y,\infty)}B_x\ dz \right) -\dfrac{\partial}{\partial y} \left( \int^{(x,y,z)}_{(x,y,\infty)}B_y\ dz \right) \right] (\hat{k})\\ &\text{{By changing the order of integration and differentiation}}\\ \end{align}
At all points where $\vec{B}$ is defined (whatever be the domain of $\vec{B}$), let's define:
$\displaystyle \vec{A}=\int^{(x,y,z)}_{(x,y,\infty)}B_y\ dz\ (\hat{i}) - \int^{(x,y,z)}_{(x,y,\infty)}B_x\ dz\ (\hat{j}) + 0\ (\hat{k}) + \vec{\nabla}f$
where $f$ is an arbitrary function of $(x,y,z)$
Therefore at all points where $\vec{B}$ is defined (whatever be the domain of $\vec{B}$):
\begin{align} \vec{B} &= \left( \dfrac{\partial A_z}{\partial y}-\dfrac{\partial A_y}{\partial z} \right) (\hat{i}) +\left( \dfrac{\partial A_x}{\partial z}-\dfrac{\partial A_z}{\partial x} \right) (\hat{j}) +\left( \dfrac{\partial A_y}{\partial x}-\dfrac{\partial A_x}{\partial y} \right) (\hat{k})\\ &= \vec{\nabla} \times \vec{A} \end{align}
You seem to know that a divergence free field $\vec V$ can be regarded as curl of some other field: There is a field $\vec U$ such that $\vec V={\rm curl}(\vec U)$. This is a consequence of the so-called Poincaré Lemma.
But there is a catch: The Poincaré Lemma guarantees the existence of such a vector potential $\vec U$ only if the domain of $\vec V$ is, e.g., a ball or star shaped. For your field $\vec V$ this is not the case. Therefore we only can say the following: Each point ${\bf p}$ in the punctured space $\dot{\mathbb R}^3:={\mathbb R}^3\setminus\{{\bf 0}\}$ is the center of an open ball $B_r({\bf p})\subset \dot{\mathbb R}^3$ such that within $B_r({\bf p})$ the field $\vec V$ can be written in the form $\vec V={\rm curl}(\vec U)$ for some $\vec U$ defined in $B_r({\bf p})$ only. These local fields $\vec U$ are not uniquely determined, and it is not at all sure whether the implied "integration constants" can be chosen in a coherent way such that we obtain a single field $\vec U_*$, which then would be defined on all of $\dot{\mathbb R}^3$.
Of course, it could be that "by coincidence" the $\vec V$ in your example possesses a global vector potential $\vec U_*$ nevertheless: not by the Poincaré Lemma per se, but because a certain integrability condition (say, the flux of $\vec V$ across a sphere around ${\bf 0}$ should be $=0$) is fulfilled. Consider as an analogue the functions $z\mapsto{1\over z}$ and $z\mapsto{1\over z^2}$ in the punctured complex plane $\dot{\mathbb C}$. Both have local primitives. But one of them does not have a global primitive in $\dot{\mathbb C}$, the other does, the reason being that $$\int_{\partial D}{1\over z}\>dz=2\pi i\ne0\>, \qquad \int_{\partial D}{1\over z^2}\>dz=0\ .$$
Update: I suggest you look at the entry Helmholtz decomposition in Wikipedia and apply Helmholtz' theorem to a large ball minus a tiny ball around the origin.