Working backwards from PSD to possible signal
There are several issues with your code:


The "Gibbs" issue is a non-issue, so you shouldn't need to window.
You are better off generating the complete $-\pi$ to $+\pi$ spectrum that you want, rather than relying on getting the symmetry correct.
The code 

FFT = magnitude-np.mean(magnitude)+phase


is changing the phase of the PSD magnitude by addition. It should be using complex multiplication:

FFT = magnitude * np.exp(1j*phase)



The code below is a modified version that seems to do the trick. The top graph shows the PSD and the (absolute value of the) FFT, showing that they overlap exactly.  The bottom graph shows the generated time domain signal for this particular realization of the phase noise.



UPDATE: I just realized that the original code probably shifted the PSD by $\pi$ because the peak was in the center, rather than at the extremities.





import math
import numpy as np
from matplotlib import pyplot as plt

    if __name__ == '__main__':
    SIGMA = 4.0
    ALPHA = 0.2
    Lc = 25

    L = 2000.0
    N = 1000

    dz = L/N

    i = np.concatenate([np.linspace(0,N/2,N/2+1), np.linspace(-N/2+1,-1,N/2-1)])

    f = i/(N*dz)
    k = 2*np.pi*f
    PSD = lambda x: (2*3*SIGMA**2*Lc)/(1+(x*Lc)**2)**(0.5+ALPHA)

    magnitude = N*np.sqrt(PSD(k))

    phase = 2*np.pi*np.random.randn(N)
    FFT = magnitude * np.exp(1j*phase)

    fig = plt.figure()


    ax1 = fig.add_subplot(211)
    ax1.plot(np.abs(FFT))
    ax1.plot(magnitude)

    ax2 = fig.add_subplot(212)
    ax2.plot(np.fft.ifft(FFT))

    plt.show()