Divergence theorem in curvilinear coordinates 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?
 A: 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}$$
