
function [umean,vmean,rhomean,ui,vi,rhoi,Yi] = BGC_DO_simsvdIC (app);

%-------------------------------------------------------------------------------
% Load realizations.
%-------------------------------------------------------------------------------

% Load ICs.

load (fullfile(app.DOicDir,'000'));

% Set sizes & working storage.

Nu = length(u);
Nv = length(v);
Nrho = length(rho);

var_mat = repmat ( [u./app.sigma_nd(1); v./app.sigma_nd(2); rho./app.sigma_nd(3)] , ...
                   [1 (app.DOnumIC+1)] );

% Load remaining realizations.

for n = 1:app.DOnumIC
    load (fullfile(app.DOicDir,num2str(n,'%3.3d')));
    var_mat(:,n+1) = [u./app.sigma_nd(1); v./app.sigma_nd(2); rho./app.sigma_nd(3)];
end;

%-------------------------------------------------------------------------------
% Create modes from SVD.
%-------------------------------------------------------------------------------

% Remove mean.

var_mean = mean(var_mat,2);

var_mat = var_mat - repmat ( var_mean , [1 (app.DOnumIC+1)] );

var_mat = var_mat;

% Take SVD.

[Umat, Svals, Vmat] = svd (var_mat, 0);

clear var_mat;

% Extract redimensionalized modes.

umean = var_mean((1:Nu)).*app.sigma_nd(1);
vmean = var_mean((Nu+(1:Nv))).*app.sigma_nd(2);
rhomean = var_mean((Nu+Nv+(1:Nrho))).*app.sigma_nd(3);

clear var_mean

mode_ind = 1:app.S(1);

ui = Umat((1:Nu),mode_ind).*app.sigma_nd(1);
vi = Umat((Nu+(1:Nv)),mode_ind).*app.sigma_nd(2);
rhoi = Umat((Nu+Nv+(1:Nrho)),mode_ind).*app.sigma_nd(3);

clear Umat Vmat;

%-------------------------------------------------------------------------------
% Create coefficients from SVD.
%-------------------------------------------------------------------------------

Yi = zeros (app.MC,app.S(1));

for n=1:app.S(1)
    Yi(:,n) = sqrt(Svals(n,n)).*normrnd(0,1,app.MC,1);
end

Yi = Yi.*(sqrt(app.varscale)/ (prod(app.sigma_nd)^(1/length(app.sigma_nd))) ); % if M is realiz matrix, then P=M*M^T/n_r.

%keyboard

