function [S]=stats_4dvar(HisFile, IniFile, RefFile, IniRec, RefRec);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                      %
%  function [S]=stats_4dvar(HisFile, IniFile, RefFile, IniRec, RefRec) %
%                                                                      %
%  This function computes several statistics between the estimated     %
%  4DVAR initial conditions file and the reference file.               %
%                                                                      %
%  On Input:                                                           %
%                                                                      %
%     HisFile    Nonlinear model history file.                         %
%     IniFile    Initial condition NetCDF file name from 4DVAR.        %
%     RefFile    Reference NetCDF file name.                           %
%     IniRec     Initial conditions NetCDF file record to process.     %
%     RefRec     Reference NetCDF file record to process.              %
%                                                                      %
%  On Output:                                                          %
%                                                                      %
%     S          Statistics (structure array):                         %
%                                                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize.

if (nargin < 4),
  RefRec=1;
end,
if (nargin < 3),
  IniRec=1;
end,
is3d=0;
got_mask=0;

ncid=mexcdf('ncopen',HisFile,'nc_nowrite');
[OBSfile,status]=mexcdf('attget',ncid,-1,'obs_file');
[NRMfile,status]=mexcdf('attget',ncid,-1,'nrm_file');
[status]=mexcdf('ncclose',ncid);

ncid=mexcdf('ncopen',IniFile,'nc_nowrite');
[BCKfile,status]=mexcdf('attget',ncid,-1,'bck_file');
if (isempty(BCKfile)),
  [BCKfile,status]=mexcdf('attget',ncid,-1,'ini_file');
end,
[status]=mexcdf('ncclose',ncid);

Nouter=nc_read(HisFile,'Nouter');
Ninner=nc_read(HisFile,'Ninner');
Hgamma=nc_read(HisFile,'Hgamma');
Vgamma=nc_read(HisFile,'Vgamma');
Hdecay=nc_read(HisFile,'Hdecay');
Vdecay=nc_read(HisFile,'Vdecay');

%-----------------------------------------------------------------------
%  Inquire input NetCDF files.
%-----------------------------------------------------------------------

[dnames,dsizes]=nc_dim(IniFile);

for n=1:length(dsizes),
  name=deblank(dnames(n,:));
  switch name
    case 'xi_rho'
      Lr=dsizes(n);
    case 'xi_u'
      Lu=dsizes(n);
    case 'xi_v'
      Lv=dsizes(n);
    case 'eta_rho'
      Mr=dsizes(n);
    case 'eta_u'
      Mu=dsizes(n);
    case 'eta_v'
      Mv=dsizes(n);
    case 's_rho'
      Nr=dsizes(n);
      is3d=1;
    case 's_rho'
      Nw=dsizes(n);
  end,
end,

disp(' ');
disp([' Initial Conditions file: ', IniFile]);
disp([' Reference file:          ', RefFile]);
disp([' Background file:         ', BCKfile]);
disp([' Normalization file:      ', NRMfile]);
disp([' Observation file:        ', OBSfile]);
disp(' ');
disp([' Hgamma = ', num2str(Hgamma)]);
disp([' Vgamma = ', num2str(Vgamma)]);
disp(' ');
disp([' Hdecay(isFsur) = ', num2str(Hdecay(1)/1000), ' km']);
disp([' Hdecay(isUvel) = ', num2str(Hdecay(4)/1000), ' km']);
disp([' Hdecay(isVvel) = ', num2str(Hdecay(5)/1000), ' km']);
disp([' Hdecay(isTemp) = ', num2str(Hdecay(6)/1000), ' km']);
disp([' Hdecay(isSalt) = ', num2str(Hdecay(7)/1000), ' km']);
disp([' Vdecay(isUvel) = ', num2str(Vdecay(4)), ' m']);
disp([' Vdecay(isVvel) = ', num2str(Vdecay(5)), ' m']);
disp([' Vdecay(isTemp) = ', num2str(Vdecay(6)), ' m']);
disp([' Vdecay(isSalt) = ', num2str(Vdecay(7)), ' m']);
disp(' ');
disp([' Nouter = ', num2str(Nouter)]);
disp([' Ninner = ', num2str(Ninner)]);
disp(' ')
disp(' Statistics with respect Background and Reference fields:');
disp(' ');
disp(' Var      RMS with        RMS with        Correlation     Correlation');
disp('          Background      Reference       Background      Reference');
disp(' ');

%-----------------------------------------------------------------------
%  Process free-surface.
%-----------------------------------------------------------------------

Vname='zeta';
if (got_mask),
  F=nc_readmasked(IniFile,Vname,IniRec);
  F=F(:);
  ind=isnan(F);
  if (~isempty(ind)),
    F(ind)=[];
  end,
  B=nc_readmasked(IniFile,Vname,IniRec+1);
  B=B(:);
  ind=isnan(B);
  if (~isempty(ind)),
    B(ind)=[];
  end,
  R=nc_readmasked(RefFile,Vname,RefRec);
  R=R(:);
  ind=isnan(R);
  if (~isempty(ind)),
    R(ind)=[];
  end,
else
  F=nc_read(IniFile,Vname,IniRec);
  F=F(:);
  B=nc_read(IniFile,Vname,IniRec+1);
  B=B(:);
  R=nc_read(RefFile,Vname,RefRec);
  R=R(:);
end,

%  Compute RMS and pattern correlation coefficien with respect
%  the background and reference field.

Fmean=mean(F);
Fvar=sqrt(mean((F-Fmean).^2));

Bmean=mean(B);
Bvar=sqrt(mean((B-Bmean).^2));

Rmean=mean(R);
Rvar=sqrt(mean((R-Rmean).^2));

S.zeta_Brms=sqrt(mean((F-B).^2));
S.zeta_Bcor=mean((F-Fmean).*(B-Bmean))/(Fvar*Bvar);

S.zeta_Rrms=sqrt(mean((F-R).^2));
S.zeta_Rcor=mean((F-Fmean).*(R-Rmean))/(Fvar*Rvar);

disp([' zeta   ', num2str(S.zeta_Brms, '%15.8e'), ...
            '  ', num2str(S.zeta_Rrms, '%15.8e'), ...
            '  ', num2str(S.zeta_Bcor, '%15.8e'), ...
            '  ', num2str(S.zeta_Rcor, '%15.8e')]);

%-----------------------------------------------------------------------
%  Process 2D U-momentum.
%-----------------------------------------------------------------------

Vname='ubar';
if (got_mask),
  F=nc_readmasked(IniFile,Vname,IniRec);
  F=F(:);
  ind=isnan(F);
  if (~isempty(ind)),
    F(ind)=[];
  end,
  B=nc_readmasked(IniFile,Vname,IniRec+1);
  B=B(:);
  ind=isnan(B);
  if (~isempty(ind)),
    B(ind)=[];
  end,
  R=nc_readmasked(RefFile,Vname,RefRec);
  R=R(:);
  ind=isnan(R);
  if (~isempty(ind)),
    R(ind)=[];
  end,
else
  F=nc_read(IniFile,Vname,IniRec);
  F=F(:);
  B=nc_read(IniFile,Vname,IniRec+1);
  B=B(:);
  R=nc_read(RefFile,Vname,RefRec);
  R=R(:);
end,

%  Compute RMS and pattern correlation coefficien with respect
%  the background and reference field.

Fmean=mean(F);
Fvar=sqrt(mean((F-Fmean).^2));

Bmean=mean(B);
Bvar=sqrt(mean((B-Bmean).^2));

Rmean=mean(R);
Rvar=sqrt(mean((R-Rmean).^2));

S.ubar_Brms=sqrt(mean((F-B).^2));
S.ubar_Bcor=mean((F-Fmean).*(B-Bmean))/(Fvar*Bvar);

S.ubar_Rrms=sqrt(mean((F-R).^2));
S.ubar_Rcor=mean((F-Fmean).*(R-Rmean))/(Fvar*Rvar);

disp([' ubar   ', num2str(S.ubar_Brms, '%15.8e'), ...
            '  ', num2str(S.ubar_Rrms, '%15.8e'), ...
            '  ', num2str(S.ubar_Bcor, '%15.8e'), ...
            '  ', num2str(S.ubar_Rcor, '%15.8e')]);

%-----------------------------------------------------------------------
%  Process 2D V-momentum.
%-----------------------------------------------------------------------

Vname='vbar';
if (got_mask),
  F=nc_readmasked(IniFile,Vname,IniRec);
  F=F(:);
  ind=isnan(F);
  if (~isempty(ind)),
    F(ind)=[];
  end,
  B=nc_readmasked(IniFile,Vname,IniRec+1);
  B=B(:);
  ind=isnan(B);
  if (~isempty(ind)),
    B(ind)=[];
  end,
  R=nc_readmasked(RefFile,Vname,RefRec);
  R=R(:);
  ind=isnan(R);
  if (~isempty(ind)),
    R(ind)=[];
  end,
else
  F=nc_read(IniFile,Vname,IniRec);
  F=F(:);
  B=nc_read(IniFile,Vname,IniRec+1);
  B=B(:);
  R=nc_read(RefFile,Vname,RefRec);
  R=R(:);
end,

%  Compute RMS and pattern correlation coefficien with respect
%  the background and reference field.

Fmean=mean(F);
Fvar=sqrt(mean((F-Fmean).^2));

Bmean=mean(B);
Bvar=sqrt(mean((B-Bmean).^2));

Rmean=mean(R);
Rvar=sqrt(mean((R-Rmean).^2));

S.vbar_Brms=sqrt(mean((F-B).^2));
S.vbar_Bcor=mean((F-Fmean).*(B-Bmean))/(Fvar*Bvar);

S.vbar_Rrms=sqrt(mean((F-R).^2));
S.vbar_Rcor=mean((F-Fmean).*(R-Rmean))/(Fvar*Rvar);

disp([' vbar   ', num2str(S.vbar_Brms, '%15.8e'), ...
            '  ', num2str(S.vbar_Rrms, '%15.8e'), ...
            '  ', num2str(S.vbar_Bcor, '%15.8e'), ...
            '  ', num2str(S.vbar_Rcor, '%15.8e')]);

%-----------------------------------------------------------------------
%  Process 3D U-momentum.
%-----------------------------------------------------------------------

if (is3d),
  Vname='u';
  if (got_mask),
    F=nc_readmasked(IniFile,Vname,IniRec);
    F=F(:);
    ind=isnan(F);
    if (~isempty(ind)),
      F(ind)=[];
    end,
    B=nc_readmasked(IniFile,Vname,IniRec+1);
    B=B(:);
    ind=isnan(B);
    if (~isempty(ind)),
      B(ind)=[];
    end,
    R=nc_readmasked(RefFile,Vname,RefRec);
    R=R(:);
    ind=isnan(R);
    if (~isempty(ind)),
      R(ind)=[];
    end,
  else
    F=nc_read(IniFile,Vname,IniRec);
    F=F(:);
    B=nc_read(IniFile,Vname,IniRec+1);
    B=B(:);
    R=nc_read(RefFile,Vname,RefRec);
    R=R(:);
  end,

%  Compute RMS and pattern correlation coefficien with respect
%  the background and reference field.

  Fmean=mean(F);
  Fvar=sqrt(mean((F-Fmean).^2));

  Bmean=mean(B);
  Bvar=sqrt(mean((B-Bmean).^2));

  Rmean=mean(R);
  Rvar=sqrt(mean((R-Rmean).^2));

  S.u_Brms=sqrt(mean((F-B).^2));
  S.u_Bcor=mean((F-Fmean).*(B-Bmean))/(Fvar*Bvar);

  S.u_Rrms=sqrt(mean((F-R).^2));
  S.u_Rcor=mean((F-Fmean).*(R-Rmean))/(Fvar*Rvar);

  disp([' u      ', num2str(S.u_Brms, '%15.8e'), ...
              '  ', num2str(S.u_Rrms, '%15.8e'), ...
              '  ', num2str(S.u_Bcor, '%15.8e'), ...
              '  ', num2str(S.u_Rcor, '%15.8e')]);
end,

%-----------------------------------------------------------------------
%  Process 3D V-momentum.
%-----------------------------------------------------------------------

if (is3d),
  Vname='v';
  if (got_mask),
    F=nc_readmasked(IniFile,Vname,IniRec);
    F=F(:);
    ind=isnan(F);
    if (~isempty(ind)),
      F(ind)=[];
    end,
    B=nc_readmasked(IniFile,Vname,IniRec+1);
    B=B(:);
    ind=isnan(B);
    if (~isempty(ind)),
      B(ind)=[];
    end,
    R=nc_readmasked(RefFile,Vname,RefRec);
    R=R(:);
    ind=isnan(R);
    if (~isempty(ind)),
      R(ind)=[];
    end,
  else
    F=nc_read(IniFile,Vname,IniRec);
    F=F(:);
    B=nc_read(IniFile,Vname,IniRec+1);
    B=B(:);
    R=nc_read(RefFile,Vname,RefRec);
    R=R(:);
  end,

%  Compute RMS and pattern correlation coefficien with respect
%  the background and reference field.

  Fmean=mean(F);
  Fvar=sqrt(mean((F-Fmean).^2));

  Bmean=mean(B);
  Bvar=sqrt(mean((B-Bmean).^2));

  Rmean=mean(R);
  Rvar=sqrt(mean((R-Rmean).^2));

  S.v_Brms=sqrt(mean((F-B).^2));
  S.v_Bcor=mean((F-Fmean).*(B-Bmean))/(Fvar*Bvar);

  S.v_Rrms=sqrt(mean((F-R).^2));
  S.v_Rcor=mean((F-Fmean).*(R-Rmean))/(Fvar*Rvar);

  disp([' v      ', num2str(S.v_Brms, '%15.8e'), ...
              '  ', num2str(S.v_Rrms, '%15.8e'), ...
              '  ', num2str(S.v_Bcor, '%15.8e'), ...
              '  ', num2str(S.v_Rcor, '%15.8e')]);
end,

%-----------------------------------------------------------------------
%  Process temperature.
%-----------------------------------------------------------------------

if (is3d),
  Vname='temp';
  if (got_mask),
    F=nc_readmasked(IniFile,Vname,IniRec);
    F=F(:);
    ind=isnan(F);
    if (~isempty(ind)),
      F(ind)=[];
    end,
    B=nc_readmasked(IniFile,Vname,IniRec+1);
    B=B(:);
    ind=isnan(B);
    if (~isempty(ind)),
      B(ind)=[];
    end,
    R=nc_readmasked(RefFile,Vname,RefRec);
    R=R(:);
    ind=isnan(R);
    if (~isempty(ind)),
      R(ind)=[];
    end,
  else
    F=nc_read(IniFile,Vname,IniRec);
    F=F(:);
    B=nc_read(IniFile,Vname,IniRec+1);
    B=B(:);
    R=nc_read(RefFile,Vname,RefRec);
    R=R(:);
  end,

%  Compute RMS and pattern correlation coefficien with respect
%  the background and reference field.

  Fmean=mean(F);
  Fvar=sqrt(mean((F-Fmean).^2));

  Bmean=mean(B);
  Bvar=sqrt(mean((B-Bmean).^2));

  Rmean=mean(R);
  Rvar=sqrt(mean((R-Rmean).^2));

  S.temp_Brms=sqrt(mean((F-B).^2));
  S.temp_Bcor=mean((F-Fmean).*(B-Bmean))/(Fvar*Bvar);

  S.temp_Rrms=sqrt(mean((F-R).^2));
  S.temp_Rcor=mean((F-Fmean).*(R-Rmean))/(Fvar*Rvar);

  disp([' temp   ', num2str(S.temp_Brms, '%15.8e'), ...
              '  ', num2str(S.temp_Rrms, '%15.8e'), ...
              '  ', num2str(S.temp_Bcor, '%15.8e'), ...
              '  ', num2str(S.temp_Rcor, '%15.8e')]);
end,

%-----------------------------------------------------------------------
%  Process salinity.
%-----------------------------------------------------------------------

if (is3d),
  Vname='salt';
  if (got_mask),
    F=nc_readmasked(IniFile,Vname,IniRec);
    F=F(:);
    ind=isnan(F);
    if (~isempty(ind)),
      F(ind)=[];
    end,
    B=nc_readmasked(IniFile,Vname,IniRec+1);
    B=B(:);
    ind=isnan(B);
    if (~isempty(ind)),
      B(ind)=[];
    end,
    R=nc_readmasked(RefFile,Vname,RefRec);
    R=R(:);
    ind=isnan(R);
    if (~isempty(ind)),
      R(ind)=[];
    end,
  else
    F=nc_read(IniFile,Vname,IniRec);
    F=F(:);
    B=nc_read(IniFile,Vname,IniRec+1);
    B=B(:);
    R=nc_read(RefFile,Vname,RefRec);
    R=R(:);
  end,

%  Compute RMS and pattern correlation coefficien with respect
%  the background and reference field.

  Fmean=mean(F);
  Fvar=sqrt(mean((F-Fmean).^2));

  Bmean=mean(B);
  Bvar=sqrt(mean((B-Bmean).^2));

  Rmean=mean(R);
  Rvar=sqrt(mean((R-Rmean).^2));

  S.salt_Brms=sqrt(mean((F-B).^2));
  S.salt_Bcor=mean((F-Fmean).*(B-Bmean))/(Fvar*Bvar);

  S.salt_Rrms=sqrt(mean((F-R).^2));
  S.salt_Rcor=mean((F-Fmean).*(R-Rmean))/(Fvar*Rvar);

  disp([' salt   ', num2str(S.salt_Brms, '%15.8e'), ...
              '  ', num2str(S.salt_Rrms, '%15.8e'), ...
              '  ', num2str(S.salt_Bcor, '%15.8e'), ...
              '  ', num2str(S.salt_Rcor, '%15.8e')]);
end,

return
