Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/computeGrad.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/computeGrad.m	(revision 27693)
+++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/computeGrad.m	(revision 27694)
@@ -1,4 +1,4 @@
 function [gradx, grady]=computeGrad(index,x,y,field)
-%COMPUTEHESSIAN - compute the gradient from a field
+%COMPUTEGRAD - compute the gradient from a field
 
 %some variables
Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/extractTransientSolutions.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/extractTransientSolutions.m	(revision 27693)
+++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/extractTransientSolutions.m	(revision 27694)
@@ -9,5 +9,7 @@
 transientSolutions.thickness = cell2mat({md.results.TransientSolution(:).Thickness});
 transientSolutions.SigmaVM = cell2mat({md.results.TransientSolution(:).SigmaVM});
-transientSolutions.smb = cell2mat({md.results.TransientSolution(:).SmbMassBalance});
+if (isfield(md.results.TransientSolution, 'SmbMassBalance'))
+	transientSolutions.smb = cell2mat({md.results.TransientSolution(:).SmbMassBalance});
+end
 transientSolutions.ice_levelset = cell2mat({md.results.TransientSolution(:).MaskIceLevelset});
 transientSolutions.calvingRate = cell2mat({md.results.TransientSolution(:).CalvingCalvingrate});
Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverDomain.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverDomain.m	(revision 27693)
+++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverDomain.m	(revision 27694)
@@ -8,11 +8,13 @@
 	weights = ones(size(data));
 	if nargin<3
-		masked = [];
+		masked = logical(zeros(size(data)));
 	end
 end
 
+masked = masked | isnan(data) | isnan(weights);
 % Set the area with masked=1 to nan
 data(masked) = nan;
 weights(masked) =nan;
+
 
 % get the mesh
@@ -28,5 +30,5 @@
 eleAreas = 1/3*eleAreas.*(weights(elements(:,1),:)+weights(elements(:,2),:)+weights(elements(:,3),:));
 
-intData = sum(eleData(:),'omitnan');
-areas = sum(eleAreas(:),'omitnan');
-meanData = intData / areas;
+intData = sum(eleData, 1, 'omitnan');
+areas = sum(eleAreas, 1, 'omitnan');
+meanData = intData ./ areas;
Index: /issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m	(revision 27693)
+++ /issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m	(revision 27694)
@@ -27,4 +27,8 @@
 elseif strcmp(glacier, 'Rink')
 	foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Rink_2008_2022/';
+elseif strcmp(glacier, 'Upernavik')
+	foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Upernavik_2008_2022/';
+elseif strcmp(glacier, 'Helheim')
+	foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Helheim_2008_2023/';
 else
 	error(['The velocity data for ', glacier, ' is not available, please download from NSIDC first.']);
Index: /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpMonthlyIceMaskGreene.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpMonthlyIceMaskGreene.m	(revision 27694)
+++ /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpMonthlyIceMaskGreene.m	(revision 27694)
@@ -0,0 +1,68 @@
+function icemask = interpMonthlyIceMaskGreene(X, Y, time, includedRockMask, ncdata)
+%INTERPMONTHLYICEMASKGREENE - interpolate monthly reconstructed ice masks onto X and Y, within the given time period
+%
+%	 Usage:
+%		distance = interpMonthlyIceMaskGreene(md.mesh.x, md.mesh.y, [md.timestepping.start_time, md.timestepping.final_time]);
+%
+%
+%   - X, Y: coordinates of the mesh or grid
+%	 - time: the starting and end point of the time series
+%   - optional 4th input argument: a flag to cover rock mask by ice and merge it into the ice mask, 
+%											so that there is no hole in the interior part of the domain, 
+%											and the mask will only represent the levelset of the ice front 
+%   - optional 5th input argument: path to the data file, by default it is the path on totten
+%
+%   Author: Cheng Gong
+%   Last modified: 2023-04-19
+
+% set icemask=-1 for the region with rocks
+if nargin < 4
+	includedRockMask = 1;
+end
+if nargin < 5
+	ncdata = '/totten_1/ModelData/Greenland/IceFrontsGreene/greenland_ice_masks_1972-2022_v1.nc';
+end
+
+x = ncread(ncdata, 'x');
+y = ncread(ncdata, 'y');
+d = ncread(ncdata, 'time');
+% convert t to decyear
+t = date2decyear(datenum(datetime('1900-01-01')+days(d)));
+
+offset=2;
+
+% get x-index covers the domain
+xmin=min(X(:)); xmax=max(X(:));
+idx = sort(find((x>=xmin) & (x<=xmax)));
+idx_min = max(idx(1)-offset, 1);
+idx_max = min(idx(end)+offset, length(x));
+x = x(idx_min:idx_max);
+
+% get y-index covers the domain
+ymin=min(Y(:)); ymax=max(Y(:));
+idy = sort(find((y>=ymin) & (y<=ymax)));
+idy_min = max(idy(1)-offset, 1);
+idy_max = min(idy(end)+offset, length(y));
+y = y(idy_min:idy_max);
+
+% get time index
+idt_min = max(find(t<=time(1), 1, 'last'), 1);
+idt_max = min(find(t>=time(end), 1, 'first'), length(t));
+t = t(idt_min:idt_max);
+
+% load icemask and rockmask from netCDF
+ice = ncread(ncdata, 'ice', [idx_min, idy_min, idt_min], [idx_max-idx_min+1, idy_max-idy_min+1, idt_max-idt_min+1], [1,1,1]);
+rock = ncread(ncdata, 'rock', [idx_min, idy_min], [idx_max-idx_min+1, idy_max-idy_min+1], [1,1]);
+
+% merge ice and rock
+if includedRockMask
+	iceall = ice + rock;
+else 
+	iceall = ice;
+end
+% Convert to ice_levelset values
+icemask = zeros(numel(X)+1, numel(t));
+icemask(end,:) = t;
+for i = 1:numel(t)
+	icemask(1:end-1, i) = InterpFromGrid(x, y, double(1-2*iceall(:,:,i)'), X, Y,'nearest');
+end
Index: /issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m	(revision 27693)
+++ /issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m	(revision 27694)
@@ -16,5 +16,5 @@
 %
 %   Author: Cheng Gong
-%   Date: 2021-12-06
+%   Date: 2023-04-10
 
 N = length(velList);
@@ -46,3 +46,7 @@
     h=colorbar;
     title(h,'m/a')
+	 % add a vertical line to indicate initial ice front position
+	 id = max(find(~isnan(cumsum(vel_flowline{p}(:,1)))));
+	 icefront = flowline.Xmain(id);
+	 plot([icefront, icefront], [2007,2020], '-r')
 end