On $\mathbb{R}^3 \setminus \{0\}$, define the vector field $$F(x,y,z) = \left(0, \frac{-2yz}{r^4}, \frac{-r^2 + 2y^2}{r^4} \right), $$ where $r = \sqrt{x^2+y^2+z^2}$. Compute the integral of $F$ over the surface $$S: x^2 + y^2 + \frac{z^2}{4} = 1, z > 0. $$ The orientation does not matter.
This question seems evil. How can one compute this in an easy way? Parametrizing the surface with spherical coordinates leads to a very difficult integral. And we cannot apply the divergence theorem, since we would have to close the surface, and then the domain enclosed by it will contain the origin.
Your divergence theorem idea should work, with one small modification. Isolate the singularity at the origin by placing it in a small sphere. The divergence can then be applied to the solid region between the (closed up version of the) surface and the small sphere.
With this setup, the normal to the sphere is a constant multiple of $\langle x,y,z\rangle$, and there is some cancellation when you take the inner product of this normal with your $F$: $$F \cdot \langle x,y,z\rangle = 0 - 2y^2z/r^4 + (-r^2+2y^2)z/r^4 = -z/r^2.$$ This shouldn't be too hard to deal with.