Studying zero-coupon bond prices in the CIR (1985) short rate model, $\text{d}r_t=\kappa(\theta-r_t)\text{d}t+\xi\sqrt{r_t}\text{d}W_t$, Hirsa (2013, Section 1.2.6.2) states that the characteristic function of the realised interest rate $R_t=\int_0^t r_s\text{d}s$ is
\begin{align*}
\varphi_{R_t}(u)=\mathbb{E}\left[e^{iuR_t}\right] = A_t(u)e^{B_t(u)r_0},
\end{align*}
where
\begin{align*}
A_t(u) &= \frac{\exp\left(\frac{\kappa^2\theta t}{\xi^2}\right)}{\left(\cosh\left(\frac{1}{2}\gamma t\right)+\frac{\kappa}{\gamma}\sinh\left(\frac{1}{2}\gamma t\right)\right)^{2\kappa\theta/\xi^2}}, \\
B_t(u) &= \frac{2iu}{\kappa+\gamma\coth\left(\frac{1}{2}\gamma t\right)},\\
\gamma &= \sqrt{\kappa^2-2\xi^2iu}.
\end{align*}
As you say, $\mathbb{E}[R_t]$ can be easily computed using Fubini's theorem but from here you also have $$\mathbb{E}[R_t]=-i\varphi_{R_t}'(0).$$ The variance is
\begin{align}
\mathbb{V}\text{ar}[R_t] &= \mathbb{E}[R_t^2] - \mathbb{E}[R_t]^2 \\
&=-\varphi_{R_t}''(0) + \varphi_{R_t}'(0)^2.
\end{align}
Computing these derivatives may be ugly. You could use finite differences instead, $$\mathbb{E}[R_t^2]\approx-\frac{\varphi_{R_t}(-h)-2\varphi_{R_t}(0)+\varphi_{R_t}(h)}{h^2}=\frac{2-\varphi_{R_t}(-h)-\varphi_{R_t}(h)}{h^2}.$$
Note A similar term to $\gamma$ appears in the characteristic function of the log stock price of the Heston (1993) model. One needs to be careful with the sign of the root (``little Heston trap''). I'm not sure whether the same applies here.
