Index: /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/interpBedmachine.m
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/interpBedmachine.m	(revision 25589)
+++ /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/interpBedmachine.m	(revision 25589)
@@ -0,0 +1,44 @@
+function output = interpBedmachine(X,Y,string,continent),
+
+if strcmpi(continent,'antarctica'),
+	ncdate='2014-09-23';
+	morlighemnc=['/Users/larour/ModelData2/MMAntarctica2013/AntarcticaMCdataset-' ncdate '.nc'];
+else
+	%ncdate='2015-10-05';
+	ncdate='2017-04-04';
+	morlighemnc=['/Users/larour/ModelData2/MMGreenland2013/MCdataset-' ncdate '.nc'];
+end
+
+disp(['   -- BedMachine version: ' ncdate]);
+xdata = double(ncread(morlighemnc,'x'));
+ydata = double(ncread(morlighemnc,'y'));
+
+offset=2;
+
+xmin=min(X(:)); xmax=max(X(:));
+posx=find(xdata<=xmax);
+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);
+id1y=max(1,find(ydata<=ymax,1)-offset);
+id2y=min(numel(ydata),posy(end)+offset);
+
+disp(sprintf('   -- BedMachine %s : loading %s',continent,string));
+data  = double(ncread(morlighemnc,string,[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))';
+xdata=xdata(id1x:id2x);
+ydata=ydata(id1y:id2y);
+data(find(data==-999999))=NaN;
+
+disp(sprintf('   -- BedMachine %s : interpolating %s',continent,string));
+if strcmp(string,'mask') | strcmp(string,'source'),
+	%Need nearest neighbor to avoid interpolation between 0 and 2
+	if strcmpi(continent,'greenland'), ydata=flipud(ydata); data=flipud(data); end
+	output = InterpFromGridToMesh(xdata,ydata,data,double(X),double(Y),0);
+else
+	ydata=flipud(ydata); data=flipud(data); 
+	output = InterpFromGridToMesh(xdata,ydata,data,double(X),double(Y),0);
+end
+
+output(find(output==-999999))=NaN;
Index: /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/interpBedmap2.m
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/interpBedmap2.m	(revision 25589)
+++ /issm/branches/trunk-larour-SLPS2020/src/m/contrib/larour/interpBedmap2.m	(revision 25589)
@@ -0,0 +1,47 @@
+function [dataout] = interpBedmap2(X,Y,string),
+%INTERPBEDMAP2 - interpolate bedmap2 data
+%
+%   Available data:
+%      1. bed                          is bed height
+%      2. surface                      is surface height
+%      3. thickness                    is ice thickness
+%      4. icemask_grounded_and_shelves is a mask file showing the grounding line and the extent of the floating ice shelves
+%      5. rockmask                     is a mask file showing rock outcrops
+%      6. lakemask_vostok              is a mask file showing the extent of the lake cavity of Lake Vostok
+%      7. bed_uncertainty              is the bed uncertainty grid shown in figure 12 of the manuscript
+%      8. thickness_uncertainty_5km    is the thickness uncertainty grid shown in figure 11 of the manuscript
+%      9. data_coverage                is a binary grid showing the dis tribution of ice thickness data used in the grid of ice thickness
+%
+%   Usage:
+%      [dataout] = interpBedmap2(X,Y,string)
+
+%reqad data
+%path = '/u/astrid-r1b/ModelData/BedMap2/bedmap2_bin/';
+%path = '~/issm-jpl/proj-group/ModelData/BedMap2/bedmap2_bin/';
+path = '/Users/larour/ModelData/BedMap2/bedmap2_bin/';
+fid=fopen([path '/bedmap2_' string '.flt'],'r','l');
+data=fread(fid,[6667,6667],'float32');
+fclose(fid);
+
+% define grid
+if strcmp(string,'thickness_uncertainty_5km'),
+	ncols    =1361;
+	nrows    =1361;
+	xll      =-3401000;
+	yll      =-3402000;
+	gridsize =5000;
+else
+	ncols    =6667;
+	nrows    =6667;
+	xll      =-3333000;
+	yll      =-3333000;
+	gridsize =1000;
+end
+x_m=xll+(0:1:ncols-1)'*gridsize;
+y_m=yll+(0:1:nrows-1)'*gridsize;
+
+%Change default to NaN
+data(find(data==-9999))=NaN;
+
+%Interpolate
+dataout = InterpFromGridToMesh(x_m,y_m,flipud(data'),double(X),double(Y),0);