Because I've gone ahead with trying to parameterize $S$ directly and learned the hard way that the resulting integral is large and annoying to work with, I'll propose a less direct approach.

Rather than compute the surface integral over $S$ straight away, let's close off the hemisphere with the disk $D$ of radius 9 centered at the origin and coincident with the plane $y=0$ . Then by the divergence theorem, since the region $S\cup D$ is closed, we have

$\displaystyle\iint_{S\cup D}\vec F\cdot\mathrm d\vec S=\iiint_R(\nabla\cdot\vec F)\,\mathrm dV$

where $R$ is the interior of $S\cup D$ . $\vec F$ has divergence

$\nabla\cdot\vec F(x,y,z)=\dfrac{\partial(xz)}{\partial x}+\dfrac{\partial(x)}{\partial y}+\dfrac{\partial(y)}{\partial z}=z$

so the flux over the closed region is

$\displaystyle\iiint_Rz\,\mathrm dV=\int_0^\pi\int_0^\pi\int_0^9\rho^3\cos\varphi\sin\varphi\,\mathrm d\rho\,\mathrm d\theta\,\mathrm d\varphi=0$

The total flux over the closed surface is equal to the flux over its component surfaces, so we have

$\displaystyle\iint_{S\cup D}\vec F\cdot\mathrm d\vec S=\iint_S\vec F\cdot\mathrm d\vec S+\iint_D\vec F\cdot\mathrm d\vec S=0$

$\implies\boxed{\displaystyle\iint_S\vec F\cdot\mathrm d\vec S=-\iint_D\vec F\cdot\mathrm d\vec S}$

Parameterize $D$ by

$\vec s(u,v)=u\cos v\,\vec\imath+u\sin v\,\vec k$

with $0\le u\le9$ and $0\le v\le2\pi$ . Take the normal vector to $D$ to be

$\vec s_u\times\vec s_v=-u\,\vec\jmath$

Then the flux of $\vec F$ across $S$ is

$\displaystyle\iint_D\vec F\cdot\mathrm d\vec S=\int_0^{2\pi}\int_0^9\vec F(x(u,v),y(u,v),z(u,v))\cdot(\vec s_u\times\vec s_v)\,\mathrm du\,\mathrm dv$