Index: /issm/trunk-jpl/src/m/contrib/dlcheng/MeshMergedOutlineToShp.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/dlcheng/MeshMergedOutlineToShp.m	(revision 25733)
+++ /issm/trunk-jpl/src/m/contrib/dlcheng/MeshMergedOutlineToShp.m	(revision 25733)
@@ -0,0 +1,21 @@
+function MeshMergedOutlineToShp(mds,shapefilename)
+%MeshMergedOutlineToShp - export multiple mesh outlines to shp file
+%
+%   Usage:
+%      MeshMergedOutlineToShp({md1,md2},'Greenland.shp');
+
+	contours= struct([]);
+    count=1;
+    for i=1:length(mds)
+        md = mds{i};
+        for j=1:md.mesh.numberofvertices
+            if md.mesh.vertexonboundary(j)
+                contours(count).x = md.mesh.x(j);
+                contours(count).y = md.mesh.y(j);
+                contours(count).id = count;
+                contours(count).Geometry = 'Point';
+                count = count + 1;
+            end
+        end
+	end
+	shpwrite(contours,shapefilename);
Index: /issm/trunk-jpl/src/m/contrib/dlcheng/MeshOutlineToShp.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/dlcheng/MeshOutlineToShp.m	(revision 25733)
+++ /issm/trunk-jpl/src/m/contrib/dlcheng/MeshOutlineToShp.m	(revision 25733)
@@ -0,0 +1,16 @@
+function MeshOutlineToShp(md,shapefilename)
+%MESHOUTLINETOSHP - export mesh outline to shp file
+%
+%   Usage:
+%      MeshOutlineToShp(md,'Greenland.shp');
+
+	contours= struct([]);
+	for i=1:md.mesh.numberofvertices
+        if md.mesh.vertexonboundary(i)
+            contours(i).x = md.mesh.x(i);
+            contours(i).y = md.mesh.y(i);
+            contours(i).id = i;
+            contours(i).Geometry = 'Point';
+        end
+	end
+	shpwrite(contours,shapefilename);
Index: /issm/trunk-jpl/src/m/contrib/dlcheng/VxVyToVv.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/dlcheng/VxVyToVv.m	(revision 25733)
+++ /issm/trunk-jpl/src/m/contrib/dlcheng/VxVyToVv.m	(revision 25733)
@@ -0,0 +1,37 @@
+name='CWGreenland.'; %new Central West Greenland project to do AD from Store to Upernavik.
+org=organizer('repository','./Models','prefix',name,'steps',steps);
+NSIDCName='./Data/Greenland';
+WestCoastName='./Data/Rink/';
+modeldatapath='/Users/larour/ModelData/';
+modeldata2path='/Users/larour/ModelData2/';
+MARName=[modeldata2path 'MARv3.11/ERA/CW_Greenland/'];
+qgispath='/Users/larour/issm-jpl/proj-group/qgis/';
+greenlandshppath=[qgispath '/GreenlandStateEstimate/'];
+GIMP='/Users/larour/ModelData/HowatDEMGreenland2012/HowatDEMGreenland2012_90m_smooth.mat';
+
+%Path to shapefiles: 
+shppath=[qgispath '/Greenland/30km/'];
+shppathslr=[qgispath '/Slr/'];
+
+if perform(org,'GenerateLandsatVV'), % {{{
+    rootname = [modeldatapath '/VelHowat'];
+    infos = dir([rootname '/*/*/*_vx_*.tif']);
+    for i=1:length(infos),
+        filename = infos(i).name; %Example: OPT_W70.55N_1986-05_vx_v02.1.tif
+        path = [infos(i).folder '/' infos(i).name];
+
+        vx_path = path;
+        vy_path = strrep(vx_path, '_vx_', '_vy_');
+        vv_path = strrep(vx_path, '_vx_', '_vv_');
+
+        [vx,Rx] = readgeoraster(vx_path);
+        [vy,Ry] = readgeoraster(vy_path);
+        vx(vx==-99999)=nan;
+        vy(vy==-99999)=nan;
+
+        vel = sqrt(vx.^2+vy.^2);
+        disp(vv_path);
+        geotiffwrite(vv_path, vel, Rx, 'CoordRefSysCode', 3413);
+    end
+end % }}}
+
Index: /issm/trunk-jpl/src/m/contrib/dlcheng/interpBedmachineGreenland.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/dlcheng/interpBedmachineGreenland.m	(revision 25733)
+++ /issm/trunk-jpl/src/m/contrib/dlcheng/interpBedmachineGreenland.m	(revision 25733)
@@ -0,0 +1,96 @@
+function output = interpBedmachineGreenland(X,Y,string,ncdate),
+
+if nargin<3, string = 'bed'; end
+if nargin<4,
+	%ncdate='2013-05-21';
+	%ncdate='2013-06-27';
+	%ncdate='2013-07-18';
+	%ncdate='2013-11-15';
+	%ncdate='2013-12-03';
+	%ncdate='2014-02-26';
+	%ncdate='2014-03-24';
+	%ncdate='2014-07-31';
+	%ncdate='2014-11-14';
+	%ncdate='2015-03-03';
+	%ncdate='2015-03-10';
+	%ncdate='2015-03-26';
+	%ncdate='2015-03-30';
+	%ncdate='2015-04-27'; %BedMachine v2
+	%ncdate='2015-07-30';
+	%ncdate='2015-10-02';
+	%ncdate='2016-03-21';
+	%ncdate='2016-05-12';
+	ncdate='2016-07-06';
+	ncdate='2016-08-04';
+	ncdate='2016-10-26';
+	ncdate='2016-11-23';
+	ncdate='2016-12-21';
+	ncdate='2017-01-19';
+	ncdate='2017-03-28';
+	ncdate='2017-05-10';
+	ncdate='2017-07-21';
+	ncdate='2017-09-25'; %BedMachine v3
+	ncdate='2018-06-01';
+	ncdate='2019-04-18';
+	ncdate='2020-04-14';
+end
+
+if exist('datetime','file') 
+	date1 = sscanf(ncdate,'%d-%d-%d');
+	date2 = datetime(date1(1),date1(2),date1(3));
+	if date2<datetime(2016,10,24),
+		basename = 'MCdataset'; 
+	else
+		basename = 'BedMachineGreenland';
+	end
+else
+  basename = 'BedMachineGreenland';
+end
+
+switch oshostname(),
+	case {'murdo','thwaites','astrid'}
+		morlighem2013nc=['/u/astrid-r1b/ModelData/ModelData/MCdataset-' ncdate '.nc']';
+	case {'ronne'}
+		morlighem2013nc=['/home/ModelData/Greenland/BedMachine/' basename '-' ncdate '.nc'];
+	otherwise
+		morlighem2013nc=['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'];
+		disp('machine not supported yet, using default path /Users/larour/ModelData/BedMachine/');
+end
+
+disp(['   -- BedMachine Greenland version: ' ncdate]);
+xdata = double(ncread(morlighem2013nc,'x'));
+ydata = double(ncread(morlighem2013nc,'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(morlighem2013nc,string,[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))';
+xdata=xdata(id1x:id2x);
+ydata=ydata(id1y:id2y);
+data(find(data==-9999))=NaN;
+
+ydata=flipud(ydata);
+data=flipud(data);
+    
+disp(['   -- BedMachine Greenland: interpolating ' string]);
+if strcmp(string,'mask') | strcmp(string,'source'),
+	%Need nearest neighbor to avoid interpolation between 0 and 2
+	output = InterpFromGridToMesh(xdata,ydata,data,double(X),double(Y),'nearest');
+else
+    
+	output = InterpFromGridToMesh(xdata,ydata,data,double(X),double(Y),NaN);
+end
+
+%TEST https://www.mathworks.com/matlabcentral/fileexchange/10772-fast-2-dimensional-interpolation
Index: /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m	(revision 25732)
+++ /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m	(revision 25733)
@@ -54,5 +54,6 @@
 		morlighem2013nc=['/home/ModelData/Greenland/BedMachine/' basename '-' ncdate '.nc'];
 	otherwise
-		error('machine not supported yet');
+		morlighem2013nc=['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'];
+		disp('machine not supported yet, using default path /Users/larour/ModelData/BedMachine/');
 end
 
@@ -81,10 +82,14 @@
 data(find(data==-9999))=NaN;
 
+ydata=flipud(ydata);
+data=flipud(data);
+    
 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');
+	output = InterpFromGridToMesh(xdata,ydata,data,double(X),double(Y),'nearest');
 else
-	output = InterpFromGrid(xdata,ydata,data,double(X),double(Y));
+    
+	output = InterpFromGridToMesh(xdata,ydata,data,double(X),double(Y),NaN);
 end