% **** WARNING ****
% This code will only work for a lat/lon model grid

clear all, clf reset, cd /home/dimitri/mitgcm/setups/cap1deg
nx=360; ny=360;
yc=readbin('YC.bin',[nx ny]); xc=readbin('XC.bin',[nx ny]);
yg=readbin('YG.bin',[nx+1 ny+1]); xg=readbin('XG.bin',[nx+1 ny+1]);
[etopo2,lon,lat]=get_etopo2;

% find median in each grid box
bathy=0*xc;
for i=1:nx, mydisp(xc(i,1))
  ix=find(lon>=xg(i,1)&lon<=xg(i+1,1));
  for j=1:ny
    iy=find(lat>=yg(i,j)&lat<=yg(i,j+1));
    tmp=etopo2(ix,iy); in=find(tmp<0);
    if ~isempty(tmp)
      bathy(i,j)=median(tmp(:));
end, end, end

% do not allow any numbers between 0 and 10
topog=bathy; topog(find(topog>-1.5))=0;
topog(find(topog<0&topog>-10))=-10;

% open a few shallow straits, just for fun
topog(1,278)=-10; topog(12,285)=-10; topog(27:28,264)=-10;
topog(43:44,226)=-10; topog(274:275,306)=-10; topog(276:278,316)=-10;
topog(286,311)=-10; topog(285,312)=-10; topog(252,312)=-10;
topog(253,313:314)=-10; topog(255,313)=-10; topog(74,320)=-10;
topog(111,329)=-10; topog(355,258)=-284;

% fill-in lakes
tmp=logical(~topog); tmp=imfill(tmp,'holes');
topog(find(tmp==1))=0; topog(60:61,360)=0; topog(270:282,356:360)=0;

% look at and save
tmp=-topog; tmp(find(~tmp))=1; clf, mypcolor(log10(tmp)'), thincolorbar
writebin('ETOPO2_360x360',topog);