I'll try to answer clearly, but the detailed calculations are pretty involved and laborious.
Firstly, the Dirac field lagrangian density is assumed to have the following form (this expression gives the usual Dirac equation in curved space, when we do an arbitrary variation of $\Psi$):
\begin{equation}\tag{1}
\mathscr{L} = i \frac{1}{2} \Bigl( \bar{\Psi} \, \Gamma^{\mu} D_{\mu} \Psi - (D_{\mu} \bar{\Psi}) \, \Gamma^{\mu} \Psi \Bigr),
\end{equation}
where $\Gamma^{\mu} = \gamma^a \, e_a^{\mu}$ are the point dependant Dirac matrices (often called "Fock-Ivanenko coefficients", and also "curvy" Dirac matrices, while $\gamma^a$ are the "flat" Dirac matrices). $D_{\mu}$ is the spinor covariant derivative (under local Lorentz transformations):
\begin{equation}\tag{2}
D_{\mu} \Psi = \partial_{\mu} \Psi + \omega_{\mu} \Psi = \partial_{\mu} \Psi - i \frac{1}{4} \, \omega_{\mu}^{\; ab} \, \sigma_{ab} \, \Psi,
\end{equation}
where $\omega_{\mu}^{\; ab}$ is the spin connexion. Notice that (1) is a real number, and it is zero for the "on-shell" spinor fields. The energy-momentum components $T_{\mu \nu}$ are defined by this expression:
\begin{equation}\tag{3}
T_{\mu \nu} \, \delta g^{\mu \nu} = 2 \, \delta \mathscr{L} - g_{\mu \nu} \, \mathscr{L} \, \delta g^{\mu \nu}.
\end{equation}
The last term is 0 for the on-shell field. So we need to compute the explicit variation $\delta\mathscr{L}$ of (1) under an arbitrary variation of the metric: $\delta g^{\mu \nu}$. Since $g^{\mu \nu} = \eta^{ab} \, e_a^{\mu} \, e_b^{\nu}$ and $g_{\mu \nu} = \eta_{ab} \, e_{\mu}^a \, e_{\nu}^b$, we have the following variation terms (these exclude variations that don't come from the metric itself. Infinitesimal local Lorentz transformations, which leave the local metric invariant, are explicitly excluded):
\begin{align}
\delta e_{\lambda}^a &= \frac{1}{2} \, e_{\mu}^a \, g^{\mu \nu} \, \delta g_{\lambda \nu} \equiv -\, \frac{1}{2} \, e_{\mu}^a \, g_{\lambda \nu} \, \delta g^{\mu \nu}, \tag{4} \\[2ex]
\delta e_a^{\mu} &= \frac{1}{2} \, e_a^{\lambda} \, g_{\lambda \nu} \, \delta g^{\mu \nu} \equiv -\, \frac{1}{2} \, e_a^{\lambda} \, g^{\mu \nu} \, \delta g_{\lambda \nu}, \tag{5} \\[2ex]
\delta \Gamma^{\mu} &= \gamma^a \, \delta e_a^{\mu} = \frac{1}{2} \, \Gamma_{\nu} \, \delta g^{\mu \nu}. \tag{6}
\end{align}
Using (2) and (6), the variation of the lagrangian density (1) is this:
\begin{equation}\tag{7}
\delta \mathscr{L} = i \frac{1}{4} \Bigl( \bar{\Psi} \, \Gamma_{\nu} \, D_{\mu} \, \Psi - (D_{\mu} \bar{\Psi}) \, \Gamma_{\nu} \, \Psi \Bigr) \, \delta g^{\mu \nu} + i \frac{1}{2} \bar{\Psi} \{ \, \Gamma^{\mu}, \, \delta\omega_{\mu} \} \Psi,
\end{equation}
where $\{ A, \, B \} = A B + B A$ (this anti-commutator comes from the covariant derivatives of $\Psi$ and $\bar{\Psi}$). The first part is symmetric, since $\delta g^{\mu \nu}$ is symmetric. As for the scalar field, the antisymmetric part is killed by the symmetric variation $\delta g^{\mu \nu}$. So, from (3), we could write the following expression (remember that $\omega_{\mu} = -\, i \frac{1}{4} \, \sigma_{ab} \, \omega_{\mu}^{ab}$):
\begin{equation}
\begin{aligned}
T_{\mu \nu} \, \delta g^{\mu \nu} \equiv 2 \, \delta \mathscr{L} = i \frac{1}{4} \Bigl( \bar{\Psi} \, \Gamma_{\mu} \, D_{\nu} \, \Psi + \bar{\Psi} \, \Gamma_{\nu} \, D_{\mu} \, \Psi &- (D_{\mu} \bar{\Psi}) \, \Gamma_{\nu} \, \Psi - (D_{\nu} \bar{\Psi}) \, \Gamma_{\mu} \, \Psi \Bigr) \, \delta g^{\mu \nu} \\
&+ \frac{1}{4} \bar{\Psi} \{ \, \gamma^c, \, \sigma_{ab} \} \Psi \, e_c^{\mu} \, \delta \omega_{\mu}^{\; ab}.
\end{aligned}
\end{equation}
The last term can be expanded with some efforts (this is the most laborious part!). Using the Dirac algebra, we find the spin density of the spinor field:
\begin{equation}\tag{8}
\frac{1}{4} \bar{\Psi} \{ \, \gamma^c, \, \sigma_{ab} \} \Psi \, e_c^{\mu} = \frac{1}{2} \, \varepsilon_{ab}^{\; \; cd} \, \bar{\Psi} \, \gamma_d \, \gamma^5 \, \Psi \, e_c^{\mu} \equiv S^c_{\; ab} \, e_c^{\mu} = S^{\mu}_{\; ab}.
\end{equation}
This expression is fully anti-symmetric in all the indices (this is crucial for the contraction on $\delta \omega_{\mu}^{\; ab}$).
Now, we have to find the variation of the spin-connexion, and this is where things are getting a bit tricky. We assume the minimal coupling to gravity, i.e no torsion. The connection is symmetric, as in classical general relativity (Levi-Civita connexion, from the symmetric Christoffel symbols):
\begin{equation}\tag{9}
{\omega_{\mu}}^a_{\; b} = e_{\lambda}^a \, \nabla_{\mu} \, e_b^{\lambda}.
\end{equation}
This thing is the simplest expression that defines the covariant derivative (2) under the local Lorentz transformations. This is the place where torsion could enter the scene, but it isn't a necessity (contrary to what many authors may be saying). The variation of (9) above is messy to do, but you could check the following result (using (4) and (5), and using the definition of the Christoffel symbols):
\begin{equation}\tag{10}
S^{\mu}_{\; ab} \, \delta \omega_{\mu}^{\; ab} \equiv {S^{\mu}}_{\lambda}^{\; \nu} \, \delta \Gamma_{\mu \nu}^{\lambda} - \frac{1}{2} \, S^{\mu \kappa \nu} \, \nabla_{\mu} \, \delta g_{\nu \kappa} = 0.
\end{equation}
This cancellation comes directly from the contraction of a fully anti-symmetric expression on two symmetric parts (${S^{\mu}}_{\lambda}^{\; \nu}$ on $\delta \Gamma_{\mu \nu}^{\lambda}$, and $S^{\mu \kappa \nu}$ on $\delta g_{\nu \kappa}$). So the net result is that (3) gives our desired energy-momentum components:
\begin{equation}\tag{11}
T_{\mu \nu} \, \delta g^{\mu \nu} = i \frac{1}{4} \Bigl( \bar{\Psi} \, \Gamma_{\mu} \, D_{\nu} \, \Psi + \bar{\Psi} \, \Gamma_{\nu} \, D_{\mu} \, \Psi - (D_{\mu} \bar{\Psi}) \, \Gamma_{\nu} \, \Psi - (D_{\nu} \bar{\Psi}) \, \Gamma_{\mu} \, \Psi \Bigr) \, \delta g^{\mu \nu},
\end{equation}
so finally:
\begin{equation}\tag{12}
T_{\mu \nu} = i \frac{1}{4} \Bigl( \bar{\Psi} \, \Gamma_{\mu} \, D_{\nu} \, \Psi + \bar{\Psi} \, \Gamma_{\nu} \, D_{\mu} \, \Psi - (D_{\mu} \bar{\Psi}) \, \Gamma_{\nu} \, \Psi - (D_{\nu} \bar{\Psi}) \, \Gamma_{\mu} \, \Psi \Bigr).
\end{equation}
This is the well-known symmetric energy-momentum of the Dirac field, as defined from the canonical approach with the fancy Belinfante-Rosenfeld correction.
Take note that, except for the vanishing of $\mathscr{L}$ in (3) (for on-shell fields), we don't need to assume that the spinor field obeys its equation of motion. The energy-momentum only gets a correction from the symmetric $\mathscr{L} \, g_{\mu \nu}$ if the field is off-shell (maybe is this related to the cosmological constant ?). Finally, $T_{\mu \nu}$ is more complicated if we let torsion to kick in the game: the covariant derivatives get a contorsion tensor, but the full $T_{\mu \nu}$ stays symmetric.
