Suppose I have a tensor
\begin{gather} \stackrel{\leftrightarrow}{A} = \begin{bmatrix} a_{11}(\vec{r}) & a_{12}(\vec{r}) & a_{13}(\vec{r})\\ a_{21}(\vec{r}) & a_{22}(\vec{r}) & a_{23}(\vec{r})\\ a_{31}(\vec{r}) & a_{32}(\vec{r}) & a_{33}(\vec{r}) \end{bmatrix} \end{gather}
where $\vec{r} = x_{1} \hat{e}_1 + x_{2} \hat{e}_2 + x_{3} \hat{e}_3$
The divergence of this tensor in general curvilinear coordinates is given by
\begin{gather} \nabla^{c} \cdot \stackrel{\leftrightarrow}{A} = \left[ \frac{\partial A_{ij}}{\partial x^{k}} - \Gamma_{ki}^{l} A_{lj} - \Gamma_{kj}^{l} A_{il} \right] g^{ik} \vec{b}^{j}\\ \vec{b}_{i} = \frac{\partial_{x_{i}} \vec{r}}{ \left| \partial_{x_{i}} \vec{r} \right| } \end{gather}
Using Mathematica, I computed the volume integral of the curvilinear divergence for cylindrical coordinates, giving \begin{gather} \iiint \nabla^{c} \cdot \stackrel{\leftrightarrow}{A} dV = \begin{bmatrix} r \int a_{11} d\theta dz + \int a_{12} dr dz + \int r a_{13} dr d\theta - \int a_{22} dr d\theta dz\\ r \int a_{21} d\theta dz + \int a_{22} dr dz + \int r a_{23} dr d\theta + \int a_{12} dr d\theta dz\\ r \int a_{31} d\theta dz + \int a_{32} dr dz + \int r a_{33} dr d\theta \end{bmatrix} \end{gather}
This does not match the traditional Divergence theorem I'm familiar with, or at least it doesn't appear so to me because of the extra triple integral terms $\vec{C}$:
\begin{gather} \iiint \nabla^{c} \cdot \stackrel{\leftrightarrow}{A} dV = \oint A_{ij} n_j \vec{b}_i dS + \vec{C}\\ \vec{C} \neq 0 \end{gather}
What is the correct transformation from the volume integral of the curvilinear divergence to some surface integral?
Let us first consider the invariant form of classical divergence theorem in vector analysis
$$\int_{\Omega}\nabla \cdot {\bf{v}} dV = \int_{\partial \Omega} {\bf{n}} \cdot {\bf{v}}dS \tag{1}$$
For the sake of memorizing, they say that the gradient operator turns into the the unit normal vector. You can choose your vector ${\bf{v}}$ to be
$${\bf{v}} = {\bf{A}} \cdot {\bf{c}} \tag{2}$$
where ${\bf{A}}$ is a second order tensor and ${\bf{c}}$ is a constant vector. Then using $(1)$ and $(2)$ you can prove that
$$\int_{\Omega} \nabla \cdot {\bf{A}} dV = \int_{\partial \Omega} {\bf{n}} \cdot {\bf{A}} dS \tag{3}$$
Note that $(1)$ is a scalar equation while $(2)$ is a vector equation.
Now, you can use $(3)$ to write the divergence theorem in a curve-linear coordinate. So the next step is to compute the divergence of a second order tensor
$$\begin{align} \nabla \cdot {\bf{A}} &= {\bf{g}}^i \partial_{i} \cdot (A_{jk} {\bf{g}}^j \otimes {\bf{g}}^k) \\ &= {\bf{g}}^i \cdot \partial_{i} (A_{jk} {\bf{g}}^j \otimes {\bf{g}}^k) \\ &= \partial_i A_{jk} ({\bf{g}}^i \cdot {\bf{g}}^j) {\bf{g}}^k + A_{jk} ({\bf{g}}^i \cdot \partial_i{\bf{g}}^j) {\bf{g}}^k + A_{jk} ({\bf{g}}^i \cdot {\bf{g}}^j) \partial_i {\bf{g}}^k \\ &= \partial_i A_{jk} g^{ij} {\bf{g}}^k - \Gamma_{il}^{j} A_{jk} ({\bf{g}}^i \cdot {\bf{g}}^l) {\bf{g}}^k - \Gamma_{il}^{k} A_{jk} g^{ij} {\bf{g}}^l \\ &= \partial_i A_{jk} g^{ij} {\bf{g}}^k - \Gamma_{il}^{j} A_{jk} g^{il} {\bf{g}}^k - \Gamma_{il}^{k} A_{jk} g^{ij} {\bf{g}}^l \\ &= \partial_i A_{jk} g^{ij} {\bf{g}}^k - \Gamma_{ij}^{l} A_{lk} g^{ij} {\bf{g}}^k - \Gamma_{ik}^{l} A_{jl} g^{ij} {\bf{g}}^k \\ &= \left[ \partial_i A_{jk} - \Gamma_{ij}^{l} A_{lk} - \Gamma_{ik}^{l} A_{jl} \right] g^{ij} {\bf{g}}^k \\ &= \left[ \partial_i A_{kj} - \Gamma_{ik}^{l} A_{lj} - \Gamma_{ij}^{l} A_{kl} \right] g^{ik} {\bf{g}}^j \end{align} \tag{4}$$
and also we have
$$\begin{align} {\bf{n}} \cdot {\bf{A}} &= n_i {\bf{g}}^i \cdot ( A_{kj} {\bf{g}}^k \otimes {\bf{g}}^j ) \\ &= n_i A_{jk} ( {\bf{g}}^i \cdot {\bf{g}}^k ) {\bf{g}}^j \\ &= n_i A_{jk} g^{ik} {\bf{g}}^j \end{align} \tag{5}$$
and so the final result is
$$\boxed{\int_{\Omega} \left[ \partial_i A_{kj} - \Gamma_{ik}^{l} A_{lj} - \Gamma_{ij}^{l} A_{kl} \right] g^{ik} {\bf{g}}^j dV = \int_{\partial \Omega} n_{i} A_{jk} g^{ik} {\bf{g}}^{j} dS} \tag{6}$$