function [S]=stats(HisFile, RefFile, HisRec, RefRec);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                      %
%  function [S]=stats(HisFile, RefFile, HisRec, RefRec)                %
%                                                                      %
%  This function computes several comparison statistics between two    %
%  ROMS NetCDF files.                                                  %
%                                                                      %
%  On Input:                                                           %
%                                                                      %
%     HisFile    History NetCDF file name.                             %
%     RefFile    Reference NetCDF file name to compare against.        %
%     HisRec     History 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),
  HisRec=1;
end,
is3d=0;
got_mask=0;

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

[dnames,dsizes]=nc_dim(HisFile);

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([' History file:   ', HisFile]);
disp([' Reference file: ', RefFile]);
disp(' ');
disp([' History Record =   ', num2str(HisRec)]);
disp([' Reference Record = ', num2str(RefRec)]);
disp(' ');
disp(' Statistics between History and Reference fields:');
disp(' ');
disp(' Var      RMS with        Correlation');
disp('          Reference       Reference');
disp(' ');

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

Vname='zeta';
if (got_mask),
  H=nc_readmasked(HisFile,Vname,HisRec);
  H=H(:);
  ind=isnan(H);
  if (~isempty(ind)),
    H(ind)=[];
  end,
  R=nc_readmasked(RefFile,Vname,RefRec);
  R=R(:);
  ind=isnan(R);
  if (~isempty(ind)),
    R(ind)=[];
  end,
else
  H=nc_read(HisFile,Vname,HisRec);
  H=H(:);
  R=nc_read(RefFile,Vname,RefRec);
  R=R(:);
end,

%  Compute RMS and pattern correlation coefficient between history and
%  reference field.

Hmean=mean(H);
Hvar=sqrt(mean((H-Hmean).^2));

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

S.zeta_rms=sqrt(mean((H-R).^2));
if (Hvar*Rvar == 0),
  S.zeta_cor=1.0;
else,
  S.zeta_cor=mean((H-Hmean).*(R-Rmean))/(Hvar*Rvar);
end,

disp([' zeta   ', num2str(S.zeta_rms, '%15.8e'), ...
            '  ', num2str(S.zeta_cor, '%15.8e')]);

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

Vname='ubar';
if (got_mask),
  H=nc_readmasked(HisFile,Vname,HisRec);
  H=H(:);
  ind=isnan(H);
  if (~isempty(ind)),
    H(ind)=[];
  end,
  R=nc_readmasked(RefFile,Vname,RefRec);
  R=R(:);
  ind=isnan(R);
  if (~isempty(ind)),
    R(ind)=[];
  end,
else
  H=nc_read(HisFile,Vname,HisRec);
  H=H(:);
  R=nc_read(RefFile,Vname,RefRec);
  R=R(:);
end,

%  Compute RMS and pattern correlation coefficient between history and
%  reference field.

Hmean=mean(H);
Hvar=sqrt(mean((H-Hmean).^2));

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

S.ubar_rms=sqrt(mean((H-R).^2));
if (Hvar*Rvar == 0),
  S.ubar_cor=1.0;
else,
  S.ubar_cor=mean((H-Hmean).*(R-Rmean))/(Hvar*Rvar);
end,

disp([' ubar   ', num2str(S.ubar_rms, '%15.8e'), ...
            '  ', num2str(S.ubar_cor, '%15.8e')]);

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

Vname='vbar';
if (got_mask),
  H=nc_readmasked(HisFile,Vname,HisRec);
  H=H(:);
  ind=isnan(H);
  if (~isempty(ind)),
    H(ind)=[];
  end,
  R=nc_readmasked(RefFile,Vname,RefRec);
  R=R(:);
  ind=isnan(R);
  if (~isempty(ind)),
    R(ind)=[];
  end,
else
  H=nc_read(HisFile,Vname,HisRec);
  H=H(:);
  R=nc_read(RefFile,Vname,RefRec);
  R=R(:);
end,

%  Compute RMS and pattern correlation coefficient between history and
%  reference field.

Hmean=mean(H);
Hvar=sqrt(mean((H-Hmean).^2));

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

S.vbar_rms=sqrt(mean((H-R).^2));
if (Hvar*Rvar == 0),
  S.vbar_cor=1.0;
else,
  S.vbar_cor=mean((H-Hmean).*(R-Rmean))/(Hvar*Rvar);
end,

disp([' vbar   ', num2str(S.vbar_rms, '%15.8e'), ...
            '  ', num2str(S.vbar_cor, '%15.8e')]);

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

if (is3d),
  Vname='u';
  if (got_mask),
    H=nc_readmasked(HisFile,Vname,HisRec);
    H=H(:);
    ind=isnan(H);
    if (~isempty(ind)),
      H(ind)=[];
    end,
    R=nc_readmasked(RefFile,Vname,RefRec);
    R=R(:);
    ind=isnan(R);
    if (~isempty(ind)),
      R(ind)=[];
    end,
  else
    H=nc_read(HisFile,Vname,HisRec);
    H=H(:);
    R=nc_read(RefFile,Vname,RefRec);
    R=R(:);
  end,

%  Compute RMS and pattern correlation coefficient between history and
%  reference field.

  Hmean=mean(H);
  Hvar=sqrt(mean((H-Hmean).^2));

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

  S.u_rms=sqrt(mean((H-R).^2));
  if (Hvar*Rvar == 0),
    S.u_cor=1.0;
  else,
    S.u_cor=mean((H-Hmean).*(R-Rmean))/(Hvar*Rvar);
  end,

  disp([' u      ', num2str(S.u_rms, '%15.8e'), ...
              '  ', num2str(S.u_cor, '%15.8e')]);
end,

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

if (is3d),
  Vname='v';
  if (got_mask),
    H=nc_readmasked(HisFile,Vname,HisRec);
    H=H(:);
    ind=isnan(H);
    if (~isempty(ind)),
      H(ind)=[];
    end,
    R=nc_readmasked(RefFile,Vname,RefRec);
    R=R(:);
    ind=isnan(R);
    if (~isempty(ind)),
      R(ind)=[];
    end,
  else
    H=nc_read(HisFile,Vname,HisRec);
    H=H(:);
    R=nc_read(RefFile,Vname,RefRec);
    R=R(:);
  end,

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

  Hmean=mean(H);
  Hvar=sqrt(mean((H-Hmean).^2));

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

  S.v_rms=sqrt(mean((H-R).^2));
  if (Hvar*Rvar == 0),
    S.v_cor=1.0;
  else,
    S.v_cor=mean((H-Hmean).*(R-Rmean))/(Hvar*Rvar);
  end,

  disp([' v      ', num2str(S.v_rms, '%15.8e'), ...
              '  ', num2str(S.v_cor, '%15.8e')]);
end,

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

if (is3d),
  Vname='temp';
  if (got_mask),
    H=nc_readmasked(HisFile,Vname,HisRec);
    H=H(:);
    ind=isnan(H);
    if (~isempty(ind)),
      H(ind)=[];
    end,
    R=nc_readmasked(RefFile,Vname,RefRec);
    R=R(:);
    ind=isnan(R);
    if (~isempty(ind)),
      R(ind)=[];
    end,
  else
    H=nc_read(HisFile,Vname,HisRec);
    H=H(:);
    R=nc_read(RefFile,Vname,RefRec);
    R=R(:);
  end,

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

  Hmean=mean(H);
  Hvar=sqrt(mean((H-Hmean).^2));

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

  S.temp_rms=sqrt(mean((H-R).^2));
  if (Hvar*Rvar == 0),
    S.temp_cor=1.0;
  else,
    S.temp_cor=mean((H-Hmean).*(R-Rmean))/(Hvar*Rvar);
  end,

  disp([' temp   ', num2str(S.temp_rms, '%15.8e'), ...
              '  ', num2str(S.temp_cor, '%15.8e')]);
end,

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

if (is3d),
  Vname='salt';
  if (got_mask),
    H=nc_readmasked(HisFile,Vname,HisRec);
    H=H(:);
    ind=isnan(H);
    if (~isempty(ind)),
      H(ind)=[];
    end,
    R=nc_readmasked(RefFile,Vname,RefRec);
    R=R(:);
    ind=isnan(R);
    if (~isempty(ind)),
      R(ind)=[];
    end,
  else
    H=nc_read(HisFile,Vname,HisRec);
    H=H(:);
    R=nc_read(RefFile,Vname,RefRec);
    R=R(:);
  end,

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

  Hmean=mean(H);
  Hvar=sqrt(mean((H-Hmean).^2));

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

  S.salt_rms=sqrt(mean((H-R).^2));
  if (Hvar*Rvar == 0),
    S.salt_cor=1.0;
  else,
    S.salt_cor=mean((H-Hmean).*(R-Rmean))/(Hvar*Rvar);
  end,

  disp([' salt   ', num2str(S.salt_rms, '%15.8e'), ...
              '  ', num2str(S.salt_cor, '%15.8e')]);
end,

return
