Allow me to reformulate. A superfluid is usually described by a complex order parameter $\Psi(x) = \rho(x) e^{i \varphi(x)}$ with $\rho(x) \geq 0$ and $\varphi(x)$ a real field. In the ordered phase, $\rho(x) = \rho_0$, and the phase fluctuates. However, this is not exactly true, because it might happen that $\rho(x_*) = 0 $ for some $x_*$, which is called a vortex, and $x_*$ is the location of the vortex. Note that the decomposition of $\Psi$ into magnitude and phase is only well-defined if the magnitude is non-zero, and thus the phase $\varphi$ is ill-defined at the vortex core!
It is thus dangerous to work with expressions involving $\varphi$ in the presence of vortices, in particular the superfluid velocity is a well-defined concept only away from the vortices. Instead of working with the velocity $v_s$, it is safer to work with the particle current $\vec{j}$. So let's see first what is the connection between the two. In mean-field theory, the superfluid is described by the Lagrangian density
$ L = \frac{i}{2} ( \overline{\Psi} \partial_t \Psi - \Psi \partial_t \overline{\Psi} ) - \frac{1}{2 m} \| \vec{\nabla} \Psi \|^2 - V(|\Psi|) $
from there we can derive the equations of motion
$ i \partial_t \Psi = \frac{1}{2 m} \nabla^2 \Psi + \frac{\partial V}{\partial \overline{\Psi}} $
$ -i \partial_t \overline{\Psi} = \frac{1}{2 m} \nabla^2 \Psi + \frac{\partial V}{\partial \Psi} $
We may now calculate the time-dependence of the particle density $\rho = \sqrt{\overline{\Psi} \Psi}$:
$\partial_t \rho = \frac{1}{2\rho}\left[ \overline{\Psi} \partial_t \Psi + \Psi \partial_t \overline{\Psi} \right] = -\frac{i}{4 m \rho}\left[ \overline{\Psi} \nabla^2 \Psi - \Psi \nabla^2 \overline{\Psi} \right] - \frac{i}{2\rho} \left[ \overline{\Psi} \frac{\partial V}{\partial \overline{\Psi}} - \Psi \frac{\partial V}{\partial \Psi} \right] $
Now in the second term vanishes since $V$ depends only on the magnitude $|\Psi|$, so that we get a current conservation equation:
$ \partial_t \rho = - \vec{\nabla} \cdot \vec{j} \ \ \ , \ \vec{j} = \frac{1}{2m i} \frac{\text{Im}\{ \overline{\Psi} \vec{\nabla} \Psi\}}{\rho} = \frac{\rho(x)}{2m} \vec{\nabla} \varphi$
So this current $\vec{j}$ has a direct physical meaning: it describes the motion of the superfluid density! Now let's assume that $\rho(x) = \rho_0$. Then the expression for the current reduces to
$\vec{j} = \frac{\rho_0}{2 m} \vec{\nabla} \varphi = \rho_0 \vec{v}_s $
So in this limit the current is just the particle density times your superfluid velocity.
If we have now a simply connected region $D$ such that $\rho(x) = \rho_0$ in $D$, then indeed $\vec{\nabla} \times \vec{v}_s = 0$ since it is a gradient. We can use Stokes theorem to infer that all line-integrals inside $D$ of $\vec{v}_s$ vanish: let $\gamma \in D$ be a closed curve bounding an area $A$, then
$ \int_\gamma \vec{v}_s \cdot d \vec{l} = \int_A \vec{\nabla} \times \vec{v}_s \cdot d \vec{\sigma} = 0$
This however does not imply that $\vec{v}_s$ is zero in $D$! for example take the case of $\vec{v}_s$ a constant adnd the path $\gamma$ enclosing a square of side length $L$, where two of the edges are perpendicular to $\vec{v}_s$, so that $\vec{v}_s \cdot d \vec{l} = 0$, and two of the edges are parallel to $\vec{v}_s$, so that $\vec{v}_s \cdot d \vec{l} = \pm \|\vec{v}_s\|$. Then:
$\int_\gamma \vec{v}_s \cdot d \vec{l} = L ( \|\vec{v}_s\| + 0 - \|\vec{v}_s\| + 0) = 0$
as expected, without saying anything about the magnitude of the superfluid velocity.
Now let's consider the case mentioned in your comment, the case that $\gamma$ encloses a circle $C$ of radius $R$ and $\vec{v}_s \cdot d \vec{l} = \|\vec{v}_s\|$ all along $\gamma$, so that
$\int_\gamma \vec{v}_s \cdot d \vec{l} = 2\pi R \|\vec{v}_s\| . $
Now this is impossible if $\vec{v}_s = \frac{\vec{\nabla}\varphi}{2m}$ for a function $\varphi$ defined in the interior of the circle. For then we may calculate the line-integral by evaluating $\varphi$ on the endpoints of $\gamma$, which coincide, and thus give zero. Hence this is only possible if the phase $\varphi$ is not well-defined inside the circle, and this implies that $\Psi(x_*) = 0 $ for some $x_*$ inside the circle, so there has to be a vortex inside! Note that in this case, we should use the well-defined $\vec{j}$ instead of $\vec{v}_s(x) = \frac{\vec{j}(x)}{\rho(x)}$ which is only defined away from the vortex core. Far away from the vortex core, $\rho$ will be constant, so we may write
$\int_\gamma \vec{j} \cdot d \vec{l} = \rho_0 \int_\gamma \vec{v}_s \cdot d \vec{l}$
but if we now use Stokes' theorem, we do not get zero since there is no reason for $\vec{\nabla} \times \vec{j} = 0 $. However, if we insist on using $\vec{v}_s$, which has $\vec{\nabla} \times \vec{v}_s$ everywhere except at the vortex core, we may now longer use Stokes' theorem
$\int_\gamma \vec{v}_s \cdot d \vec{l} \neq \int_C \vec{\nabla} \times \vec{v}_s d \vec{\sigma} $
simply because the right-hand side does not exist! (It is like trying to make sense of the integral $\int_{0}^1 \frac{d x}{x}$.)
