
function [umean,vmean,rhomean,ui,vi,rhoi,Yi] = ...
         BGC_DO_simsvdIC (srcdir,num_realizations,recip_var,num_mode_keep, ...
                          num_MC_realizations,IC_scale)

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

% Load ICs.

load (fullfile(srcdir,'000'));

% Set sizes & working storage.

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

var_mat = repmat ( [u.*sqrt(recip_var(1)); v.*sqrt(recip_var(2)); rho.*sqrt(recip_var(3))] , ...
                   [1 (num_realizations+1)] );

% Load remaining realizations.

for n = 1:num_realizations
    load (fullfile(srcdir,num2str(n,'%3.3d')));
    var_mat(:,n+1) = [u.*sqrt(recip_var(1)); v.*sqrt(recip_var(2)); rho.*sqrt(recip_var(3))];
end;

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

% Remove mean.

var_mean = mean(var_mat,2);

var_mat = var_mat - repmat ( var_mean , [1 (num_realizations+1)] );

var_mat = var_mat .* IC_scale;

% Take SVD.

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

clear var_mat;

% Extract redimensionalized modes.

umean = var_mean((1:Nu))./sqrt(recip_var(1));
vmean = var_mean((Nu+(1:Nv)))./sqrt(recip_var(2));
rhomean = var_mean((Nu+Nv+(1:Nrho)))./sqrt(recip_var(3));

clear var_mean

mode_ind = 1:num_mode_keep;

ui = Umat((1:Nu),mode_ind)./sqrt(recip_var(1));
vi = Umat((Nu+(1:Nv)),mode_ind)./sqrt(recip_var(2));
rhoi = Umat((Nu+Nv+(1:Nrho)),mode_ind)./sqrt(recip_var(3));

clear Umat Vmat;

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

Yi = zeros (num_MC_realizations,num_mode_keep);

for n=1:num_mode_keep
    Yi(:,n) = Svals(n,n).*normrnd(0,1,num_MC_realizations,1);
end

Yi = Yi./sqrt(num_mode_keep);

