function [H,Z,x,y]=extract_brios(lon_lo,lon_hi,lat_lo,lat_hi)
%function [H,Z,x,y]=extract_brios(lon_lo,lon_hi,lat_lo,lat_hi)
%
% extract topography data from BRIOS data set
%
  
  ncquiet;
  
  load DATAPATH datapath

  fname = fullfile(datapath,'shelfice.mat');
  theResult = load(fname);
  if isempty(theResult);
    error([fname ' not found'])
  end
% $$$   lon = theResult.x+180;
% $$$   lat = theResult.y;
% $$$   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 = -theResult.h(ii,jj);
% $$$   Z = -theResult.z(ii,jj);

  x = theResult.x(:)+180;
  y = theResult.y(:);
  H = -theResult.h;
  Z = -theResult.z;
  
  return