function [Fmod,Fobs,Ferr]=obs_read(fname,vname,obsvar,record);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2005 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function [Fmod,Fobs,Ferr]=obs_read(fname,vname,obsvar,record)             %
%                                                                           %
% This function reads requested field from 4DVAR observation NetCDF used    %
% in the twin experiments.                                                  &
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    fname       NetCDF file name (character string).                       %
%    vname       NetCDF state variable name to read (character string).     %
%    obsvar      Observation variable to process:                           %
%                  obsvar = 1,   nonlinear model values.                    %
%                  obsvar = 2,   tangent linear model values.               %        
%    record      Optional, time record to read (integer).                   %
%                                                                           %
% On Output:                                                                %
%                                                                           %
%    Fmod        Model field (array).                                       %
%    Fobs        Observation field (array).                                 %
%    Ferr        Observation error (array).                                 %
%    status      Error flag.                                                %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (obsvar == 2),
  model='TLmodel_value';
else,
  model='NLmodel_value';
end,

%----------------------------------------------------------------------------
%  Determine model state variable identifier.
%----------------------------------------------------------------------------

is2d=0;
is3d=0;

switch deblank(vname)
  case 'zeta'
    vid=1;
    is2d=1;
    Ioff=1;
    Joff=1;
  case 'ubar'
    vid=2;
    is2d=1;
    Ioff=0;
    Joff=1;
  case 'vbar'
    vid=3;
    is2d=1;
    Ioff=1;
    Joff=0;
  case 'u'
    vid=4;
    is3d=1;
    Ioff=0;
    Joff=1;
  case 'v'
    vid=5;
    is3d=1;
    Ioff=1;
    Joff=0;
  case 'temp'
    vid=6;
    is3d=1;
    Ioff=1;
    Joff=1;
  case 'salt'
    vid=7;
    is3d=1;
    Ioff=1;
    Joff=1;
end,

%----------------------------------------------------------------------------
%  Set index to extract from the datum dimension.
%----------------------------------------------------------------------------
  
type=nc_read(fname,'obs_type');
ind=find(type == vid);
clear type

%----------------------------------------------------------------------------
%  Extract coordinates.
%----------------------------------------------------------------------------

x=nc_read(fname,'obs_Xgrid');
Im=max(x(ind))+Ioff;
clear x

y=nc_read(fname,'obs_Ygrid');
Jm=max(y(ind))+Joff;
clear y

if (is3d),
  z=nc_read(fname,'obs_Zgrid');
  Km=max(z(ind));
  clear z
end,

%----------------------------------------------------------------------------
%  Extract requested model field.
%----------------------------------------------------------------------------

F=nc_read(fname,model);
f=squeeze(F(record,ind));

clear F

if (is2d),
 Fmod=reshape(f,[Im,Jm]);
else,
 Fmod=reshape(f,[Im,Jm,Km]);
end,
clear f

%----------------------------------------------------------------------------
%  Extract requested observation field.
%----------------------------------------------------------------------------

F=nc_read(fname,'obs_value');
f=F(ind);
clear F

if (is2d),
 Fobs=reshape(f,[Im,Jm]);
else,
 Fobs=reshape(f,[Im,Jm,Km]);
end,
clear f

%----------------------------------------------------------------------------
%  Extract requested observation error.
%----------------------------------------------------------------------------

F=nc_read(fname,'obs_error');
f=F(ind);
clear F

if (is2d),
 Ferr=reshape(f,[Im,Jm]);
else,
 Ferr=reshape(f,[Im,Jm,Km]);
end,
clear f

return
