| 1 | dimitri | 1.2 | % **** WARNING **** | 
| 2 |  |  | % This code will only work for a lat/lon model grid | 
| 3 |  |  |  | 
| 4 | heimbach | 1.1 | clear all, clf reset, cd /home/dimitri/mitgcm/setups/cap1deg | 
| 5 |  |  | nx=360; ny=360; | 
| 6 |  |  | yc=readbin('YC.bin',[nx ny]); xc=readbin('XC.bin',[nx ny]); | 
| 7 |  |  | yg=readbin('YG.bin',[nx+1 ny+1]); xg=readbin('XG.bin',[nx+1 ny+1]); | 
| 8 |  |  | [etopo2,lon,lat]=get_etopo2; | 
| 9 |  |  |  | 
| 10 |  |  | % find median in each grid box | 
| 11 |  |  | bathy=0*xc; | 
| 12 |  |  | for i=1:nx, mydisp(xc(i,1)) | 
| 13 |  |  | ix=find(lon>=xg(i,1)&lon<=xg(i+1,1)); | 
| 14 |  |  | for j=1:ny | 
| 15 |  |  | iy=find(lat>=yg(i,j)&lat<=yg(i,j+1)); | 
| 16 |  |  | tmp=etopo2(ix,iy); in=find(tmp<0); | 
| 17 |  |  | if ~isempty(tmp) | 
| 18 |  |  | bathy(i,j)=median(tmp(:)); | 
| 19 |  |  | end, end, end | 
| 20 |  |  |  | 
| 21 |  |  | % do not allow any numbers between 0 and 10 | 
| 22 |  |  | topog=bathy; topog(find(topog>-1.5))=0; | 
| 23 |  |  | topog(find(topog<0&topog>-10))=-10; | 
| 24 |  |  |  | 
| 25 |  |  | % open a few shallow straits, just for fun | 
| 26 |  |  | topog(1,278)=-10; topog(12,285)=-10; topog(27:28,264)=-10; | 
| 27 |  |  | topog(43:44,226)=-10; topog(274:275,306)=-10; topog(276:278,316)=-10; | 
| 28 |  |  | topog(286,311)=-10; topog(285,312)=-10; topog(252,312)=-10; | 
| 29 |  |  | topog(253,313:314)=-10; topog(255,313)=-10; topog(74,320)=-10; | 
| 30 |  |  | topog(111,329)=-10; topog(355,258)=-284; | 
| 31 |  |  |  | 
| 32 |  |  | % fill-in lakes | 
| 33 |  |  | tmp=logical(~topog); tmp=imfill(tmp,'holes'); | 
| 34 |  |  | topog(find(tmp==1))=0; topog(60:61,360)=0; topog(270:282,356:360)=0; | 
| 35 |  |  |  | 
| 36 |  |  | % look at and save | 
| 37 |  |  | tmp=-topog; tmp(find(~tmp))=1; clf, mypcolor(log10(tmp)'), thincolorbar | 
| 38 |  |  | writebin('ETOPO2_360x360',topog); |