This can be just a typo, but you have a wrong parametrization of the sphere, it should be
$$
\begin{align}
x & = a r \sin \theta \cos \varphi, \\
y & = b r \sin \theta \sin \varphi, \\
z & = c r \cos \theta.
\end{align}
$$
Your limits for each variable are correct though. Your Jacobian is incorrect because you forgot to take in account the factors $a,b,c$. It should be
$$\frac{\partial (x,y,z)}{\partial (r, \theta, \varphi)} = - abc r^2 \sin \theta.$$
The $-$ sign is because this parametrization of the sphere reverses orientation.
When I edited your post I made sure to clarify some things but I didn't edit a couple of mistakes, which I intend to explain now:
1) I turned $d$'s into $\partial$'s for the Jacobian to correct your notation.
2) The notations $d(x,y,z)$ and $d(r, \theta, \varphi)$ don't make sense, it is best to stick to $dx \, dy \, dz$ and $dr \, d \theta \, d \varphi$.
Your set up for the mass is correct if you fix the Jacobian and add $abc$. The calculation seems to be too (I haven't checked that thoroughly).
I don't understand what you mean by $x_s$. If you want to compute the $x$ coordinate of the center of mass, I assume you are using
$$x_s = \frac{1}{M} \int\limits_{E} x \mu \, dV, \text{ or } x_s M = \int\limits_{E} x \mu \, dV.$$
As you have seen, this is zero, just like the others will be. This has to do with mjqxxxx's comment that the ellipsoid has symmetry about all axis, therefore its center of mass has to be at the origin.
Your divergence integral looks fine, as $\operatorname{div}F=3.$ I don't think you've set up your surface integral correctly. And yes, brackets would aid in clarity! I would use spherical coordinates for the top of the semi-sphere, and polar (cylindrical) coordinates for the base.
You have:
\begin{align*}
\iint_S(\mathbf{F}\cdot\mathbf{n})\,dS&=\iint_{\text{top}}(\mathbf{F}\cdot\mathbf{n})\,dS+\iint_{\text{base}}(\mathbf{F}\cdot\mathbf{n})\,dS\\
&=\int_0^{2\pi}\int_0^{\pi/2}\left[\langle x,y,z\rangle\cdot\underbrace{\frac{\langle x,y,z\rangle}{\sqrt{x^2+y^2+z^2}}}_{\mathbf{n}=\hat{\mathbf{r}}}\right]r^2\sin(\varphi)\,d\varphi\,d\theta \\
&+\int_0^{2\pi}\int_0^2\left[\langle x,y,z\rangle\cdot\underbrace{\langle 0,0,-1\rangle}_{\mathbf{n}}\right]\,r\,dr\,d\theta\\
&=r^3\int_0^{2\pi}\int_0^{\pi/2}\sin(\varphi)\,d\varphi\,d\theta+\underbrace{\int_0^{2\pi}\int_0^2(-z\,r)\,dr\,d\theta}_{=0,\;\text{because}\; z=0}\\
&=16\pi,
\end{align*}
as before. Note that I used the simplification
$$\langle x,y,z\rangle\cdot\frac{\langle x,y,z\rangle}{\sqrt{x^2+y^2+z^2}}=\frac{x^2+y^2+z^2}{\sqrt{x^2+y^2+z^2}}=\frac{r^2}{r}=r,$$
and that $r=2$ on the top surface of the semi-sphere.
Best Answer
Intuivitely you can think of $r^2\sin\varphi$ as the product of $r$ (which is the arc length of a radian of latitude along a meridian of a sphere with radius $r$), and $r\sin\varphi$ (which is the arc length of a radian of longitude along a parallel with the given $\varphi$ of a sphere with radius $r$).
Formally it's just the Jacobian determinant (modulo a sign that I don't care to work out) of the spherical coordinate map $$ \begin{pmatrix} r \\ \varphi \\ \theta \end{pmatrix} \mapsto \begin{pmatrix} r\sin\varphi\cos\theta \\ r \sin\varphi\sin\theta \\ r \cos\varphi \end{pmatrix} $$