function [S]=checkpointing(Cname, Hname, FULL);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Copyright (c) 2005 ROMS/TOMS Adjoint Group                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function [S]=checkpointing(Cname, Hname, FULL);                           %
%                                                                           %
% This routine extract state variables from GST checkpointing file.         %
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    Cname       Checkpointing NetCDF file name (character string).         %
%    Hname       History NetCDF file name (character string).               %
%    FULL        Switch to process full grid or interior point (optional).  %
%                                                                           %
% On Output:                                                                %
%                                                                           %
%    S           State structure array with several quantities.             %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (nargin < 3),
  FULL=0;
end,
SOLVE3D=0;

%----------------------------------------------------------------------------
%  Inquire about grid dimensions.
%----------------------------------------------------------------------------

[dnames,dsizes]=nc_dim(Hname);

ndims=length(dsizes);
for n=1:ndims,
  name=deblank(dnames(n,:));
  switch (name),
    case 'xi_rho',
      Lp=dsizes(n);
    case 'eta_rho',
      Mp=dsizes(n);
    case 's_rho',
      N=dsizes(n);
      SOLVE3D=1;
    case 'tracer',
      NT=dsizes(n);
  end,
end,
L=Lp-1;
Lm=Lp-2;
M=Mp-1;
Mm=Mp-2;

if (FULL),
  S.lenZeta=Lp*Mp;        Lr=Lp;    Mr=Mp;
  S.lenUbar=Lm*Mp;        Lu=Lm;    Mu=Mp;
  S.lenVbar=Lp*Mm;        Lv=Lp;    Mv=Mm;
  S.lenUvel=Lm*Mp*N;
  S.lenVvel=Lp*Mm*N;
  S.lenTsrc=Lp*Mp*N*NT;
else,
  S.lenZeta=Lm*Mm;        Lr=Lm;    Mr=Mm;
  S.lenUbar=(Lm-1)*Mm;    Lu=Lm-1;  Mu=Mm;
  S.lenVbar=Lm*(Mm-1);    Lv=Lm;    Mv=Mm-1;
  S.lenUvel=(Lm-1)*Mm*N;
  S.lenVvel=Lm*(Mm-1)*N;
  S.lenTrcs=Lm*Mm*N*NT;
end,

%  Get checkpointing pointer array.

ipntr=nc_read(Cname,'ipntr');

%----------------------------------------------------------------------------
%  Extract tangent state variables from WorkD ARPACK vector.
%----------------------------------------------------------------------------

WorkD=nc_read(Cname, 'SworkD');

%  Free surface.

Istr=ipntr(1); Iend=Istr+S.lenZeta;
f(1:Iend-Istr)=WorkD(Istr:Iend-1);
S.tl_zeta=reshape(f,[Lr,Mr]);
clear f;

% Momentum.

if (SOLVE3D),
  Istr=Iend; Iend=Istr+S.lenUvel;
  f(1:Iend-Istr)=WorkD(Istr:Iend-1);
  S.tl_u=reshape(f,[Lu,Mu,N]);
  clear f;

  Istr=Iend; Iend=Istr+S.lenVvel;
  f(1:Iend-Istr)=WorkD(Istr:Iend-1);
  S.tl_v=reshape(f,[Lv,Mv,N]);
  clear f;
else,
  Istr=Iend; Iend=Istr+S.lenUbar;
  f(1:Iend-Istr)=WorkD(Istr:Iend-1);
  S.tl_ubar=reshape(f,[Lu,Mu]);
  clear f;

  Istr=Iend; Iend=Istr+S.lenVbar;
  f(1:Iend-Istr)=WorkD(Istr:Iend-1);
  S.tl_vbar=reshape(f,[Lv,Mv,N]);
  clear f;
end,

% Tracers.

if (SOLVE3D),
  Istr=Iend; Iend=Istr+S.lenTrcs;
  f(1:Iend-Istr)=WorkD(Istr:Iend-1);
  S.tl_t=reshape(f,[Lr,Mr,N,NT]);
  clear f;
end,


%----------------------------------------------------------------------------
%  Extract adjoint state variables from WorkD ARPACK vector.
%----------------------------------------------------------------------------

%  Free surface.

Istr=ipntr(2); Iend=Istr+S.lenZeta;
f(1:Iend-Istr)=WorkD(Istr:Iend-1);
S.ad_zeta=reshape(f,[Lr,Mr]);
clear f;

% Momentum.

if (SOLVE3D),
  Istr=Iend; Iend=Istr+S.lenUvel;
  f(1:Iend-Istr)=WorkD(Istr:Iend-1);
  S.ad_u=reshape(f,[Lu,Mu,N]);
  clear f;

  Istr=Iend; Iend=Istr+S.lenVvel;
  f(1:Iend-Istr)=WorkD(Istr:Iend-1);
  S.ad_v=reshape(f,[Lv,Mv,N]);
  clear f;
else,
  Istr=Iend; Iend=Istr+S.lenUbar;
  f(1:Iend-Istr)=WorkD(Istr:Iend-1);
  S.ad_ubar=reshape(f,[Lu,Mu]);
  clear f;

  Istr=Iend; Iend=Istr+S.lenVbar;
  f(1:Iend-Istr)=WorkD(Istr:Iend-1);
  S.ad_vbar=reshape(f,[Lv,Mv,N]);
  clear f;
end,

% Tracers.

if (SOLVE3D),
  Istr=Iend; Iend=Istr+S.lenTrcs;
  f(1:Iend-Istr)=WorkD(Istr:Iend-1);
  S.ad_t=reshape(f,[Lr,Mr,N,NT]);
  clear f;
end,

%----------------------------------------------------------------------------
%  Extract adjoint state variables from WorkD ARPACK vector.
%----------------------------------------------------------------------------

resid=nc_read(Cname, 'resid');

%  Free surface.

Istr=1; Iend=Istr+S.lenZeta;
f(1:Iend-Istr)=resid(Istr:Iend-1);
S.r_zeta=reshape(f,[Lr,Mr]);
clear f;

% Momentum.

if (SOLVE3D),
  Istr=Iend; Iend=Istr+S.lenUvel;
  f(1:Iend-Istr)=resid(Istr:Iend-1);
  S.r_u=reshape(f,[Lu,Mu,N]);
  clear f;

  Istr=Iend; Iend=Istr+S.lenVvel;
  f(1:Iend-Istr)=resid(Istr:Iend-1);
  S.r_v=reshape(f,[Lv,Mv,N]);
  clear f;
else,
  Istr=Iend; Iend=Istr+S.lenUbar;
  f(1:Iend-Istr)=resid(Istr:Iend-1);
  S.r_ubar=reshape(f,[Lu,Mu]);
  clear f;

  Istr=Iend; Iend=Istr+S.lenVbar;
  f(1:Iend-Istr)=resid(Istr:Iend-1);
  S.r_vbar=reshape(f,[Lv,Mv,N]);
  clear f;
end,

% Tracers.

if (SOLVE3D),
  Istr=Iend; Iend=Istr+S.lenTrcs;
  f(1:Iend-Istr)=resid(Istr:Iend-1);
  S.r_t=reshape(f,[Lr,Mr,N,NT]);
  clear f;
end,

return

