Find a vector field ${\bf F}$ on $ {\bf R}^3$ with $$\int_S {\bf F}\cdot{\bf n}\ dS=c > 0 \tag{1} $$ where $S$ is any closed surface containing $0$ and ${\bf n}$ is normal
Here there is a solution $\frac{k}{r^3} (x,y,z)$. Note that from divergence therem we know that ${\bf F}$ has divergence $0$ so that $(x,y,-2z)$ is possible. But this solution does not satisfy $(1)$.
So the solution is unique ?
On
Take this one: $\iint_Sr\cdot ndS$ where S is a closed surface as you wish. Indeed, we then have $$\iint_Sr\cdot ndS=\iiint_V\nabla\cdot r dV=\iiint_V\left(\partial_x\text{i}+\partial_y\text{j}+\partial_z\text{k}\right)\cdot(x\text{i}+y\text{j}+z\text{k})dV=3V$$ where $V$ is the volume enclosed by $S$.
On
It is not unique, in fact there are infinitely many possible solutions for any choice of $c>0$. Simply take a vector field $\mathbf{F}$ satisfying: $$\int_S\mathbf{F}\cdot\mathbf{n}dS=k\neq0$$ and then define your new vector field by $\mathbf{G}=\frac{1}{k}\mathbf{F}$.
Now there are infinitely many such $\mathbf{F}$, since by Stokes theorem you we that: $$\int_S\mathbf{F}\cdot\mathbf{n}dS=\int_V\operatorname{div}\mathbf{F}dV$$ where $V$ is the volume contained in $S$, and it's easy to find vector fields with e.g. always positive divergence, such as $\mathbf{F}(x,y,z)=(x,0,0)$.
Also, the vector field you gave in your answer is not defined at the origin.
On
$\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$ $\ds{}$ \begin{align} &\color{#66f}{\large\int_{V}\nabla\cdot\pars{\vec{r} \over r^{3}}\,\dd V} =-\int_{V}\nabla\cdot\nabla\pars{1 \over r}\,\dd V =-\int_{V}\ \overbrace{\ \nabla^{2}\pars{1 \over r}\ } ^{\color{#c00000}{\ds{-\,4\pi\,\delta\pars{\vec{r}}}}}\,\dd V\ =\ 4\pi\int_{V}\delta\pars{\vec{r}}\,\dd V \\[3mm]&=\color{#66f}{\large\left\lbrace\begin{array}{lcl} 4\pi & \mbox{if} & \vec{0}\ \in\ V \\ 0 & \mbox{if} & \vec{0}\ \not\in\ V \end{array}\right.} \end{align}
When they talk about "closed surfaces $S$ containing $0$" they tacitly mean that such $S$ should bound a compact body $B\subset{\mathbb R}^3$ which contains $0$ in its interior. Now we cannot have arbitrarily tiny such surfaces giving a fixed value $c>0$ for the integral in question unless something terrible happens at $0$.
You have remarked that the flow field $${\bf G}(x,y,z):=\left({x\over r^3},{y\over r^3},{z\over r^3}\right)$$ could play a rôle in this question, and you are right: This field is divergence free outside of $0$ and has flow integral $$\int_{S_R}{\bf G}\cdot {\bf n}\>{\rm d}\omega=4\pi$$ when integrated over the sphere $S_R$ of radius $R$ centered at $0$. Using Gauss' theorem, applied to $B\setminus S_R$ with $R\ll1$ it follows that the field $${\bf F}:={c\over4\pi}{\bf G}\tag{1}$$ has flow $c$ through any surface of the kind described in the first paragraph, whence is a solution of the problem.
But this ${\bf F}$ is not the only solution of the problem at hand: Add to ${\bf F}$ an arbitrary $C^1$-vector field ${\bf v}$ which has divergence $0$ in all of ${\mathbb R}^3$. Then ${\bf F}+{\bf v}$ again solves the problem; and it is not difficult to show that all solutions are obtained in this way.