We consider the linear system of viscoelasticity with the homogeneous Dirichlet boundary condition. First we prove a Carleman estimate with boundary values of solutions of the viscoelasticity system. Since a solution u under consideration is not assumed to have compact support, in the decoupling of the Lame operator by introducing the divergence and the n-dimensional rotation of u, we have no boundary condition for them, so that we have to carry out arguments by a pseudodifferential operator. Second we apply the Carleman estimate to an inverse source problem of determining a spatially varying factor of the external source in the linear viscoelastitiy by extra Neumann data on a suitable lateral subboundary over a sufficiently long time interval and establish a stability estimate.