Solution (source code)

= Solution

For a nondegenerate surface take $R>0$, and abbreviate $C=\cosh R$, $H=\sinh R$. Close the upper <oblate spheroid> by its equatorial disk $D$, obtaining a half-ellipsoid $V$. The <divergence> is
$$
\nabla\cdot\mathbf A=y^2+2y.
$$
The integral of $2y$ vanishes by reflection symmetry. At height $z$, the cross-section is a disk of radius $q=C\sqrt{1-z^2/H^2}$. Using planar <polar coordinates>,
$$
\int_{x^2+y^2\leq q^2}y^2\,dx\,dy=\int_0^{2\pi}\sin^2\vartheta\,d\vartheta\int_0^q r^3dr=\frac{\pi q^4}{4}.
$$
Consequently
$$
\int_Vy^2\,dV=\frac{\pi C^4}{4}\int_0^H\left(1-\frac{z^2}{H^2}\right)^2dz=\frac{2\pi}{15}C^4H.
$$
The outward <unit normal vector> on $D$ is $-\mathbf e_z$, so its <flux> is $-2\pi C^2$. Applying the <divergence theorem> to the closed surface and subtracting this disk contribution gives
$$
\boxed{\int_S\mathbf A\cdot\mathbf n\,dS=2\pi\cosh^2R+\frac{2\pi}{15}\cosh^4R\sinh R.}
$$
This is a special case of <polynomial flux through a half-ellipsoid>. Including the closing disk is essential: the constant vertical component has zero <divergence> but nonzero flux through the open curved surface.