Index: /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/SESAWslr.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/SESAWslr.m	(revision 25442)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/SESAWslr.m	(revision 25443)
@@ -1,3 +1,3 @@
-function [eust,rsl,vlm] = SESAWslr(index,lat,long,greens,para) 
+function [eust,rsl,vlm,geoid] = SESAWslr(index,lat,long,greens,para) 
 %SESAWslr :: computes GRD slr due to applied surface loads based on SESAW method. 
 % 
@@ -5,5 +5,10 @@
 %
 %   Usage:
-%      [eus,rsl]=greensfunctions(md.mesh.elements,md.mesh.lat,md.mesh.long,greens,para) 
+%      [eus,rsl,vlm,geoid]=greensfunctions(md.mesh.elements,md.mesh.lat,md.mesh.long,greens,para) 
+% 
+%      eus = eustatic aka barystatic sea level
+%      rsl = relative sea level 
+%      vlm = vertical land motion 
+%      geoid = change in geoid height 
 %      
 %      greens.Grigid = Gravitational Green's function for the rigid planet. 
@@ -105,7 +110,8 @@
 end 
 
-% compute bedrock motions (VLM: vertical land motion) 
+% compute bedrock motion aka VLM (vertical land motion) and geoid height change 
+rsl_no_mass_conservation = term1+term2; % eus-term3-term4 is excluded. See Tamisiea, 2011. 
 if strcmpi(para.solidearth,'rigid')
-	vlm = 0.0*rsl;   
+	vlm = 0.0*rsl;
 elseif strcmpi(para.solidearth,'elastic') 
 	% reset Green's function for VLM, i.e. h_l
@@ -115,4 +121,5 @@
 	vlm = term1+term2; 
 end
+% now add VLM. 
+geoid = rsl_no_mass_conservation+vlm; 
 
-