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.