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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2006 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% [X,Y,Z,F]=interp_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.Xgrid,         Model X-grid locations at RH0-points.   %
%                  P.Ygrid,         Model Y-grid locations at RH0-points.   %
%                  P.kstr,          Starting k-level to process.            %
%                  P.kend,          Ending k-level to process.              %
%                  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           Extracter 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 (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,'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,

%----------------------------------------------------------------------------
%  Interpolate requested variable data to provided cruise track(s).
%----------------------------------------------------------------------------

if (is3d),

  [x,y,z]=ndgrid(1:Lp,0:Mp-1,Kstr:Kend);

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

  if (Zlimit),

    F=f(1:I:Lp,1:J:Mp,Kstr:Kend);
    X=x(1:I:Lp,1:J:Mp,Kstr:Kend);
    Y=y(1:I:Lp,1:J:Mp,Kstr:Kend);
    Z=Zvar;

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

    if (got_mask),
      mask=mask(:);
      imsk=find(mask == 0);
      F(imsk)=[];
      X(imsk)=[];
      Y(imsk)=[];
      Z(imsk)=[];
    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(1:I:Lp,1:J:Mp,Kstr:Kend);
      X=x(1:I:Lp,1:J:Mp);
      Y=y(1:I:Lp,1:J:Mp);
      Z=z(1:I:Lp,1:J:Mp);

      if (got_mask),
        mask=mask(:);
        imsk=find(mask == 0);
        F(imsk)=[];
        X(imsk)=[];
        Y(imsk)=[];
        Z(imsk)=[];
      end,

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

    else,

      N=Kend-Kstr+1;
      F=f(1:I:Lp,1:J:Mp,1:N);
      X=x(1:I:Lp,1:J:Mp,1:N);
      Y=y(1:I:Lp,1:J:Mp,1:N);
      Z=z(1:I:Lp,1:J:Mp,1:N);

      if (got_mask),
        mask=repmat(mask,[1,1,N]);
        mask=mask(:);
        imsk=find(mask == 0);
        F(imsk)=[];
        X(imsk)=[];
        Y(imsk)=[];
        Z(imsk)=[];
      end,

    end,

  end,

else,

%  Extract 2D field.

    [y,x]=meshgrid(0:Mp-1,0:Lp-1);

    F=f(1:I:Lp,1:J:Mp);
    X=x(1:I:Lp,1:J:Mp);
    Y=y(1:I:Lp,1:J:Mp);
    Z=zeros(size(F));

    if (got_mask),
      mask=mask(:);
      imsk=find(mask == 0);
      F(imsk)=[];
      X(imsk)=[];
      Y(imsk)=[];
      Z(imsk)=[];
    end,

end,

%  Convert extracted observation to a vector.

F=F(:);
X=X(:);
Y=Y(:);
Z=Z(:);

return
