function [H,x,y]=extract_s2004(lon_lo,lon_hi,lat_lo,lat_hi)
%function [H,x,y]=extract_s2004(lon_lo,lon_hi,lat_lo,lat_hi)
%
% extract topography data from GEBCO 1min data set 
%
  
  ncquiet;
  
  load DATAPATH datapath

% $$$   fname = fullfile(datapath,'gebcosubmap.mat');
% $$$   load(fname)
  fname = fullfile(datapath,'s2004_2min.grd');
  ncload(fname)
  dmap = z;
  clear z
% $$$   nc = netcdf(fname,'r');
% $$$   if isempty(nc);
% $$$     error([fname ' not found'])
% $$$   end
  lon = x;
  lat = y;
  
  if lon(1)==lon(end)-360
    % remove last longitude
    lon(end) = [];
    dmap(:,end) = [];
  end
  nx = length(lon);
  ny = length(lat);
  if size(lon,1) > 1; lon = lon'; end
  if size(lat,1) > 1; lat = lat'; end
  
  if lon_lo<0
    i0=find( lon(1:end-1)<=lon_lo+360 & lon(2:end)>lon_lo+360);
  else
    i0=find( lon(1:end-1)<=lon_lo & lon(2:end)>lon_lo);
  end
  if lon_hi<0
    i1=find( lon<lon_hi+360 & [lon(2:end) 360]>=lon_hi+360)+1;
  else
    i1=find( lon<lon_hi & [lon(2:end) 360]>=lon_hi)+1;
  end
  j0=find( [-90 lat(1:end-1)]<=lat_lo & lat>lat_lo)-1; j0=max(1,j0);
  j1=find( lat<lat_hi & [lat(2:end) 90]>=lat_hi)+1; j1=min(ny,j1);
  if i0<i1
    ii=mod((i0:i1)-1,nx)+1;
  else
    ii=mod([i0:nx 1:i1]-1,nx)+1;
  end
  jj=j0:j1;
  x = lon(ii);
  y = lat(jj);
  H = dmap(jj,ii)';
  
  return