Frequency Domain Answer
The flipside of your question would be "What does the value of the spectral density at $0$ represent?"
More generally, the spectral measure $\mu$ of a (mean-zero) covariance-stationary time series is a positive measure on $[-\frac12, \frac12]$.
In cases where it assigns positive mass to ${0}$, the same question can be asked about $\mu({0})$ instead of $f(0)$.
Any covariance-stationary time series $x_t$ has frequency domain representation
$$
x_t = \int_{-\frac12}^{\frac12} e^{-2 \pi i \, \lambda\, t} d \nu (\lambda) \quad \quad (*)
$$
as a stochastic integral in frequency domain. Here $\nu$ is a vector-valued measure on $[-\frac12, \frac12]$ that assigns mean-zero $L^2$-random variables to Borel sets, such that
$$
Var(\nu(B)) = \mu(B), \; \forall \mbox{ Borel } B \subset [-\frac12, \frac12].
$$
Moreover, if $B_1, B_2 \subset [-\frac12, \frac12]$ are disjoint, then $Cov(\nu(B_1), \nu(B_2)) = 0$.
For simplicity, suppose the spectral measure $\mu$ consists of point masses at frequency $\lambda = 0$ and $\lambda = \frac13$. Then the frequency domain representation (*) immediately implies that $\mu(0) = Var(\nu(0))$ is the long run variance of $x_t$.
According to the stochastic integral representation,
$$
x_t = \nu(0) + ( e^{-2 \pi i \, \frac13 \, t} + e^{2 \pi i \, \frac13 \, t} ) \cdot \nu(\frac13) = \nu(0) + 2 \cos \frac{2 \pi}{3} t \cdot \nu(\frac13),
$$
where $\nu(0)$ and $\nu(\frac13)$ are uncorrelated random variables.
Therefore
\begin{align*}
Var(\frac{1}{n} \sum_{t=1}^n x_t)
&= Var(\nu (0) \cdot \frac{1}{n} \sum_{t=1}^n 1 )
+ Var (\nu(\frac13) \cdot \frac{2}{n} \sum_{t=1}^n \cos \frac{2 \pi}{3} t ) \\
&\rightarrow Var( \nu (0) ), \, \mbox{ as } n \rightarrow \infty.
\end{align*}
Empirically, $Var( \nu (0) ) = \mu(0)$ is the intensity, as measured by variance, of the constant part of the series (at frequency $\lambda = 0$). In the long-run, it dominates contributions from higher frequencies, making it the long-run variance.
(When the spectral measure is not absolutely continuous, the ACF is not summable and the normalization factor for long-run variance is $\frac{1}{n}$, rather than $\frac{1}{\sqrt{n}}$.)
The situation where there is a spectral density $f$ is similar. Take a discrete approximation $\mu(0, \epsilon) \approx f(0)\cdot \epsilon$, with $\epsilon = \frac{1}{\sqrt{n}}$. Then the dominating term in
$$
Var(\frac{1}{\sqrt{n}} \sum_{t=1}^n x_t )
$$
would be
$$
Var( \frac{1}{\sqrt{n}} \cdot n \nu (0, \epsilon)) = f(0),
$$
i.e. $f(0)$ is equal to long-run variance.
This is not surprising. From the frequency domain perspective, the long run variance is the variance of the series obtained by applying the two-sided filter $(\cdots,1,1,1,\cdots)$, but the Fourier transform of this filter is precisely point-mass at $\lambda = 0$.
Time Domain Answer
Suppose $x_t$ has causal MA representation
$$
x_t = \Psi (L) \epsilon_t,
$$
where $\Psi(L) = \sum_{h \geq 0} \psi_h L^h$ with $\sum_{h \geq 0} |\psi_h| < \infty.$
Then $f(0) = \Psi(1)^2$.
On the other hand, we can write the lag polynomial $\Psi(L)$ as
$$
\Psi(L) = \Psi(1) + \Psi(L) - \Psi(1) = \Psi(1) + (1 - L)\alpha(L)
$$
where $\Psi(L) - \Psi(1) = (1 - L)\alpha(L)$. This is the so-called "Beveridge-Nelson decomposition". Therefore
$$
x_t = \Psi (1) \epsilon_t + \eta_t
$$
where $\eta_t = (1 - L)\alpha(L) \epsilon_t$ is an over-differenced $I(-1)$ series.
Any over-differenced series has long-run variance $0$. Therefore
the long-run variance of $x_t$ is same as that of $\Psi (1) \epsilon_t$, which is
$\Psi (1)^2 = f(0)$.
Comment
The frequency domain answer says that the long run variance is equal to the intensity of the zero-frequency (constant) part of the series, which is $f(0)$.
The time domain answer is obtained via a different decomposition. The series is decomposed into a white noise part and an over-differenced part. Only the white noise part contributes to long run variance, and it has variance $\sum_h \gamma_h$, which is $\Psi(1)^2 = f(0)$.
