  We study a variational Ginzburg-Landau type model depending on a small
parameter $\varepsilon>0$ for (tangent) vector fields on a $2$-dimensional
Riemannian manifold $S$. As $\varepsilon\to 0$, these vector fields tend to
have unit length so they generate singular points, called vortices, of a
(non-zero) index if the genus $\mathfrak{g}$ of $S$ is different than $1$. Our
first main result concerns the characterization of canonical harmonic unit
vector fields with prescribed singular points and indices. The novelty of this
classification involves flux integrals constrained to a particular
vorticity-dependent lattice in the $2\mathfrak{g}$-dimensional space of
harmonic $1$-forms on $S$ if $\mathfrak{g}\geq 1$. Our second main result
determines the interaction energy (called renormalized energy) between vortex
points as a $\Gamma$-limit (at the second order) as $\varepsilon\to 0$. The
renormalized energy governing the optimal location of vortices depends on the
Gauss curvature of $S$ as well as on the quantized flux. The coupling between
flux quantization constraints and vorticity, and its impact on the renormalized
energy, are new phenomena in the theory of Ginzburg-Landau type models. We also
extend this study to two other (extrinsic) models for embedded hypersurfaces
$S\subset \mathbb{R}^3$, in particular, to a physical model for non-tangent
maps to $S$ coming from micromagnetics.
