function [time,X,Y,Z,F]=extract_obs(Fnames,P,Vname,rec);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2006 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% [X,Y,Z,F]=extract_obs(Fnames,P,Vname,mask,rec,is3d)                       %
%                                                                           %
% This function extracts requested variable from grided ROMS history        %
% NetCDF file record for 4DVAR twin experiments.                            %
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    Fnames      IO NetCDF file names (structure):                          %
%                  Fnames.inp,      Input history NetCDF grided data        %
%                  Fnames.grd,      Input grid NetCDF file                  %
%    P           Input parameters (structure):                              %
%                  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.full,          extract full grid (interior + boundary) %
%                  P.Zmin,          If < 0, extract from surface to Zmin.   %
%    Vname       NetCDF variable name (string).                             %
%    rec         NetCDF time record to process (integer).                   %
%                                                                           %
% On Output:                                                                %
%                                                                           %
%    time        Extracted time (days).                                     %
%    X           Extracted I-grid coordinates (fractional grid units).      %
%    Y           Extracted J-grid coordinates (fractional grid units).      %
%    Z           Extracted K-grid (>0) or depths (<0) coordinates.          %
%    F           Extracted variable values.                                 %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%----------------------------------------------------------------------------
%  Set input file names.
%----------------------------------------------------------------------------

if (isfield(Fnames,'inp')),
  OBS_inp=deblank(Fnames.inp);
end,
if (isfield(Fnames,'grd')),
  GRD_inp=deblank(Fnames.grd);
  grid=1;
else,
  GRD_inp=OBS_inp;
end,

%----------------------------------------------------------------------------
%  Inquire information of variable to process.
%----------------------------------------------------------------------------

[dnames,dsizes,igrid]=nc_vinfo(Fnames.inp,deblank(Vname));

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

%----------------------------------------------------------------------------
%  Check for Land/Sea masking.
%----------------------------------------------------------------------------

MaskName=[];
if (igrid == 1),
  MaskName='mask_rho';
elseif (igrid == 3),
  MaskName='mask_u';
elseif (igrid == 4),
  MaskName='mask_v';
end,

[varnam,nvars]=nc_vname(GRD_inp);
got_mask=strmatch(MaskName,varnam);

if (got_mask),
  mask=nc_read(GRD_inp,MaskName);
  if (min(min(mask)) == 1)
    got_mask=0;
  end,
end,

%----------------------------------------------------------------------------
%  Read in requested variable from input NetCDF.
%----------------------------------------------------------------------------

f=nc_read(OBS_inp,deblank(Vname),rec);
time=nc_read(OBS_inp,'ocean_time',rec);
time=time/86400;

if (is3d),
  [Lp,Mp,N]=size(f);
else,
  [Lp,Mp]=size(f);
end,

%----------------------------------------------------------------------------
%  Set extraction parameters.
%----------------------------------------------------------------------------

if (isfield(P,'isample'))
  I=P.isample;
else,
  I=1;
end,
if (isfield(P,'jsample'))
  J=P.jsample;
else,
  J=1;
end,

if (is3d),
  if (isfield(P,'kstr'))
    Kstr=P.kstr;
  else,
    Kstr=1;
  end,
  if (isfield(P,'kend'))
    Kend=P.kend;
  else,
    Kend=N;
  end,
end,

if (isfield(P,'full'))
  full=P.full;
else,
  full=0;                % Only process interior points
end,

if (isfield(P,'zmin'))
  Zmin=P.zmin;
else,
  Zmin=0;
end,

%  Check if limiting the extraction by depth, Zmin.

Zlimit=0;
if (is3d),
  if (Zmin < 0),
    Kstr=1;
    Kend=N;
    Zlimit=1;
    Zvar=depths(OBS_inp,GRD_inp,igrid,0,rec);
  end,
end,

%----------------------------------------------------------------------------
%  Extract requested variable data.
%----------------------------------------------------------------------------

%  Set (I,J) ranges for the data extract.

if (full),                % interior and boundary points
  Istr=1;
  Iend=Lp;
  Jstr=1;
  Jend=Mp;
else,                     % interior points
  Istr=2;
  Iend=Lp-1;
  Jstr=2;
  Jend=Mp-1;
end,

if (is3d),

  if (igrid == 1),
    [x,y,z]=ndgrid(0:Lp-1,0:Mp-1,Kstr:Kend);
  elseif (igrid == 3),
    [x,y,z]=ndgrid(1:Lp  ,0:Mp-1,Kstr:Kend);
  elseif (igrid == 4),
    [x,y,z]=ndgrid(0:Lp-1,1:Mp  ,Kstr:Kend);
  end,

%  Extract 3D field from the surface to specified Zmin.

  if (Zlimit),

    F=f(Istr:I:Iend,Jstr:J:Jend,Kstr:Kend);
    X=x(Istr:I:Iend,Jstr:J:Jend,Kstr:Kend);
    Y=y(Istr:I:Iend,Jstr:J:Jend,Kstr:Kend);
    Z=Zvar(Istr:I:Iend,Jstr:J:Jend,Kstr:Kend);

    idpth=find(Zvar < Zmin);
    if (~isempty(idpth)),
      F(idpth)=NaN;
    end,

    if (got_mask),
      Mask=repmat(mask,[1,1,Kstr:Kend]);
      Mask=Mask(Istr:I:Iend,Jstr:J:Jend,Kstr:Kend);
      Mask=Mask(:);
      imsk=find(Mask == 0);
      if (~isempty(imsk)),
        F(imsk)=[];
        X(imsk)=[];
        Y(imsk)=[];
        Z(imsk)=[];
      end,
    end,

    ind=(isnan(F));

    if (~isempty(ind)),
      F(ind)=[];
      X(ind)=[];
      Y(ind)=[];
      Z(ind)=[];
    end,

  else,

%  Extract 3D field at specified vertical level.

    if (Kstr == Kend),

      F=f(Istr:I:Iend,Jstr:J:Jend,Kstr:Kend);
      X=x(Istr:I:Iend,Jstr:J:Jend);
      Y=y(Istr:I:Iend,Jstr:J:Jend);
      Z=z(Istr:I:Iend,Jstr:J:Jend);

      if (got_mask),
        Mask=mask(Istr:I:Iend,Jstr:J:Jend);
        Mask=Mask(:);
        imsk=find(Mask == 0);
         if (~isempty(imsk)),
          F(imsk)=[];
          X(imsk)=[];
          Y(imsk)=[];
          Z(imsk)=[];
        end,
      end,

%  Extract 3D field from vertical levels Kstr to Kend.

    else,

      N=Kend-Kstr+1;
      F=f(Istr:I:Iend,Jstr:J:Jend,1:N);
      X=x(Istr:I:Iend,Jstr:J:Jend,1:N);
      Y=y(Istr:I:Iend,Jstr:J:Jend,1:N);
      Z=z(Istr:I:Iend,Jstr:J:Jend,1:N);

      if (got_mask),
        Mask=repmat(mask,[1,1,N]);
        Mask=Mask(Istr:I:Iend,Jstr:J:Jend,1:N);
        Mask=Mask(:);
        imsk=find(Mask == 0);
        if (~isempty(imsk)),
          F(imsk)=[];
          X(imsk)=[];
          Y(imsk)=[];
          Z(imsk)=[];
        end,
      end,

    end,

  end,

else,

%  Extract 2D field.

    if (igrid == 1),
      [y,x]=meshgrid(0:Mp-1,0:Lp-1);
    elseif (igrid == 3),
      [y,x]=meshgrid(0:Mp-1,1:Lp  );
    elseif (igrid == 4),
      [y,x]=meshgrid(1:Mp  ,0:Lp-1);
    end,

    F=f(Istr:I:Iend,Jstr:J:Jend);
    X=x(Istr:I:Iend,Jstr:J:Jend);
    Y=y(Istr:I:Iend,Jstr:J:Jend);
    Z=zeros(size(F));

    if (got_mask),
      Mask=mask(Istr:I:Iend,Jstr:J:Jend);
      Mask=Mask(:);
      imsk=find(Mask == 0);
      if (~isempty(imsk)),
        F(imsk)=[];
        X(imsk)=[];
        Y(imsk)=[];
        Z(imsk)=[];
      end,
    end,

end,

%  Convert extracted observation to a vector.

F=F(:);
if (igrid == 1),
  X=X(:);
  Y=Y(:);
elseif (igrid == 3),
  X=X(:)-0.5;
  Y=Y(:);
elseif (igrid == 4),
  X=X(:);
  Y=Y(:)-0.5;
end,
Z=Z(:);

return
