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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2006 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function [status]=create_obs(StateVars,Fnames,P)                          %
%                                                                           %
% This function extract and creates observations NetCDF file for the 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.inp,      Input observations gridded data         %
%                  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,'inp')),
  OBS_inp=deblank(Fnames.inp);
end,
if (isfield(Fnames,'out')),
  OBS_out=deblank(Fnames.out);
end,

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

Nvars=length(StateVars);
is3d=1;
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,

%----------------------------------------------------------------------------
%  Create NetCDF file.
%----------------------------------------------------------------------------

%  Get title attribute.

[ncid]=mexcdf('ncopen',OBS_inp,'nc_nowrite');
if (ncid == -1),
  error(['CREATE_OBS: ncopen - unable to open file: ', OBS_inp]);
  return
end,

[ncglobal]=mexcdf('parameter','nc_global');
[title]=mexcdf('ncattget',ncid,ncglobal,'title');
lstr=length(title);
if (lstr == 0),
  error('CREATE_OBS: error while reading global attribute: title');
end,

[status]=mexcdf('ncclose',ncid);
if (status == -1),
  error(['CREATE_OBS: ncclose - unable to close NetCDF file', OBS_inp]);
end

%  Get spherical switch.

  spherical=nc_read(OBS_inp,'spherical');

%  Create 4DVAR observations NetCDF file.

[status]=c_observation(OBS_out,is3d,Nrecs,Nstate,title);

%----------------------------------------------------------------------------
%  Extract observations.
%----------------------------------------------------------------------------

time=zeros([Nrecs 1]);
Nobs=zeros([Nrecs 1]);

ObsVal=[];
ObsType=[];
ObsTime=[];

Xobs=[];
Yobs=[];
Zobs=[];

for n=1:Nrecs,
  for ivar=1:Nvars,
    
    [day,X,Y,Z,F]=extract_obs(Fnames,P,deblank(Vname(ivar,:)),P.rec(n));

    Nobs(n)=Nobs(n)+length(F);
    time(n)=day;

    ObsVal=[ObsVal; F];
    ObsType=[ObsType; ones(size(F)).*idfld(ivar)];
    ObsTime=[ObsTime; ones(size(F)).*time(n)];

    Xobs=[Xobs; X];
    Yobs=[Yobs; Y];
    Zobs=[Zobs; Z];

    disp([sprintf('%s',Vname(ivar,:)), ':  ', ...
	  '  F size = ', sprintf('%6g',length(F)) ...
	  '  X size = ', sprintf('%6g',length(X)) ...
	  '  Y size = ', sprintf('%6g',length(Y)) ...
	  '  Z size = ', sprintf('%6g',length(Z))])
  
  end,
end,

%  Set observation error variance to unity. Its values will
%  be determine elsewhere.

ObsErr=ones(size(ObsVal));
ObsVar=ones([1 Nstate]);

% Write observations. Notice that in 3D applications the
% variable "obs_Zgrid" is set to zero sin this will be
% computed by the model. Notice also that "obs_depth"
% can be set in terms of depths (negative values) or
% levels (positive values).

disp(' ')
disp(['obs_type   size = ', sprintf('%10g',length(ObsType))])
disp(['obs_time   size = ', sprintf('%10g',length(ObsTime))])
if (is3d),
  disp(['obs_depth  size = ', sprintf('%10g',length(Zobs))])
  disp(['obs_Zgrid  size = ', sprintf('%10g',length(Zobs))])
end,
disp(['obs_Xgrid  size = ', sprintf('%10g',length(Xobs))])
disp(['obs_Ygrid  size = ', sprintf('%10g',length(Yobs))])
disp(['obs_value  size = ', sprintf('%10g',length(ObsVal))])
disp(['obs_error  size = ', sprintf('%10g',length(ObsErr))])

if (WRITE),

  status=nc_write(OBS_out,'spherical',spherical);
  status=nc_write(OBS_out,'Nobs',Nobs);
  status=nc_write(OBS_out,'survey_time',time);

  status=nc_write(OBS_out,'obs_variance',ObsVar);

  status=nc_write(OBS_out,'obs_type',ObsType);
  status=nc_write(OBS_out,'obs_time',ObsTime);
  if (is3d),
    status=nc_write(OBS_out,'obs_depth',Zobs);
    f=zeros(size(Zobs));
    status=nc_write(OBS_out,'obs_Zgrid',f);
  end,
  status=nc_write(OBS_out,'obs_Xgrid',Xobs);
  status=nc_write(OBS_out,'obs_Ygrid',Yobs);
  status=nc_write(OBS_out,'obs_error',ObsErr);
  status=nc_write(OBS_out,'obs_value',ObsVal);
  
end,

return
