Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/SHslr.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/SHslr.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/SHslr.m	(revision 25369)
@@ -0,0 +1,103 @@
+function [barystatic, sealevel, geoid, bed] = SHslr(sh,para,force,ocean,rot_flag); 
+
+% slrlm :: a function to compute SH coefficients of slr fields. 
+% reference. Adhikari et al., 2019, ESSD, https://doi.org/10.5194/essd-11-629-2019 
+	
+	% retrieve SH parameters 
+	[nPix num_lm] = size(sh); 
+	lMax = sqrt(num_lm)-1; 
+
+	% SH coefficients of forcing function and ocean function. 
+	force_lm = funlm(lMax,nPix,force,sh);  
+	oce_lm = funlm(lMax,nPix,ocean,sh); 
+
+	% allocate. {{{
+	%sealevel = zeros(num_time,nPix);		% relative sea-level
+	%geoid = zeros(num_time,nPix);		% geoid height  
+	%bed = zeros(num_time,nPix);		% bed height 
+	%sealevel_lm = zeros(num_time,num_lm);	% relative sea-level SH coefficients 
+	%geoid_lm = zeros(num_time,num_lm);	% geoid height SH coefficients 
+	%bed_lm = zeros(num_time,num_lm);		% bed height SH coefficients 
+	% }}} 
+
+	% Compute SH coefficients for shat  {{{ 
+	ehat_lm = ehatlm(force_lm,oce_lm);  % starting solution for shat_lm (eq. B17 of Adhikari&Ivins, 2018, ESSD) 
+	shat_lm = ehat_lm;  % initialize shat_lm 
+	norm_diff = 1.0;
+	p = 0;
+	while norm_diff > para.rel_tol  
+		
+		load_lm = force_lm + shat_lm; 
+
+		norm_old = sqrt(sum(shat_lm.^2)); 
+		xhat_lm = xhatlm(load_lm,ocean,lMax,nPix,sh,para); 
+		phat_lm = phatlm(load_lm,ocean,lMax,nPix,sh,para);  
+		chat_lm = chatlm(force_lm,shat_lm,oce_lm,lMax,para,rot_flag); 
+		
+		if (rot_flag==0)
+			shat_lm = xhat_lm + phat_lm + chat_lm; 
+		elseif (rot_flag==1) 
+			yhat_lm = yhatlm(load_lm,ocean,lMax,nPix,sh,para); 
+			qhat_lm = qhatlm(load_lm,ocean,lMax,nPix,sh,para);  
+			shat_lm = xhat_lm + phat_lm + chat_lm + yhat_lm + qhat_lm; 
+		else 
+			error('rotational flag should be set either to 0 or 1'); 
+		end
+		
+		norm_new  = sqrt(sum(shat_lm.^2));
+		norm_diff = abs(norm_new-norm_old)/abs(norm_old); % relative change in solutions...  
+		p = p+1; 
+		if norm_diff > para.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 
+	%}}}
+
+	% Compute SH coefficients for slr {{{ 
+	x_lm = xpq(load_lm,lMax,para);  
+	p_lm = ppq(load_lm,lMax,para); 
+	c_lm = cpq(force_lm,shat_lm,oce_lm,lMax,para,rot_flag); 
+
+	if (rot_flag==0); 
+		slr_lm = x_lm + p_lm + c_lm; 
+		geoid_lm = x_lm;	% \Phi /g [m]  
+		bed_lm = -(p_lm);	% U [m] 
+	elseif (rot_flag==1) 
+		y_lm = ypq(load_lm,lMax,para); 
+		q_lm = qpq(load_lm,lMax,para); 
+		slr_lm = x_lm + p_lm + c_lm + y_lm + q_lm; 
+		geoid_lm = x_lm+y_lm;	% \Phi /g [m]  
+		bed_lm = -(p_lm+q_lm);	% U [m] 
+	else 
+		error('rotational flag should be set either to 0 or 1'); 
+	end
+	%}}}
+
+	% Compute final solutions {{{ 
+	p = 0;
+	sealevel = zeros(1,nPix);  % self-gravitating (relative) sea level 
+	geoid = zeros(1,nPix);  % self-gravitating geoid height 
+	bed = zeros(1,nPix);  % self-gravitating bed height 
+	for l=0:lMax 
+		for m = -l:l
+			sealevel = sealevel + slr_lm(1+p).*sh(:,1+p)'; 
+			geoid = geoid + geoid_lm(1+p).*sh(:,1+p)'; 
+			bed = bed + bed_lm(1+p).*sh(:,1+p)'; 
+			p = p+1;
+		end
+	end 
+	%}}}
+
+	% Test whether mass is conserved (to ensure the solution is trustworthy!) {{{ 
+	% ocean-average sea level and eustatic value must be equal! 
+	avg_eustatic = -force_lm(1)/oce_lm(1); % see equation B16 
+	barystatic = sum(sealevel.*ocean)/nnz(ocean); % we can do this because of equal-area pixelization 
+	sol_diff = abs(abs(barystatic) - abs(avg_eustatic))/abs(barystatic)*100; % per cent 
+	if (sol_diff > 1) % per cent. 
+		error('Ocean-average (eustatic - sealevel) is NOT negligible. Returning.'); 
+	end
+	%}}} 
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/chatlm.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/chatlm.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/chatlm.m	(revision 25369)
@@ -0,0 +1,16 @@
+function chat_lm = chatlm(force_pq,shat_pq,oce_pq,lMax,para,rotation) 
+
+%---------------------------------------------------------------------
+% chatlm :: a function to compute SH coefficients for \hat{C} 
+% Adhikari & Ivins, ESSD, 2018: Equation B17 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+% obtain c_pq 
+c_pq = cpq(force_pq,shat_pq,oce_pq,lMax,para,rotation); 
+
+chat_lm = c_pq.*oce_pq; 
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/cpq.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/cpq.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/cpq.m	(revision 25369)
@@ -0,0 +1,30 @@
+function c_pq = cpq(force_pq,shat_pq,oce_pq,lMax,para,rotation) 
+
+%---------------------------------------------------------------------
+% cpq :: a function to compute SH coefficient for C 
+% Adhikari & Ivins, ESSD, 2018: Equation B17 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+q1 = 0; 
+
+load_pq = force_pq+shat_pq; 
+x_pq = xpq(load_pq,lMax,para); 
+p_pq = ppq(load_pq,lMax,para); 
+y_pq = ypq(load_pq,lMax,para); 
+q_pq = qpq(load_pq,lMax,para);  
+
+if (rotation==0); 
+	c_pq_0 = x_pq+p_pq; 
+elseif (rotation==1); 
+	c_pq_0 = x_pq+p_pq + y_pq+q_pq; 
+else 
+	error('rotational flag should be set either to 0 or 1'); 
+end
+
+c_pq = zeros(1,(lMax+1)^2);
+c_pq(1) = -force_pq(1)/oce_pq(1) -sum(c_pq_0.*oce_pq)/oce_pq(1);  
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/ehatlm.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/ehatlm.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/ehatlm.m	(revision 25369)
@@ -0,0 +1,13 @@
+function ehat_lm = ehatlm(force_lm,oce_lm) 
+
+%---------------------------------------------------------------------
+% ehatlm :: a function to compute SH coefficients for eustatic terms \hat{E}
+% Adhikari & Ivins, ESSD, 2018: Equation B17 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+ehat_lm = -(force_lm(1)/oce_lm(1))*oce_lm;
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/funlm.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/funlm.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/funlm.m	(revision 25369)
@@ -0,0 +1,23 @@
+function fun_lm = funlm(lMax,nPix,fun,sh) 
+
+%---------------------------------------------------------------------
+% funlm :: a function to compute SH coefficients of "fun" function 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+p1 = 0; 
+
+fun_lm = zeros(1,(lMax+1)^2); 
+
+for l=0:lMax
+   for m=-l:l
+      fun_lm(1+p1) = sum(sh(:,1+p1).*fun');
+      p1 = p1+1;
+   end
+end
+
+fun_lm = fun_lm/nPix;
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/gmtOcean.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/gmtOcean.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/gmtOcean.m	(revision 25369)
@@ -0,0 +1,43 @@
+function gmtOce(num) 
+
+%---------------------------------------------------------------------
+% gmtOce :: a function to write a GMT code for selecting pixels (e.g., ocean function) 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+text0 = ''; 
+text1 = '#!/bin/sh'; 
+text2 = '# GMT code for selecting ocean and continental data (from Selen)'; 
+text3 = '# This code is written as a part of the ISSM-PSL project'; 
+text4 = '# (c) Surendra Adhikari'; 
+text5 = '#     Jet Propulsion Laboratory, Caltech'; 
+text6 = '#     November 3, 2014'; 
+text7 = sprintf('pixAll="./gmtdir/pixels/pix%dall.txt"',num); 
+text8 = sprintf('pixOce="./gmtdir/pixels/pix%doce.txt"',num); 
+text9 = sprintf('pixCon="./gmtdir/pixels/pix%dcon.txt"',num); 
+text10 = '# Select ocean (and continent) pixels'; 
+text11 = 'gmt gmtselect $pixAll -h0 -Df -R0/360/-90/90  -A0 -JQ180/200 -Nk/s/k/s/k > $pixOce'; 
+text12 = '#gmt gmtselect $pixAll -h0 -Df -R              -A0 -JQ        -Ns/k/s/k/s > $pixCon'; 
+% 
+fid = fopen(sprintf('./gmtdir/oceanFunc.sh'),'w');  % "./" = ISSM_PSL directory! 
+fprintf(fid,'%s\n',text1); 
+fprintf(fid,'%s\n',text0); 
+fprintf(fid,'%s\n',text2); 
+fprintf(fid,'%s\n',text3); 
+fprintf(fid,'%s\n',text4); 
+fprintf(fid,'%s\n',text5); 
+fprintf(fid,'%s\n',text6); 
+fprintf(fid,'%s\n',text0); 
+fprintf(fid,'%s\n',text7); 
+fprintf(fid,'%s\n',text8); 
+fprintf(fid,'%s\n',text9); 
+fprintf(fid,'%s\n',text0); 
+fprintf(fid,'%s\n',text10); 
+fprintf(fid,'%s\n',text11); 
+fprintf(fid,'%s\n',text12); 
+fclose(fid); 
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/lonLat.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/lonLat.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/lonLat.m	(revision 25369)
@@ -0,0 +1,20 @@
+function lon_lat = lonLat(nSide,nPix)
+
+%---------------------------------------------------------------------
+% lonLat :: a function to compute (long, lat) in degrees for GMT 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+thetaLambda = pix2ang(nSide); 
+data = cell2mat(thetaLambda); 
+data = data*180/pi;  % in degrees 
+data(1,:) = -data(1,:)+90;  % lat \in [-90, 90] in GMT 
+% 
+lon_lat = zeros(2,nPix); 
+lon_lat(1,:) = data(2,:); 
+lon_lat(2,:) = data(1,:); 
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/oceFun.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/oceFun.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/oceFun.m	(revision 25369)
@@ -0,0 +1,22 @@
+function ocean = oceFun(nSide,nPix,oce_pix) 
+
+%---------------------------------------------------------------------
+% oceFun :: a function to compute ocean function 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+oce_fun = oce_pix';
+oce_fun(2,:) = -(oce_fun(2,:)-90);  % lat in [0,180] from North pole 
+oce_fun = oce_fun*pi/180;
+oce_fun_latLon = zeros(size(oce_fun));  % write in lat,long now 
+oce_fun_latLon(1,:) = oce_fun(2,:);
+oce_fun_latLon(2,:) = oce_fun(1,:);
+ocean = zeros(1,nPix);  % initialize the ocean function 
+oce_array = arrayfun(@(x,y)([x;y]),oce_fun_latLon(1,:),oce_fun_latLon(2,:),'UniformOutput',false); 
+ocean_res = ang2pix(nSide,oce_array); 
+ocean(ocean_res) = 1.0; 
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/phatlm.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/phatlm.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/phatlm.m	(revision 25369)
@@ -0,0 +1,36 @@
+function phat_lm = phatlm(load_pq,ocean,lMax,nPix,sh,para)
+
+%---------------------------------------------------------------------
+% phatlm :: a function to compute SH coefficients for \hat{P} 
+% Adhikari & Ivins, ESSD, 2018: Equation B17 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+q1 = 0; 
+q2 = 0; 
+
+% obtain p_pq 
+p_pq = ppq(load_pq,lMax,para); 
+
+% compute xp, i.e. x' 
+xp = zeros(1,nPix);
+for p=0:lMax
+   for q = -p:p
+      xp = xp + p_pq(1+q1).*sh(:,1+q1)';
+      q1 = q1 + 1;
+   end
+end
+
+phat_lm = zeros(1,(lMax+1)^2);
+for l=0:lMax
+   for m=-l:l
+      phat_lm(1+q2) = sum((xp'.*sh(:,1+q2)).*ocean');
+      q2 = q2+1;
+   end
+end
+% 
+phat_lm = phat_lm/nPix;
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/ppq.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/ppq.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/ppq.m	(revision 25369)
@@ -0,0 +1,24 @@
+function p_pq = ppq(load_pq,lMax,para) 
+
+%---------------------------------------------------------------------
+% ppq :: a function to compute SH coefficient for P
+% Adhikari & Ivins, ESSD, 2018: Equation B17 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+q1 = 0; 
+
+for p=0:lMax
+   for q=-p:p 
+      if (p<1) 
+         p_pq(1+q1) = 0; 
+      else 
+         p_pq(1+q1) = -3*(para.rho_ocean/para.rho_earth)*load_pq(1+q1)*para.loveH(p+1)/(2*p+1);
+      end
+      q1 = q1+1; 
+   end 
+end 
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/qhatlm.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/qhatlm.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/qhatlm.m	(revision 25369)
@@ -0,0 +1,36 @@
+function qhat_lm = qhatlm(load_pq,ocean,lMax,nPix,sh,para) 
+
+%---------------------------------------------------------------------
+% qhatlm :: a function to compute SH coefficients for \hat{Q} 
+% Adhikari & Ivins, ESSD, 2018: Equation B17 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+q1 = 0; 
+q2 = 0; 
+
+% obtain q_pq 
+q_pq = qpq(load_pq,lMax,para); 
+
+% compute xp, i.e. x' 
+xp = zeros(1,nPix);
+for p=0:lMax
+   for q = -p:p
+      xp = xp + q_pq(1+q1).*sh(:,1+q1)';
+      q1 = q1 + 1;
+   end
+end
+
+qhat_lm = zeros(1,(lMax+1)^2);
+for l=0:lMax
+   for m=-l:l
+      qhat_lm(1+q2) = sum((xp'.*sh(:,1+q2)).*ocean');
+      q2 = q2+1;
+   end
+end
+% 
+qhat_lm = qhat_lm/nPix;
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/qpq.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/qpq.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/qpq.m	(revision 25369)
@@ -0,0 +1,56 @@
+function q_pq = qpq(load_pq,lMax,para) 
+
+%---------------------------------------------------------------------
+% qpq :: a function to compute SH coefficient for Q 
+% Adhikari & Ivins, ESSD, 2018: Equation B21 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+a = para.earth_radius; 
+rho_o = para.rho_ocean;
+Omega = para.Omega; 
+loveK = para.loveK;
+A = para.A;
+C = para.C;
+g = para.g; 
+h2 = para.h2;
+k2 = para.k2;
+ks = para.ks; 
+
+% Chandler wobble frequency 
+sigma_0 = (1-k2/ks)*Omega*(C-A)/A; 
+
+dI_13 = -(4*pi/sqrt(15)) *rho_o *a^4 *load_pq(8); 
+dI_23 = -(4*pi/sqrt(15)) *rho_o *a^4 *load_pq(6); 
+dI_33 = (8*pi/3) *rho_o *a^4 *(load_pq(1) - (1/sqrt(5))*load_pq(7)); 
+
+m1 = dI_13 * Omega*(1+loveK(3))/(A*sigma_0);  % loveK(3) is degree-2 load Love number 
+m2 = dI_23 * Omega*(1+loveK(3))/(A*sigma_0); 
+m3 = dI_23 * -(1+loveK(3))/C; 
+
+lambda_00 = (2/3) * a^2 * Omega^2 * m3;  
+lambda_20 = (-2/(3*sqrt(5))) *a^2 *Omega^2 *m3;  
+lambda_21p = (-1/sqrt(15)) *a^2 *Omega^2 *m1;  
+lambda_21m = (-1/sqrt(15)) *a^2 *Omega^2 *m2;  
+
+q_pq = zeros(1,(lMax+1)^2);
+
+q1 = 0; 
+for p=0:lMax
+   for q=-p:p 
+		if (p==2)
+			if (q==-1)
+				q_pq(1+q1) = -h2*lambda_21m/g; 
+			elseif (q==0)
+				q_pq(1+q1) = -h2*lambda_20/g; 
+			elseif (q==1)
+				q_pq(1+q1) = -h2*lambda_21p/g; 
+			end
+      end
+      q1 = q1+1; 
+   end 
+end 
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/xhatlm.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/xhatlm.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/xhatlm.m	(revision 25369)
@@ -0,0 +1,37 @@
+function xhat_lm = xhatlm(load_pq,ocean,lMax,nPix,sh,para)
+
+%---------------------------------------------------------------------
+% xhatlm :: a function to compute SH coefficients for \hat{X} 
+% Adhikari & Ivins, ESSD, 2018: Equation B17 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+q1 = 0; 
+q2 = 0; 
+
+% obtain x_pq 
+x_pq = xpq(load_pq,lMax,para); 
+
+% compute xp, i.e. x' 
+xp = zeros(1,nPix);
+for p=0:lMax
+   for q = -p:p
+      xp = xp + x_pq(1+q1).*sh(:,1+q1)';
+      q1 = q1 + 1;
+   end
+end
+
+xhat_lm = zeros(1,(lMax+1)^2);
+for l=0:lMax
+   for m=-l:l
+      xhat_lm(1+q2) = sum((xp'.*sh(:,1+q2)).*ocean');
+      q2 = q2+1;
+   end
+end
+% 
+xhat_lm = xhat_lm/nPix;
+
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/xpq.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/xpq.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/xpq.m	(revision 25369)
@@ -0,0 +1,26 @@
+function x_pq = xpq(load_pq,lMax,para) 
+
+%---------------------------------------------------------------------
+% xpq :: a function to compute SH coefficient for X 
+% Adhikari & Ivins, ESSD, 2018: Equation B17 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+q1 = 0;
+rho_o2e = para.rho_ocean/para.rho_earth; 
+loveK = para.loveK; 
+
+for p=0:lMax
+   for q=-p:p 
+      if (p<1) 
+         x_pq(1+q1) = 3*rho_o2e*load_pq(1+q1)/(2*p+1);
+      else 
+         x_pq(1+q1) = 3*rho_o2e*load_pq(1+q1)*(1+loveK(p+1))/(2*p+1);
+      end
+      q1 = q1+1; 
+   end 
+end 
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/yhatlm.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/yhatlm.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/yhatlm.m	(revision 25369)
@@ -0,0 +1,36 @@
+function yhat_lm = yhatlm(load_pq,ocean,lMax,nPix,sh,para) 
+
+%---------------------------------------------------------------------
+% yhatlm :: a function to compute SH coefficients for \hat{Y} 
+% Adhikari & Ivins, ESSD, 2018: Equation B17 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+q1 = 0; 
+q2 = 0; 
+
+% obtain y_pq 
+y_pq = ypq(load_pq,lMax,para); 
+
+% compute xp, i.e. x' 
+xp = zeros(1,nPix);
+for p=0:lMax
+   for q = -p:p
+      xp = xp + y_pq(1+q1).*sh(:,1+q1)';
+      q1 = q1 + 1;
+   end
+end
+
+yhat_lm = zeros(1,(lMax+1)^2);
+for l=0:lMax
+   for m=-l:l
+      yhat_lm(1+q2) = sum((xp'.*sh(:,1+q2)).*ocean');
+      q2 = q2+1;
+   end
+end
+% 
+yhat_lm = yhat_lm/nPix;
+
Index: /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/ypq.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/ypq.m	(revision 25369)
+++ /issm/trunk-jpl/src/m/contrib/adhikari/sphericalharmonic_slr_solver/ypq.m	(revision 25369)
@@ -0,0 +1,57 @@
+function y_pq = ypq(load_pq,lMax,para)
+
+%---------------------------------------------------------------------
+% ypq :: a function to compute SH coefficient for Y
+% Adhikari & Ivins, ESSD, 2018: Equation B21 
+%---------------------------------------------------------------------
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     October 29, 2018
+%---------------------------------------------------------------------
+
+a = para.earth_radius; 
+rho_o = para.rho_ocean;
+Omega = para.Omega; 
+loveK = para.loveK;
+A = para.A;
+C = para.C;
+g = para.g; 
+k2 = para.k2;
+ks = para.ks; 
+
+% Chandler wobble frequency 
+sigma_0 = (1-k2/ks)*Omega*(C-A)/A; 
+
+dI_13 = -(4*pi/sqrt(15)) *rho_o *a^4 *load_pq(8); 
+dI_23 = -(4*pi/sqrt(15)) *rho_o *a^4 *load_pq(6); 
+dI_33 = (8*pi/3) *rho_o *a^4 *(load_pq(1) - (1/sqrt(5))*load_pq(7)); 
+
+m1 = dI_13 * Omega*(1+loveK(3))/(A*sigma_0);  % loveK(3) is degree-2 load Love number 
+m2 = dI_23 * Omega*(1+loveK(3))/(A*sigma_0); 
+m3 = dI_23 * -(1+loveK(3))/C; 
+
+lambda_00 = (2/3) * a^2 * Omega^2 * m3;  
+lambda_20 = (-2/(3*sqrt(5))) *a^2 *Omega^2 *m3;  
+lambda_21p = (-1/sqrt(15)) *a^2 *Omega^2 *m1;  
+lambda_21m = (-1/sqrt(15)) *a^2 *Omega^2 *m2;  
+
+y_pq = zeros(1,(lMax+1)^2);
+
+q1 = 0; 
+for p=0:lMax
+   for q=-p:p 
+      if (p==0) 
+         y_pq(1+q1) = lambda_00/g; 
+		elseif (p==2)
+			if (q==-1)
+				y_pq(1+q1) = (1+k2)*lambda_21m/g; 
+			elseif (q==0)
+				y_pq(1+q1) = (1+k2)*lambda_20/g; 
+			elseif (q==1)
+				y_pq(1+q1) = (1+k2)*lambda_21p/g; 
+			end
+      end
+      q1 = q1+1; 
+   end 
+end 
+