Karhunen-Loeve expansions (KLE) of stochastic processes are important tools in mathematics, the sciences, economics, and engineering. However, the KLE is primarily useful for those processes for which we can identify the necessary components, i.e., a set of basis functions, and the distribution of an associated set of stochastic coefficients. Our ability to derive these components explicitly is limited to a handful processes. In this paper we derive all the necessary elements to implement the KLE for a square-integrable Levy process. We show that the eigenfunctions are sine functions, identical to those found in the expansion of a Wiener process. Further, we show that stochastic coefficients have a jointly infinitely divisible distribution, and we derive the generating triple of the first d coefficients. We also show, that, in contrast to the case of the Wiener process, the coefficients are not independent unless the process has no jumps. Despite this, we develop a series representation of the coefficients which allows for simulation of any process with a strictly positive Levy density. We implement our theoretical results by simulating the KLE of a variance gamma process.