  Let $(M,g)$ be a compact Riemann surface with no boundary and
$u=(u_1,...,u_n)$ be a solution of the following singular Liouville system:
\begin{equation*} \Delta_g
u_i+\sum_{j=1}^na_{ij}\rho_j(\frac{h_je^{u_j}}{\int_M h_j
e^{u_j}dV_g}-\frac{1}{vol_g(M)})=\sum_{t=1}^N4\pi \gamma_t(
\delta_{p_t}-\frac{1}{vol_g(M)}), \end{equation*} where $i=1,...,n$,
  $h_1,...,h_n$ are positive smooth functions, $p_1,...,p_N$ are distinct
points on $M$, $\delta_{p_t}$ are Dirac masses, $\rho=(\rho_1,...,\rho_n)$
($\rho_i\ge 0)$ and $(\gamma_{1},...,\gamma_{N})$ ($\gamma_{t}>-1$ ) are
constant vectors. If the coefficient matrix $A=(a_{ij})_{n\times n}$ satisfies
standard assumptions we identify a family of critical hyper-surfaces $\Gamma_k$
for $\rho=(\rho_1,..,\rho_n)$ so that a priori estimate of $u$ holds if $\rho$
is not on any of the $\Gamma_k$s. Thanks to the a priori estimate, a
topological degree for $u$ is well defined for $\rho $ staying between every
two consecutive $\Gamma_k$s. In this article we establish this degree counting
formula which depends only on the Euler Characteristic of $M$ and the location
of $\rho$. Finally if the Liouville system is defined on a bounded domain in
$\mathbb R^2$ with Dirichlet boundary condition, a similar degree counting
formula that depends only on the topology of the domain and the location of
$\rho$ is also determined.
