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 25417)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/SESAWslr.m	(revision 25417)
@@ -0,0 +1,75 @@
+function [eust,rsl] = SESAWslr(index,greens,para) 
+%SESAWslr :: computes GRD slr due to applied surface loads based on SESAW method. 
+% 
+%Reference: Adhikari et al., 2016, GMD: https://doi.org/10.5194/gmd-9-1087-2016 
+%
+%   Usage:
+%      [eus,rsl]=greensfunctions(md.mesh.elements,greens,para) 
+%      
+%      greens.Grigid = Gravitational Green's function for the rigid planet. 
+%      greens.Gelast = Gravitational Green's function for the elastic planet.
+%      greens.Uelast = Deformational (radial displacement) Green's function.
+% 
+%      para.ocean_element = ocean funnction mapped onto elemental centroid. 
+%      para.loads_element = land loads (ice or water) [m] 
+%      para.area_element = area of elements [m^2]  
+%      para.earth_density = averae density of the solid earth [kg/m^3] 
+%      para.ocean_density = ocean water density [kg/m^3] 
+%      para.loads_density = land loads (ice or freshwater) density [kg/m^3] 
+%      para.rel_tol = relative tolerance for iterative SLR solver 
+%      para.solidearth = rheological model for solid Earth: 'rigid' or 'elastic' 
+%      para.isgrd_R = Do we want to account for Rotational feedback? 0 or 1. 
+
+area_element = para.area_element; 
+ocean_element = para.ocean_element; 
+loads_element = para.loads_element; 
+rho_o = para.ocean_density; 
+rho_l = para.loads_density; % land loads, ice or water... 
+rho_e = para.earth_density; 
+rel_tol = para.rel_tol; 
+
+% total Green's function for elastic earth, i.e. (1+k_l-h_l)
+if strcmpi(para.solidearth,'rigid') 
+	Galpha = greens.Grigid; 
+elseif strcmpi(para.solidearth,'elastic') 
+	Galpha = greens.Grigid + greens.Gelast- greens.Uelast;		
+else
+	error(['Unknown solidearth model:' para.solidearth, '; Should be rigid or elastic']); 
+end
+
+density_o_e = rho_o/rho_e; 
+density_l_e = rho_l/rho_e; 
+density_l_o = rho_l/rho_o; 
+
+% eustatic term 
+eust=-density_l_o*sum(loads_element.*area_element)/sum(ocean_element.*area_element); 
+
+term1=3*density_l_e.*sum(bsxfun(@times,Galpha,(loads_element.*area_element)'),2)./sum(area_element); 
+func3=mean(term1(index),2).*ocean_element;
+term3=sum(func3.*area_element)./sum(ocean_element.*area_element); 
+
+% computation of sea level change 
+rsl=eust+term1-term3;
+
+norm_diff = 10;			p = 0; 
+
+while norm_diff > rel_tol
+	norm_old = sqrt(sum(rsl.^2)); 
+	% 
+	term2=3*density_o_e.*sum(bsxfun(@times,Galpha,(mean(rsl(index),2).*ocean_element.*area_element)'),2)./sum(area_element); 
+	func4=mean(term2(index),2).*ocean_element;
+	term4=sum(func4.*area_element)./sum(ocean_element.*area_element); 
+	% 
+	rsl=eust+term1+term2-term3-term4; 
+	norm_new  = sqrt(sum(rsl.^2));
+	norm_diff = abs(norm_new-norm_old)./norm_old;
+	p = p+1;
+	if norm_diff > rel_tol
+		disp(['     iteration # ', num2str(p), ' :: difference in norm = ', num2str(norm_diff)]);
+	else
+		disp(['     iteration # ', num2str(p), ' :: difference in norm = ', num2str(norm_diff)]);
+		disp(['     solution converged! ']);
+	end
+end 
+
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/greensfunctions.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/greensfunctions.m	(revision 25417)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/greensfunctions.m	(revision 25417)
@@ -0,0 +1,50 @@
+function [Grigid,Gelastic,Uelastic] = greensfunctions(md_mesh,love_numbers) 
+%greensfunctions :: computes Greens funtions for rigid and elastic Earth. 
+%
+%   Usage:
+%      [Grigid,Gelastic,Uelastic]=greensfunctions(md_mesh,love_numbers)
+
+nv=md_mesh.numberofvertices; 
+ne=md_mesh.numberofelements; 
+
+% compute lat,long at elemental centroids. 
+[late,longe] = latelonge(md_mesh); 
+
+% love numbers. 
+loveH = love_numbers.h;  % radial displacement (height) 
+loveK = love_numbers.k;  % gravitational potential (phi) 
+love_num=length(loveH)-1; % maximum Legendre degree 
+
+Grigid=zeros(nv,ne);			Gelastic=zeros(nv,ne);			Uelastic=zeros(nv,ne);  
+loveH_inf=loveH-ones(length(loveH),1).*loveH(love_num+1); 
+loveK_inf=loveK-ones(length(loveK),1).*loveK(love_num+1); 
+
+phi1=zeros(ne,1);        lambda1=zeros(ne,1); 
+phi2=late/180*pi;			 lambda2=longe/180*pi; 
+delPhi=zeros(nv,1);		 delLambda=zeros(nv,1); 
+% refine in the nearfield.
+theta_rad_nearfield = [0:0.001:0.01]; 
+theta_rad_farfield = [0.01:0.01:180]; 
+theta_rad = unique([theta_rad_nearfield theta_rad_farfield])*pi/180;  
+xx=cos(theta_rad); 
+legendreP=p_polynomial_value(length(xx),love_num,xx'); 
+elast_loveK=(loveK(love_num+1).*0.5./sin(0.5.*theta_rad))'...
+	+ sum(bsxfun(@times,legendreP,loveK_inf'),2); 
+elast_loveH=(loveH(love_num+1).*0.5./sin(0.5.*theta_rad))'...
+	+ sum(bsxfun(@times,legendreP,loveH_inf'),2); 
+
+qq=1; 
+for j=1:nv 
+	phi1(:,1)=md_mesh.lat(j)./180.*pi;	lambda1(:,1)=md_mesh.long(j)./180.*pi; % size #vertices 
+	delPhi=abs(phi2-phi1);					delLambda=abs(lambda2-lambda1);
+	alpha=2.*asin(sqrt(sin(delPhi./2).^2+cos(phi1).*cos(phi2).*sin(delLambda./2).^2)); 
+	Grigid(j,:)=0.5./sin(0.5.*alpha); % analytical soln 
+	Gelastic(j,:)=interp1(xx,elast_loveK,cos(alpha)); 
+	Uelastic(j,:)=interp1(xx,elast_loveH,cos(alpha)); 
+	if (j==500*qq)
+		display([num2str(j),' of ', num2str(md_mesh.numberofvertices),' vertices done!']);
+		qq=qq+1; 
+	end 
+end  
+
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/latelonge.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/latelonge.m	(revision 25417)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/latelonge.m	(revision 25417)
@@ -0,0 +1,49 @@
+function [late,longe] = latelonge(md_mesh) 
+%latelonge :: computes lat,long at elemental centroids on spherical surface
+%
+%   Usage:
+%      [late,longe]=latelonge(md.mesh)
+
+nv = md_mesh.numberofvertices; 
+if length(md_mesh.lat)~=nv | length(md_mesh.long)~=nv 
+	error('lat,long not defined properly at vertices.');
+end
+
+% lat -> [0,180]; long -> [0,360] to compute centroids 
+lat=90-md_mesh.lat;		lon=md_mesh.long; 
+lon(lon<0)=180+(180+lon(lon<0)); 
+
+ax_0=lat(md_mesh.elements(:,1)); ay_0=lon(md_mesh.elements(:,1)); 
+bx_0=lat(md_mesh.elements(:,2)); by_0=lon(md_mesh.elements(:,2)); 
+cx_0=lat(md_mesh.elements(:,3)); cy_0=lon(md_mesh.elements(:,3)); 
+% find whether long is 0 or 360! This is important to compute centroids as well as elemental area 
+for ii=1:md_mesh.numberofelements
+	if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180)
+		if ay_0(ii)==0
+			ay_0(ii)=360;
+		end 
+		if by_0(ii)==0
+			by_0(ii)=360; 
+		end 
+		if cy_0(ii)==0 
+			cy_0(ii)=360; 
+		end
+	end 
+end
+% correction at the north pole 
+ay_0(ax_0==0)=(by_0(ax_0==0)+cy_0(ax_0==0))./2; 
+by_0(bx_0==0)=(cy_0(bx_0==0)+ay_0(bx_0==0))./2; 
+cy_0(cx_0==0)=(ay_0(cx_0==0)+by_0(cx_0==0))./2; 
+% correction at the south pole 
+ay_0(ax_0==180)=(by_0(ax_0==180)+cy_0(ax_0==180))./2; 
+by_0(bx_0==180)=(cy_0(bx_0==180)+ay_0(bx_0==180))./2; 
+cy_0(cx_0==180)=(ay_0(cx_0==180)+by_0(cx_0==180))./2; 
+% 
+late=(ax_0+bx_0+cx_0)/3; 
+longe=(ay_0+by_0+cy_0)/3;
+
+% back to [-90 90] [-180 180] ranges. 
+late = 90-late; 
+longe(longe>180) = longe(longe>180)-360; 
+
+
Index: sm/trunk-jpl/src/m/contrib/adhikari/latelonge.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/latelonge.m	(revision 25416)
+++ 	(revision )
@@ -1,49 +1,0 @@
-function [late,longe] = latelonge(md_mesh) 
-%latelonge :: computes lat,long at elemental centroids on spherical surface
-%
-%   Usage:
-%      [late,longe]=latelonge(md.mesh)
-
-nv = md_mesh.numberofvertices; 
-if length(md_mesh.lat)~=nv | length(md_mesh.long)~=nv 
-	error('lat,long not defined properly at vertices.');
-end
-
-% lat -> [0,180]; long -> [0,360] to compute centroids 
-lat=90-md_mesh.lat;		lon=md_mesh.long; 
-lon(lon<0)=180+(180+lon(lon<0)); 
-
-ax_0=lat(md_mesh.elements(:,1)); ay_0=lon(md_mesh.elements(:,1)); 
-bx_0=lat(md_mesh.elements(:,2)); by_0=lon(md_mesh.elements(:,2)); 
-cx_0=lat(md_mesh.elements(:,3)); cy_0=lon(md_mesh.elements(:,3)); 
-% find whether long is 0 or 360! This is important to compute centroids as well as elemental area 
-for ii=1:md_mesh.numberofelements
-	if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180)
-		if ay_0(ii)==0
-			ay_0(ii)=360;
-		end 
-		if by_0(ii)==0
-			by_0(ii)=360; 
-		end 
-		if cy_0(ii)==0 
-			cy_0(ii)=360; 
-		end
-	end 
-end
-% correction at the north pole 
-ay_0(ax_0==0)=(by_0(ax_0==0)+cy_0(ax_0==0))./2; 
-by_0(bx_0==0)=(cy_0(bx_0==0)+ay_0(bx_0==0))./2; 
-cy_0(cx_0==0)=(ay_0(cx_0==0)+by_0(cx_0==0))./2; 
-% correction at the south pole 
-ay_0(ax_0==180)=(by_0(ax_0==180)+cy_0(ax_0==180))./2; 
-by_0(bx_0==180)=(cy_0(bx_0==180)+ay_0(bx_0==180))./2; 
-cy_0(cx_0==180)=(ay_0(cx_0==180)+by_0(cx_0==180))./2; 
-% 
-late=(ax_0+bx_0+cx_0)/3; 
-longe=(ay_0+by_0+cy_0)/3;
-
-% bak to [-90 90] [-180 180] ranges. 
-late = 90-late; 
-longe(longe>180) = longe(longe>180)-360; 
-
-