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 25439)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/SESAWslr.m	(revision 25440)
@@ -1,3 +1,3 @@
-function [eust,rsl] = SESAWslr(index,lat,long,greens,para) 
+function [eust,rsl,vlm] = SESAWslr(index,lat,long,greens,para) 
 %SESAWslr :: computes GRD slr due to applied surface loads based on SESAW method. 
 % 
@@ -105,3 +105,14 @@
 end 
 
+% compute bedrock motions (VLM: vertical land motion) 
+if strcmpi(para.solidearth,'rigid')
+	vlm = 0.0*rsl;   
+elseif strcmpi(para.solidearth,'elastic') 
+	% reset Green's function for VLM, i.e. h_l
+	Galpha = greens.Uelast;		
+	term1 = 3*density_l_e.*sum(bsxfun(@times,Galpha,(loads_element.*area_element)'),2)./earth_area; 
+	term2 = 3*density_o_e.*sum(bsxfun(@times,Galpha,(mean(rsl(index),2).*ocean_element.*area_element)'),2)./earth_area; 
+	vlm = term1+term2; 
+end
 
+