function h=oct_add_topo(grdname,toponame)
% OCT_ADD_TOPO  Read grid and topology using octave-netcdf and interpolate
ncid = netcdf.open(grdname,'NC_NOWRITE');
lon = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'lon_rho'));
lat = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'lat_rho'));
pm  = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'pm'));
pn  = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'pn'));
netcdf.close(ncid);

% Get CROCO averaged resolution

dx=mean(mean(1./pm));
dy=mean(mean(1./pn));
dx_croco=mean([dx dy]);
disp(['   CROCO resolution : ',num2str(dx_croco/1000,3),' km'])

dl=max([1 2*(dx_croco/(60*1852))]);
lonmin=min(min(lon))-dl;
lonmax=max(max(lon))+dl;
latmin=min(min(lat))-dl;
latmax=max(max(lat))+dl;

% open the topo file
ncid = netcdf.open(toponame,'NC_NOWRITE');
tl = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'lon'));
tt = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'lat'));
% Try standard variable names for topo
try
  topo_var = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'topo'));   % Octave added transpose
catch
  % fallback: try 'z' or 'bath'
  try topo_var = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'z')); catch; end
  try topo_var = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'bath')); catch; end
end
netcdf.close(ncid);
% adapt longitude wrapping
j=find(tt>=latmin & tt<=latmax);
i1=find(tl-360>=lonmin & tl-360<=lonmax);
i2=find(tl>=lonmin & tl<=lonmax);
i3=find(tl+360>=lonmin & tl+360<=lonmax);
x=cat(1,tl(i1)-360,tl(i2),tl(i3)+360);
y=tt(j);

% Read data (assemble wrapped sections)
if ~isempty(i2)
%%  topo=-topo_var(j,i2);
  topo=-topo_var(i2,j);
else
  topo=[];
end
if ~isempty(i1)
  topo=cat(2,-topo_var(j,i1),topo);
end
if ~isempty(i3)
  topo=cat(2,topo,-topo_var(j,i3));
end

% Degrade topography if needed (same as original)
R=6367442.76;
deg2rad=pi/180;
dg=mean(x(2:end)-x(1:end-1));
dphi=y(2:end)-y(1:end-1);
dy=R*deg2rad*dphi;
dx=R*deg2rad*dg*cos(deg2rad*y);
dx_topo=mean([dx ;dy]);
disp(['   Topography data resolution : ',num2str(dx_topo/1000,3),' km'])

n=0;
while dx_croco>(dx_topo)
  n=n+1;
  x=0.5*(x(2:end)+x(1:end-1));
  x=x(1:2:end);
  y=0.5*(y(2:end)+y(1:end-1));
  y=y(1:2:end);
  topo=0.25*(topo(2:end,1:end-1)  +topo(2:end,2:end)+... 
             topo(1:end-1,1:end-1)+topo(1:end-1,2:end));
  topo=topo(1:2:end,1:2:end);
  dg=mean(x(2:end)-x(1:end-1));
  dphi=y(2:end)-y(1:end-1);
  dy=R*deg2rad*dphi;
  dx=R*deg2rad*dg*cos(deg2rad*y);
  dx_topo=mean([dx ;dy]);
end

disp(['   Topography resolution halved ',num2str(n),' times'])
disp(['   New topography resolution : ',num2str(dx_topo/1000,3),' km'])

% interpolate the topo
keyboard
h=interp2(x,y,topo',lon,lat,'cubic');
return
