Evaluating the flow rate of water across the net 
A net is dipped in a river. Determine the flow rate of water across
  the net if the velocity vector field for the river is given by
  $\vec{v}=(x-y)\vec{i}+(z+y)\vec{j}+z^2\vec{k}$ and the net is given by
  $y=\sqrt{1-x^2-z^2}, y\geq 0$, oriented in the positive $y$-direction.

My attempt:
We need to find $\iint_{S}\vec{v}\,d\vec{S}$
The surface describes a hemisphere whose bottom is the disk $x^2+z^2\leq 1$ in the $xz$-plane. This sphere can be parameterized as $\vec{r}(\theta,\phi)=\cos\theta \sin\phi\vec{i}+\sin\theta \sin\phi\vec{j}+\cos\phi\vec{k}$ with $0\leq\theta\leq \pi$ and $0\leq\phi\leq \pi$.
Then I calculated $\vec{r}_{\theta}\times\vec{r}_{\phi}$ which turns out to be $-\cos\theta \sin^2\phi\vec{i}-\sin\theta \sin^2\phi\vec{j}-\sin\phi \cos\phi\vec{k}$.
Since we are interested in the positive $y$-direction, our normal vector $\vec{n}=-\vec{r}_{\theta}\times\vec{r}_{\phi}=\cos\theta \sin^2\phi\vec{i}+\sin\theta \sin^2\phi\vec{j}+\sin\phi \cos\phi\vec{k}$ because the $y$-component, which is $\sin\theta \sin^2\phi$ is always positive in the  interval $0\leq\theta\leq \pi$.
Rewriting $\vec{v}$ in terms of our parameters we get, $\vec{v}=(\cos\theta \sin\phi-\sin\theta \sin\phi)\vec{i}+(\cos\phi+\sin\theta \sin\phi)\vec{j}+\cos^2\phi\vec{k}$
So we can express what we are seeking as the double integral $\int^\pi_{0}\int^\pi_{0}\vec{v}\cdot\vec{n}\,d\phi\,d\theta$ which equals $\int^\pi_{0}\int^\pi_{0}((\cos\theta \sin\phi-\sin\theta \sin\phi)\vec{i}+(\cos\phi+\sin\theta \sin\phi)\vec{j}+\cos^2\phi\vec{k})\cdot(\cos\theta \sin^2\phi\vec{i}+\sin\theta \sin^2\phi\vec{j}+\sin\phi \cos\phi\vec{k})\,d\phi\,d\theta$
This looks extremely tedious to deal with. I expanded the brackets but could not see any similar terms that could be simplifed. Am I miscalculating something or is my approach not the best one? This is a freshman year multivariable calculus problem so I don't know any other fancy methods to solve this.
 A: Great News: your calculus is correct top-to-bottom. And yes, this integration is complicated, but in the end, you get a cute answer.
Given:
\begin{align}
Q &= \iint\limits_{S}\vec{v} \cdot d\vec{S}\\
\end{align}
The flow through the net is, after solving that equation using Mathematica,...
SPOILER ALLERT

 $$Q = \dfrac{4\pi}{3}$$

A: The divergence theorem works well here, provided you close the surface beforehand with the disk $S_1$ of equation $x^2+z^2=1$ in the $xz$ plane.
Your flux equals
$$
\iint_S \vec{v}\cdot d\vec{S}+\iint_{S_1} \vec{v}\cdot d\vec{S}-\iint_{S_1} \vec{v}\cdot d\vec{S}\\
=\iint_{S\cup S_1} \vec{v}\cdot d\vec{S}-\iint_{S_1} \vec{v}\cdot d\vec{S}\\
=\iiint_E \nabla \cdot\vec{v} \;dV-\iint_{S_1} \vec{v}\cdot d\vec{S}\\
=\iiint_E (1+1+2z)\; dV + \iint_{S_1}z \;dS
$$
Computations are easier once you are there.
A: I suppose you might try to integrate over the disk, which is the base of your hemisphere...
Take the disc element $dx\,dz$ at coordinates $x,z$ – the corresponding element of the hemisphere is $dS = \frac{dx\,dz}y$ at $(x,y,z)$. The vector $d\vec S$ is parallel to vector $[x, y, z]$, so it is $$d\vec S = x\,dS\,\vec i +  y\,dS\,\vec j + z\,dS\,\vec k$$ because $(x,y,z)$ is on a unit hemisphere, so $x^2+y^2+z^2=1$ is unit length.
