Well, two and a half years later, I stumbled across an answer. The following is a paraphrasing of Petersen's "Riemannian Geometry," page 40.
We work on a Riemannian manifold $(M,\langle \cdot, \cdot \rangle)$ with Levi-Civita connection $\nabla$. Let $\text{Ric}$ denote the symmetric $(1,1)$-tensor $\text{Ric}(v) = \sum R(v,e_i)e_i$, where $\{e_i\}$ is an orthonormal basis. Let $S$ denote the scalar curvature, i.e. $S = \text{tr}(\text{Ric}) = \sum \langle \text{Ric}(e_i), e_i \rangle$. We note that $\nabla S = dS$.
Prop: $dS = 2\,\text{div}(\text{Ric})$.
Proof: Let $\{E_i\}$ be an orthonormal frame at $p \in M$ with $\nabla E_i|_p = 0$. Let $X$ be a vector field with $\nabla X |_p = 0$. From the second Bianchi identity:
$$\begin{align*}
dS(X)(p) = XS(p) & = X \sum \langle \text{Ric}(E_i), E_i \rangle \\
& = X \sum \langle R(E_i, E_j)E_j, E_i\rangle \\
& = \sum \langle \nabla_X [R(E_i, E_j)E_j], E_i \rangle \\
& = \sum \langle (\nabla_X R)(E_i, E_j)E_j, E_i \rangle \\
& = -\sum \langle (\nabla_{E_j}R)(X, E_i)E_j, E_i \rangle - \sum \langle (\nabla_{E_i}R)(E_j,X)E_j, E_i \rangle \\
& = -\sum (\nabla_{E_j}R)(X, E_i, E_j, E_i) - \sum (\nabla_{E_i}R)(E_j, X, E_j, E_i) \\
& = \sum (\nabla_{E_j}R)(E_j, E_i, E_i, X) + \sum (\nabla_{E_j} R)(E_i, E_j, E_j, X) \\
& = 2 \sum (\nabla_{E_j}R)(E_j, E_i, E_i, X) \\
& = 2 \sum \nabla_{E_j}(R(E_j, E_i, E_i, X) \\
& = 2 \sum \nabla_{E_j} \langle \text{Ric}(E_j), X \rangle \\
& = 2 \sum \nabla_{E_j} \langle \text{Ric}(X), E_j \rangle \\
& = 2 \sum \langle \nabla_{E_j}(\text{Ric}(X)), E_j) \rangle \\
& = 2 \sum \langle (\nabla_{E_j}\text{Ric})(X), E_j \rangle \\
& = 2\,\text{div}(\text{Ric})(X)_p.
\end{align*}$$
Best Answer
First note that
$$\nabla^l A_{jl}=\nabla^l R_{jl}-\nabla^l \frac{R}{2(n-1)}g_{jl}=\nabla^l R_{jl}- \frac{1}{2(n-1)}\nabla_j R=\frac{1}{2}\nabla_j R- \frac{1}{2(n-1)}\nabla_j R=\frac{n-2}{2(n-1)}\nabla_j R \qquad(1)$$
the last equality is obtained by twice contracted second Bianchi identity ( that is $2\nabla^jR_{jk}=\nabla_k R$)
Now we calculate the divergence of the Weyl tensor:
$$g^{ls}\nabla _sW_{ijkl}=g^{ls}\nabla_s R_{ijkl}-\frac{1}{n-2}(g_{ik}g^{ls}\nabla_s A_{jl}-g_{il}g^{ls}\nabla_s A_{jk}-g_{jk}g^{ls}\nabla_s A_{il}+g_{jl}g^{ls}\nabla_s A_{ik})$$
$$=\nabla_iR_{jk}-\nabla_jR_{ik}-\frac{1}{n-2}(g_{ik}\nabla^l A_{jl}-g_{il}\nabla^l A_{jk}-g_{jk}\nabla^l A_{il}+g_{jl}\nabla^l A_{ik})$$
$$=\nabla_iR_{jk}-\nabla_jR_{ik}-\frac{1}{n-2}(g_{ik}\nabla^l A_{jl}-\nabla_i A_{jk}-g_{jk}\nabla^l A_{il}+\nabla_j A_{ik})$$
We know $C_{ijk}=\nabla_{i}A_{jk}-\nabla_{j}A_{ik}$, then
$$=\nabla_iR_{jk}-\nabla_jR_{ik}-\frac{1}{n-2}(g_{ik}\nabla^l A_{jl}-g_{jk}\nabla^l A_{il}-C_{ijk})$$
Now by considering (1):
$$=\nabla_iR_{jk}-\nabla_jR_{ik}-\frac{1}{n-2}(g_{ik}\frac{n-2}{2(n-1)}\nabla_j R-g_{jk}\frac{n-2}{2(n-1)}\nabla_i R-C_{ijk})$$
$$=\nabla_iR_{jk}-\nabla_jR_{ik}-\frac{1}{2(n-1)}g_{ik}\nabla_j R+\frac{1}{2(n-1)}g_{jk}\nabla_i R+\frac{1}{n-2}C_{ijk}$$
$$=\nabla_i(R_{jk}+\frac{R}{2(n-1)}g_{jk})-\nabla_j(R_{ik}+\frac{R}{2(n-1)}g_{ik})+\frac{1}{n-2}C_{ijk}$$
By using $\frac{R}{2(n-1)}g_{ij}=R_{ij}-A_{ij}$:
$$=\nabla_i(2R_{jk}-A_{jk})-\nabla_j(2R_{ik}-A_{ik})+\frac{1}{n-2}C_{ijk}$$
$$=2\nabla_i R_{jk}-2\nabla_j R_{ik}-\nabla_iA_{ik}+\nabla_jA_{ik}+\frac{1}{n-2}C_{ijk}=-C_{ijk}+\frac{1}{n-2}C_{ijk}=-\frac{n-3}{n-2}C_{ijk}$$
Note that $2\nabla_i R_{jk}-2\nabla_j R_{ik}=0$. So we obtain