Synthesize PSD in MatLab

I want to recover a time signals from a given power spectral density, assuming a normal distribution of the original signal:

PSD;                                 % [(m/s)^2/Hz] given spectrum
T = 60;                              % [s] length of original signal
dt = 0.005;                          % [s] time step of original signal
N = T/dt;                            % [-] number of samples
NFFT = 2^nextpow2(N);                % [-] number of bins for FFT
fs = 1/dt;                           % [Hz] sampling frequency
ASD = sqrt(PSD);                     % [(m/s)/sqrt(Hz)] get amplitude spectrum 
omega = 2*pi*rand(NFFT/2,1);         % [rad] generate phase vector
Z = ASD.*exp(1i*omega);              % create complex amplitude vector
Z = [0;Z;flipud(conj(Z))];           % extend to satisfy symmetry
Y = real(ifft(Z));                   % inverse FFT
[PSDY,f] = pwelch(Y,[],[],NFFT,fs);  % generate PSD from Y to compare

The results show a power spectrum several orders of magnitude lower than the original, but the shape matches very good. I guess there is something wrong with the units or there might be a scaling factor missing. I'm not sure about the units of the time signal after ifft, since the amplitude has [(m/s)/sqrt(Hz)].

I believe there are two problems here.
First, I think that the PSD as you define it (or rather, as you use it) is in the wrong units.
When you define the signal as

Z = ASD.*exp(1i*omega); 

then ASD should be in m/s and not in (m/s)/Hz.
So you should do something like that:

ASD = sqrt(PSD*fs/2)

Now, since PSD is in (m/s)^2/Hz, ASD is in units of m/s.

Next, the ifft should be normalised. That is, you should define the Y as

Y = ifft(Z)*sqrt(NFFT);

One more thing, I am not sure if this is on purpose, but the following line

[PSDY,f] = pwelch(Y,[],[],NFFT,fs);

results in Y being divided into 8 parts (with length <NFFT) with 50% overlap. Each part is zero padded to length of NFFT.
A better practice would be to use something like

[PSDY,f] = pwelch(Y,L,L/2,L,fs);

for some L or

    [PSDY,f] = pwelch(Y,NFFT,[],NFFT,fs);

if you insist. To find out more, go to http://www.mathworks.com/help/signal/ref/pwelch.html

In conclusion, this is your (modified) code:

PSD = 5;                                 % [(m/s)^2/Hz] given spectrum
T = 60;                              % [s] length of original signal
dt = 0.005;                          % [s] time step of original signal
N = T/dt;                            % [-] number of samples
NFFT = 2^nextpow2(N);                % [-] number of bins for FFT
fs = 1/dt;                           % [Hz] sampling frequency
ASD = sqrt(PSD*fs/2);                     % [(m/s)] get amplitude spectrum 
omega = 2*pi*rand(NFFT/2,1);         % [rad] generate phase vector
Z = ASD.*exp(1i*omega);              % create complex amplitude vector
Z = [0;Z;flipud(conj(Z))];           % extend to satisfy symmetry
Y = ifft(Z)*sqrt(NFFT);                   % inverse FFT
[PSDY,f] = pwelch(Y,256,128,256,fs);  % generate PSD from Y to compare

which results in

where the blue line is the estimate PSD.