function [F]=extract_var(Fname,Vname,Fvar,P);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2006 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% [F]=extract_var(Fname,Vname,Fvar,P)                                       %
%                                                                           %
% This function extracts field variance at requested observation locations. %
% It is used in 4DVAR twin experiments.                                     %
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    Fname       Input history NetCDF name (string).                        %
%                  Fnames.grd,      Input grid NetCDF file                  %
%                  Fnames.his,      Input history NetCDF file               %
%    Vname       NetCDF variable to process (string).                       %
%    Fvar        Field variance at model grid points (real array).          %
%    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.Zmin,          If < 0, extract from surface to Zmin.   %
%                                                                           %
% On Output:                                                                %
%                                                                           %
%    F           Extracted field variance.                                  %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

if (isfield(Fnames,'his')),
  OBS_his=deblank(Fnames.his);
end,
if (isfield(Fnames,'grd')),
  OBS_grd=deblank(Fnames.grd);
else,
  OBS_grd=OBS_his;
end,

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

[dnames,dsizes,igrid]=nc_vinfo(OBS_his,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(Fname);
got_mask=strmatch(MaskName,varnam);

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

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

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

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,'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_his,OBS_grd,igrid,0,rec);
  end,
end,

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

if (is3d),

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

  if (Zlimit),

    F=Fvar(1:I:Lp,1:J:Mp,Kstr:Kend);

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

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

    ind=(isnan(F));

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

  else,

%  Extract 3D field at specified vertical level.

    if (Kstr == Kend),

      F=Fvar(1:I:Lp,1:J:Mp,Kstr:Kend);

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

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

    else,

      N=Kend-Kstr+1;
      F=Fvar(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)=[];
      end,

    end,

  end,

else,

%  Extract 2D field.

    F=Fvar(1:I:Lp,1:J:Mp);

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

end,

%  Convert extracted observation to a vector.

F=F(:);

return
