For the gradient, your mistake is that the components of the gradient vary contravariantly. On top of that, there is a issue with normalisation that I discuss below. I don't know if you are familiar with differential geometry and how it works, but basically, when we write a vector as $v^\mu$ we really are writing its components with respect to a basis.
In differential geometry, vectors are entities which act on functions $f : M \rightarrow \mathbb{R}$ defined on the manifold. Tell me if you want me to elaborate, but this implies that the basis vectors given by some set of coordinates are $\partial_\mu = \frac{\partial}{\partial x^\mu}$ and vary covariantly. Let's name those basis vectors $e_\mu$ to go back to the "familiar" linear algebra notation.
Knowing that, any vector is an invariant which can be written as $\vec{V} = V^\mu \partial_\mu$. The key here is that it is invariant, so it will be the same no matter which coordinate basis you choose.
Now, the gradient is defined in Euclidean space simply as the vector with coordinates $\partial_i f = \partial^if$ where $i = \{x,y,z\}$. Note that in cartesian coordinates covariant and contravariant components are the same. So, the invariant quantity is $\vec{\nabla}f = \partial^ife_i$. Note that, from what we did before, the components of a vector are to be treated as contravariant.
Now, since this expression is invariant we get, in general coordinates $\vec{\nabla}f = \partial^\mu fe_\mu$. So what you are looking for when computing the components is $\partial^\mu f = g^{\mu\nu}\frac{\partial f}{\partial x^\nu}$. This gives $\vec{\nabla}f = \frac{\partial f}{\partial r} e_r + \frac{1}{r^2}\frac{\partial f}{\partial \theta} e_\theta +\frac{1}{r^2 \sin^2\theta}\frac{\partial f}{\partial \phi} e_\phi$. This is still not what we're looking for.
This is due to the fact that the basis vectors are not normalised. Indeed, take a specific vector $e_I$. His components are $\delta^\mu_I$ by definition (it's a basis vector). Then, $|e_I|^2 = e_I^\mu e_I^\nu g_{\mu\nu} = g_{II}$. Great then ! Using $e_i' = e_i/\sqrt{g_{ii}}$ as normalised base vectors, we get the right answer : $\vec{\nabla}f = \frac{\partial f}{\partial r} e_r' + \frac{1}{r}\frac{\partial f}{\partial \theta} e_\theta' +\frac{1}{r \sin\theta}\frac{\partial f}{\partial \phi} e_\phi'$
Let's move onto the divergence. It seems easier since it's a scalar, there's no basis vector to mess around with. Well, actually there are still some problem attached to that. Your formula is right, again, except that when you write the invariant formula $\nabla_\mu V^{\mu}$ you implicitly use the basis that we defined earlier. This means that you are not working on a normalised basis. We know that the $\vec{V}$ you used is $\vec{V} = V^\mu e_\mu = V^\mu \sqrt{g_{\mu\mu}}e_\mu'$. So to compare formula you have to introduce the vector with respect to the normalised coordinates, $A^\mu= V^\mu\sqrt{g_{\mu\mu}}$. I'll let you plug it in your (correct) formula to see that it works.
To conclude, your formula for the curl should be right. Just be careful to use the right normalisations for the vectors and you should be fine (also be careful of the tensorial form of the levi-civita tensor, which involves the determinant of the metric). I have not the faith to do the calculations for you, but you definitely should try it to ensure yourself you understood it well.
P.S: Just for completeness, for the divergence there is a quite useful formula which is also used in Sean Caroll book : $\vec{\nabla}\cdot\vec{V} = \frac{1}{\sqrt{g}}\partial_\mu(\sqrt{g}V^\mu)$, useful when you're too lazy to compute Christoffels.
Magnetic potentials are nowhere near unique, as you have conclusively shown; for more details, look up 'gauge freedom' in your favourite EM textbook or in wikipedia. Imposing the Coulomb gauge condition $\nabla\cdot\mathbf A=0$ reduces the gauge freedom but you can still transform the potential to
$$\mathbf A\mapsto \mathbf A'=\mathbf A+\nabla\psi$$
for any harmonic $\psi$ such that $\nabla^2\psi$, and obtain a different potential which (i) returns the same magnetic field, and (ii) also satisfies the Coulomb gauge condition.
For regular-enough magnetic fields, you can often introduce additional requirements on the vector potential - regularity conditions, and suitable decay at infinity - which can specify it uniquely, but their feasibility depends on the niceness of the magnetic field.
To make this a bit more explicit, you've shown that $\mathbf A_1=\dfrac{\mu_0Iz}{2\pi s}\hat s$ works as a vector potential. However, it is equally easy to check that $\mathbf A_2=-\dfrac{\mu_0I}{2\pi}\ln(s)\hat z$ works equally well: it returns the same field, and it's also in the Coulomb gauge. What gives? Well, the two gauges are related by the transformation
$$
\mathbf A_1
= \mathbf A_2+\nabla\psi
= \mathbf A_2+\nabla\left(\dfrac{\mu_0I}{2\pi}z\ln(s)\right),
$$
where $\psi\propto z\ln(s)$ obeys the Laplace equation. Which one is preferable? Neither, really - they are both singular at the wire ($\mathbf A_1$ more than $\mathbf A_2$), and while $\mathbf A_2$ grows at infinity, the $\sim 1/s$ decay of $\mathbf A_1$ there is probably too slow to be much help. In this situation, the magnetic field is too singular (infinitely thin wire) and contains too much energy (infinitely long wire) for regularity and boundedness conditions to be much help in specifying a unique vector potential.
In general, gauge freedom is something which you can often fix to a large extent but which will always be there lurking in the background. Moreover, there are simply no universally-optimal ways to fix the gauge (so, for instance, the Coulomb gauge $\nabla\cdot \mathbf A=0$ is not Lorentz invariant but the Lorenz gauge $\nabla\cdot \mathbf A+\frac{1}{c^2}\frac{\partial \varphi}{\partial t}=0$ is awkward for non-relativistic work, and so on) so you always need to think of gauge-dependent constructs as temporary, non-unique, and non-physical. The broader-picture answer is to simply let go of the unicity of the magnetic potential.
Best Answer
The vector $\hat \varphi$ is not defined at the origin, because the coordinate transformation $$(x,y) \mapsto (r,\varphi) = \left(\sqrt{x^2 + y^2}, \arctan(y/x)\right)$$ is singular there. Hence your field $\mathbf B$ is singular at the origin.
The theorem that $$\nabla \times \mathbf B = 0 \Rightarrow \oint_C \mathbf B \cdot d\mathbf r = 0$$ requires that the curve $C$ in the line integral be contractible to a point without passing through any singularities. This is not the case for the plane with the origin excluded when the curve winds around the origin.
The singularity of course arises because you have an infinitely thin wire. Try finding the magnetic field for a wire of thickness $R_1$ with uniform current density and taking the limit of $R_1 \to 0$ while the total current remains constant. The curl will be zero outside the wire, but diverge inside the wire, as Maxwell's equations dictate.