  We present a construction of harmonic functions on bounded domains for the
spectral fractional Laplacian operator and we classify them in terms of their
divergent profile at the boundary. This is used to establish and solve boundary
value problems associated with nonhomogeneous boundary conditions. We provide a
weak-$L^1$ theory to show how problems with measure data at the boundary and
inside the domain are well-posed. We study linear and semilinear problems,
performing a sub- and supersolution method, and we finally show the existence
of large solutions for some power-like nonlinearities.
