function [status]=obs_error(StateVars,Fnames,P)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2006 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function [status]=obs_error(StateVars,Fnames,P)                           %
%                                                                           %
% This function computes the observation error covariance from gridded      %
% average and history NetCDF files. It is used in 4DVAR twin experiments.   %
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    StateVars   State variable(s) to process, a vector of field indices:   %
%                  StateVar: 1      free-surface  (isFsur)                  %
%                  StateVar: 2      2D u-momentum (isUbar)                  %
%                  StateVar: 3      2D v-momentum (isVbar)                  %
%                  StateVar: 4      3D u-momentum (isUvel)                  %
%                  StateVar: 5      3D v-momentum (isVvel)                  %
%                  StateVar: 6      Temperature   (isTemp)                  %
%                  StateVar: 7      Salinity      (isSalt)                  %
%    Fnames      IO NetCDF file names (structure):                          %
%                  Fnames.avg,      Input average NetCDF file               %
%                  Fnames.grd,      Input grid NetCDF file                  %
%                  Fnames.his,      Input history NetCDF file               %
%                  Fnames.out,      Output 4DVAR observations data          %
%    P           Input parameter structure.                                 %
%                  P.rec,           Time records to process.                %
%                  P.isample,       I-index sampling step.                  %
%                  P.jsample,       J-index sampling step.                  %
%                  P.kstr,          Starting k-level to process.            %
%                  P.kend,          Ending k-level to process.              %
%                  P.Zmin,          If < 0, extract from surface to Zmin.   %
%                                                                           %
% On Output:                                                                %
%                                                                           %
%    status      NetCDF status flag.                                        %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%----------------------------------------------------------------------------
%  Initialize.
%----------------------------------------------------------------------------

WRITE=1;

isFsur=1;
isUbar=2;
isVbar=3;
isUvel=4;
isVvel=5;
isTemp=6;
isSalt=7;

Nstate=7;

if (isfield(P,'rec'))
  Nrecs=length(P.rec);
else,
  Nrecs=1;
end,

%----------------------------------------------------------------------------
% Get file names.
%----------------------------------------------------------------------------

if (isfield(Fnames,'avg')),
  OBS_avg=deblank(Fnames.avg);
end,
if (isfield(Fnames,'his')),
  OBS_his=deblank(Fnames.his);
end,
if (isfield(Fnames,'out')),
  OBS_out=deblank(Fnames.out);
end,

%----------------------------------------------------------------------------
%  Set variables names to process.
%----------------------------------------------------------------------------

Nvars=length(StateVars);
is3d=0;
if (max(StateVars) > 3 ),
  is3d=1;
end,

%  Sort variables to extract in ascending state variable index.

StateVars=sort(StateVars,'ascend');

for ivar=1:Nvars,
  if (StateVars(ivar) == isFsur),
    Vname(ivar,:)='zeta';
    idfld(ivar)=isFsur;
  elseif (StateVars(ivar) == isUbar),
    Vname(ivar,:)='ubar';
    idfld(ivar)=isUbar;
  elseif (StateVars(ivar) == isVbar),
    Vname(ivar,:)='vbar';
    idfld(ivar)=isVbar;
  elseif (StateVars(ivar) == isUvel),
    Vname(ivar,:)='u   ';
    idfld(ivar)=isUvel;
  elseif (StateVars(ivar) == isVvel),
    Vname(ivar,:)='v   ';
    idfld(ivar)=isVvel;
  elseif (StateVars(ivar) == isTemp),
    Vname(ivar,:)='temp';
    idfld(ivar)=isTemp;
  elseif (StateVars(ivar) == isSalt),
    Vname(ivar,:)='salt';
    idfld(ivar)=isSalt;
  end,
end,

%----------------------------------------------------------------------------
%  Compute observation error covariance.
%----------------------------------------------------------------------------

ObsErr=[];

for n=1:Nrecs,
  for ivar=1:Nvars,

    Favg=average(OBS_avg,deblank(Vname(ivar,:)));
    Fvar=variance(OBS_his,deblank(Vname(ivar,:)),Favg);
    Ferr=extract_var(Fnames,deblank(Vname(ivar,:)),Fvar,P);

    ObsErr=[ObsErr; Ferr];

  end,
end,

%----------------------------------------------------------------------------
%  Write observation error covariance into 4DVAR file.
%----------------------------------------------------------------------------

f=nc_read(OBS_out,'obs_error');

if (length(f) == length(ObsErr)),
  status=nc_write(OBS_out,'obs_error',ObsErr);
else,
  error(['OBS_ERROR: inconsitent Ndatum dimension for obs_error. ']);
end,

return
