Index: /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpArcticdem.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpArcticdem.m	(revision 27547)
+++ /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpArcticdem.m	(revision 27547)
@@ -0,0 +1,68 @@
+function sout = interpArcticdem(X,Y),
+
+switch oshostname(),
+	case {'ronne'}
+		path='/home/ModelData/Greenland/ArcticDemMosaic/arcticdem_mosaic_100m_v3.0.tif';
+	case {'totten'}
+		path='/totten_1/ModelData/Greenland/ArcticDemMosaic/arcticdem_mosaic_100m_v3.0.tif';
+	otherwise
+		error('machine not supported yet');
+end
+
+usemap = 0;
+if license('test','map_toolbox')==0,
+	disp('WARNING: map toolbox not installed, trying house code');
+	usemap = 0;
+elseif license('checkout','map_toolbox')==0
+	disp('WARNING: map toolbox not available (checkout failed), trying house code');
+	usemap = 0;
+end
+
+if usemap,
+	[data,R] = geotiffread(path);
+	data=double(flipud(data));
+	xdata=R.XLimWorld(1):R.DeltaX:R.XLimWorld(2); xdata=xdata(:);
+	xdata =(xdata(1:end-1)+xdata(2:end))/2;
+	ydata=R.YLimWorld(2):R.DeltaY:R.YLimWorld(1); ydata=flipud(ydata(:));
+	ydata =(ydata(1:end-1)+ydata(2:end))/2;
+else
+
+	%Get image info
+	Tinfo = imfinfo(path);
+	N     = Tinfo.Width;
+	M     = Tinfo.Height;
+	dx    = Tinfo.ModelPixelScaleTag(1);
+	dy    = Tinfo.ModelPixelScaleTag(2);
+	minx  = Tinfo.ModelTiepointTag(4);
+	maxy  = Tinfo.ModelTiepointTag(5);
+
+	%Generate vectors
+	xdata = minx + dx/2 + ((0:N-1).*dx);
+	ydata = maxy - dy/2 - ((M  -1:-1:0).*dy);
+	ydata = fliplr(ydata);
+
+	%Get pixels we are interested in
+	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);
+
+	if 0,
+		ymin=min(Y(:)); ymax=max(Y(:));
+		posy=find(ydata<=ymax);
+		id1y=max(1,find(ydata>=ymin,1)-offset);
+		id2y=min(numel(ydata),posy(end)+offset);
+	else
+		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);
+	end
+
+	data  = double(imread(path,'PixelRegion',{[id1y,id2y],[id1x,id2x]}));
+	xdata=xdata(id1x:id2x);
+	ydata=ydata(id1y:id2y);
+end
+
+sout = InterpFromGrid(xdata,ydata,data,X,Y);
Index: /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBamber2009.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBamber2009.m	(revision 27547)
+++ /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBamber2009.m	(revision 27547)
@@ -0,0 +1,21 @@
+function demout = interpBamber2009(X, Y)
+%INTERPBAMBER2009 - interpolate surface dem of Bamber 2009
+%
+%   Surface dem Nominal year 2004 (WGS84, no firn correction)
+%
+%   Usage:
+%      demout = interpBamber2009(X, Y)
+
+switch oshostname(),
+	case {'totten'}
+		bamber2009path ='/totten_1/ModelData/Antarctica/Bamber2009DEM/krigged_dem_nsidc.mat';
+	otherwise
+		error('machine not supported yet');
+end
+
+%Convert to Bamber's projections
+%disp('   -- Bamber2009: loading dem'); 
+load(bamber2009path);
+
+disp('   -- Bamber2009: interpolating dem (WGS84)');
+demout = InterpFromGrid(x, y, surfacedem, X, Y);
Index: /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpGrIMP.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpGrIMP.m	(revision 27547)
+++ /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpGrIMP.m	(revision 27547)
@@ -0,0 +1,66 @@
+function sout = interpGimpdem(X,Y),
+
+switch oshostname(),
+	case {'totten'}
+		howatpath='/totten_1/ModelData/Greenland/GrIMP/GrIMP_100m.tif';
+	otherwise
+		error('machine not supported yet');
+end
+
+usemap = 0;
+if license('test','map_toolbox')==0,
+	disp('WARNING: map toolbox not installed, trying house code');
+	usemap = 0;
+elseif license('checkout','map_toolbox')==0
+	disp('WARNING: map toolbox not available (checkout failed), trying house code');
+	usemap = 0;
+end
+
+if usemap,
+	[data,R] = geotiffread(howatpath);
+	data=double(flipud(data));
+	xdata=R.XLimWorld(1):R.DeltaX:R.XLimWorld(2); xdata=xdata(:);
+	xdata =(xdata(1:end-1)+xdata(2:end))/2;
+	ydata=R.YLimWorld(2):R.DeltaY:R.YLimWorld(1); ydata=flipud(ydata(:));
+	ydata =(ydata(1:end-1)+ydata(2:end))/2;
+else
+
+	%Get image info
+	Tinfo = imfinfo(howatpath);
+	N     = Tinfo.Width;
+	M     = Tinfo.Height;
+	dx    = Tinfo.ModelPixelScaleTag(1);
+	dy    = Tinfo.ModelPixelScaleTag(2);
+	minx  = Tinfo.ModelTiepointTag(4);
+	maxy  = Tinfo.ModelTiepointTag(5);
+
+	%Generate vectors
+	xdata = minx + dx/2 + ((0:N-1).*dx);
+	ydata = maxy - dy/2 - ((M  -1:-1:0).*dy);
+	ydata = fliplr(ydata);
+
+	%Get pixels we are interested in
+	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);
+
+	if 0,
+		ymin=min(Y(:)); ymax=max(Y(:));
+		posy=find(ydata<=ymax);
+		id1y=max(1,find(ydata>=ymin,1)-offset);
+		id2y=min(numel(ydata),posy(end)+offset);
+	else
+		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);
+	end
+
+	data  = double(imread(howatpath,'PixelRegion',{[id1y,id2y],[id1x,id2x]}));
+	xdata=xdata(id1x:id2x);
+	ydata=ydata(id1y:id2y);
+end
+
+sout = InterpFromGrid(xdata,ydata,data,X,Y);