clear all
format compact
more off

nc_file='Permian-bathy-shallow500land300.grd';
Xp=ncread(nc_file,'lon');
Yp=ncread(nc_file,'lat');
Zp=ncread(nc_file,'z');

dlon=2;
dlat=2;
lat_lo=-90;
lat_hi=90-dlat;
lon_hi=360-dlon;
lon_lo=0;

% sNx=45;
% sNy=40;
% nSx=2;
% nSy=1;
% nPx=1;
% nPy=1;
% nxc=sNx*nSx*nPx;
% nyc=sNy*nSy*nPy;

% % Round to nearest "res" metres
% res=5;
% He1=res*round(He1/res);

dz= -[14.8461 17.7881 21.3130 25.5365 30.5969 36.6601 43.9247...
     52.6290 63.0581 75.5539 90.5259 108.4648 129.9585 155.7114...
     186.5677 223.5385 267.8356 320.9107 384.5033 460.6977 551.9909...
     661.3751 792.4352 949.4666];
%dz = -[50 70 100 140 190 240 290 340 390 440 490 540 590 640 690];
zf=[0 cumsum(dz)];

%load VGRID.mat dz zf

% Limit deepest point
Zp(Zp<zf(end))=zf(end);
% Level land
Zp(Zp>0)=0;

xn=lon_lo:dlon:lon_hi;
yn=lat_lo:dlat:lat_hi;

%interpolate data to grid
[X,Y]=meshgrid(Xp,Yp);
[xq,yq]=meshgrid(xn,yn);
Hq=interp2(X,Y,Zp',xq,yq);

Hq=Hq';
prec='real*4';
ieee='b';
fid = fopen('bathy_2deg_nowall.bin', 'w','b');
fwrite(fid, Hq, prec);
fclose(fid);