  We consider the one dimensional boundary driven harmonic model and its
continuous version, both introduced in \cite{FGK}. By combining duality and
integrability the authors of \cite{FG} obtained the invariant measures in a
combinatorial representation. Here we give an integral representation of the
invariant measures which turns out to be a convex combination of inhomogeneous
product of geometric distributions for the discrete model and a convex
combination of inhomogeneous product of exponential distributions for the
continuous one. The mean values of the geometric and of the exponential
variables are distributed according to the order statistics of i.i.d. uniform
random variables on a suitable interval fixed by the boundary sources. The
result is obtained solving exactly the stationary condition written in terms of
the joint generating function. The method has an interest in itself and can be
generalized to study other models. We briefly discuss some applications.
