function [lon,lat,topo,status] = depth_doit (ncid,londat,latdat,lonlim,latlim);
%
% function [lon,lat,topo,status] = depth_doit (ncid,londat,latdat,lonlim, ...
%                                                                 latlim);
%
% This function extracts the topographic depth bracketting the
% requested locations at the resolution of the gridded data in
% the supplied netCDF file.  If the requested range contains the
% topography file plus some, as much topography as possible is extracted.
%
% ------
% Input:
% ------
%
%    ncid.....Identifier to netCDF topography file.
%    londat...Topography longitude data.
%                londat(1):  number of longitudinal points.
%                londat(2):  starting longitude.
%                londat(3):  longitudinal grid spacing.
%                londat(4):  cyclic flag.
%    latdat...Topography latitude data.
%                latdat(1):  number of latitudinal points.
%                latdat(2):  starting latitude.
%                latdat(3):  latitudinal grid spacing.
%    lonlim...Desired longitude range for topographic depth.
%    latlim...Desired latitude range for topographic depth.
%
% -------
% Output:
% -------
%
%    lon......Extracted longitudes.
%    lat......Extracted latitudes.
%    topo.....Topographic depths.
%    status...Exit status.  [0] no errors.  [<0] error.

%========================================================%
%                                                        %
% WARNING:   Does not yet take advantage of cyclic flag. %
%                                                        %
%========================================================%

%-------------------------------------------------------------------------------
% Initialize output.
%-------------------------------------------------------------------------------

topo = [];

%-------------------------------------------------------------------------------
% Extract position data.
%-------------------------------------------------------------------------------

nlon     = round(londat(1));
lonst    = londat(2);
dlon     = londat(3);
periodic = londat(4);

nlat  = round(latdat(1));
latst = latdat(2);
dlat  = latdat(3);

%-------------------------------------------------------------------------------
% Ensure all on same branch.  Inner remainder ensures [-360, 360],
% outer remainder ensures [0, 360].
%-------------------------------------------------------------------------------

%lonstwk = rem(360+rem(lonst,360),360);
lonstwk = lonst;

%lonwk = rem(360+rem(lonlim,360),360);
lonwk = lonlim;

%-------------------------------------------------------------------------------
% Calculate desired indices.
%-------------------------------------------------------------------------------

%---------------------------
%--- Longitudinal index. ---
%---------------------------

istL = [];
indR = [];

%x   = rem(360+rem((lonwk(1)-lonstwk),360),360)/dlon;
x = (lonwk(1)-lonstwk)/dlon;
ist = floor (x);
if (ist==(nlon-1)),
   ist = ist-1;
end;
if (ist<0)
   istL = ist;
   ist  = 0;
end;
%ist = max([ist 0]);

%x   = rem(360+rem((lonwk(2)-lonstwk),360),360)/dlon;
x = (lonwk(2)-lonstwk)/dlon;
ind = ceil (x);
if (ind==(nlon-1)),
   ind = ind-1;
end;
if (ind>(nlon-1))
   indR = ind;
   ind  = nlon - 1;
end;
%ind = min([ind (nlon-1)]);

if ( (ist>(nlon-1)) |  (ind<0) ),
   disp(['***Error:  DEPTH - requested range not in file.']);
   disp(['           requested longitudes:  [',num2str(lonwk(1)),', ', ...
                     num2str(lonwk(2)),']']);
   disp(['           longitude range:  [',num2str(lonst), ...
         ', ',num2str(lonst+(nlon-1)*dlon),']']);
   status = -1;
   return;
end;

nx = ind - ist + 1;

lon = lonst + (ist:ind).*dlon;

lon = lon + 360*round((lonlim(1)-lon(1))/360);

%--------------------------
%--- Latitudinal index. ---
%--------------------------

y   = (latlim(1)-latst)/dlat;
jst = floor (y);
if (jst==(nlat-1)),
   jst = jst-1;
end;
jst = max([jst 0]);

y   = (latlim(2)-latst)/dlat;
jnd = ceil (y);
if (jnd==(nlat-1)),
   jnd = jnd-1;
end;
jnd = min([jnd (nlat-1)]);

if ( (jst>(nlat-1)) | (jnd<0) ),
   disp(['***Error:  DEPTH - requested range not in file.',setstr(7)]);
   disp(['           requested latitudes:  [',num2str(latlim(1)),', ', ...
              num2str(latlim(1)),']']);
   disp(['           latitude range:  [',num2str(latst), ...
         ', ',num2str(latst+(nlat-1)*dlat),']']);
   status = -1;
   return;
end;

ny = jnd - jst + 1;

lat = latst + (jst:jnd).*dlat;

%-------------------------------------------------------------------------------
% Read topographic square.
%-------------------------------------------------------------------------------

[topo, status] = get_var (ncid,'topo',[jst ist],[ny nx]);

if (status < 0),
   disp(['***Error:  DEPTH - unable to read variable ', ...
                      setstr(34),'topo',setstr(34),'.']);
   return
end;

%----------------------------
%--- Check left wrapping. ---
%----------------------------

%if (periodic & ~isempty(istL))

topo = topo';
