function output = interpBedmachineGreenland(X,Y,string,ncdate),

if nargin<3, string = 'bed'; end
if nargin<4,
	%ncdate='2015-04-27'; %BedMachine v2
	ncdate='2017-09-25'; %BedMachine v3
	ncdate='2020-04-14';
end

basename = 'BedMachineGreenland';

%List of common paths to try
paths = {...
	['/u/astrid-r1b/ModelData/ModelData/MCdataset-' ncdate '.nc'],...
	['/home/ModelData/Greenland/BedMachine/' basename '-' ncdate '.nc'],...
	['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'],...
	['./' basename '-' ncdate '.nc'],...
	};

found = 0;
for i=1:numel(paths)
	if exist(paths{i},'file')
		ncfile = paths{i};
		found = 1;
		break;
	end
end

if ~found
	error(['Could not find ' basename '-' ncdate '.nc, you can add the path to the list']);
end


disp(['   -- BedMachine Greenland version: ' ncdate]);
xdata = double(ncread(ncfile,'x'));
ydata = double(ncread(ncfile,'y'));

offset=2;

xmin=min(X(:)); xmax=max(X(:));
posx=find(xdata<=xmax);
if isempty(posx), posx=numel(xdata); end
id1x=max(1,find(xdata>=xmin,1)-offset);
id2x=min(numel(xdata),posx(end)+offset);

ymin=min(Y(:)); ymax=max(Y(:));
posy=find(ydata>=ymin);
if isempty(posy), posy=numel(ydata); end
id1y=max(1,find(ydata<=ymax,1)-offset);
id2y=min(numel(ydata),posy(end)+offset);

disp(['   -- BedMachine Greenland: loading ' string]);
data  = double(ncread(ncfile,string,[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))';
xdata=xdata(id1x:id2x);
ydata=ydata(id1y:id2y);
data(find(data==-9999))=NaN;

disp(['   -- BedMachine Greenland: interpolating ' string]);
if strcmp(string,'mask') | strcmp(string,'source'),
	%Need nearest neighbor to avoid interpolation between 0 and 2
	output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest');
else
	output = InterpFromGrid(xdata,ydata,data,double(X),double(Y));
end

%TEST https://www.mathworks.com/matlabcentral/fileexchange/10772-fast-2-dimensional-interpolation
