Index: /issm/trunk-jpl/src/m/psl/gmtOcean.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/gmtOcean.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/gmtOcean.m	(revision 19296)
@@ -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="./data_gmt/pixels/pix%dall.txt"',num); 
+text8 = sprintf('pixOce="./data_gmt/pixels/pix%doce.txt"',num); 
+text9 = sprintf('pixCon="./data_gmt/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('./data_gmt/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/psl/gmtResults.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/gmtResults.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/gmtResults.m	(revision 19296)
@@ -0,0 +1,126 @@
+function gmtResults(data_dir,seaLevel,ice_func,displace,farrell) 
+
+%---------------------------------------------------------------------
+% gmtResults :: a function to write a GMT code for plotting solutions 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+min_d = ceil(min(displace)); % for negative 
+max_d = floor(max(displace)); 
+ten_d = floor((max_d-min_d)/10);
+%
+min_s = ceil(min(seaLevel)); % for negative 
+max_s = floor(max(seaLevel)); 
+ten_s = floor((max_s-min_s)/10);
+%
+min_i = ceil(min(ice_func)); % for negative 
+max_i = floor(max(ice_func)); 
+ten_i = floor((max_i-min_i)/10);
+%
+ten_o = 1/10; 
+
+text0 = ''; 
+text00 = '   ;;'; 
+text01 = '   #'; 
+text1 = '#!/bin/sh'; 
+text2 = '# GMT code for plotting results (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('dG="./results_%s"',data_dir); 
+% 
+text20 = 'step=1'; 
+text21 = 'case $step in'; 
+text22 = '   1) # plot RSL for rigid Earth'; 
+text23 = sprintf('   gmt makecpt -Cno_green -T%d/%d/%d > $dG/pale_s.cpt',min_s,max_s,ten_s); 
+text24 = '   gmt psbasemap -Y19 -P -Bf180a90/f90a60WSEN -R0/360/-85/85 -JQ180/16 -K > $dG/seaLevel.ps'; 
+text25 = '   gmt pscontour -I -JQ -O -K -R0/360/-85/85 $dG/seaLevel.txt -C$dG/pale_s.cpt  >> $dG/seaLevel.ps'; 
+text25a = '   gmt pscontour -W0.5p -JQ -O -K -R0/360/-85/85 $dG/seaLevel.txt -C$dG/pale_s.cpt  >> $dG/seaLevel.ps'; 
+text26 = '   gmt pscoast -B -R -O -K -W0.1p -JQ -Df -A1000 >> $dG/seaLevel.ps'; 
+text26b = '   gmt pscoast -B -R -O -K -W0.1p -JQ -Df -Ggray -A1000 >> $dG/seaLevel.ps'; 
+text27 = '   gmt psscale -U/.2/.2/"ISSM-PSL" -E -K -C$dG/pale_s.cpt -Ba:relative_sea_level:/:m: -D8/-1/10/0.75h -O >> $dG/seaLevel.ps'; 
+% 
+text280 = sprintf('   gmt makecpt -Cno_green -T%d/%d/%d > $dG/pale_d.cpt',min_d,max_d,ten_d); 
+text281 = '   gmt psbasemap -Y-13 -P -Bf180a90/f90a60WSEN -R0/360/-85/85 -JQ180/16 -K -O >> $dG/seaLevel.ps'; 
+text282 = '   gmt pscontour -I -JQ -O -K -R0/360/-85/85 $dG/displace.txt -C$dG/pale_d.cpt  >> $dG/seaLevel.ps'; 
+text283 = '   gmt pscoast -B -R -O -K -W0.1p -JQ -Df -A1000 >> $dG/seaLevel.ps'; 
+text284 = '   gmt psscale -U/.2/.2/"ISSM-PSL" -E -K -C$dG/pale_d.cpt -Baf:vertical_displacement:/:m: -D8/-1/10/0.75h -O >> $dG/seaLevel.ps'; 
+text28 = '   open $dG/seaLevel.ps'; 
+% 
+text50 = '   2) # plot reconstructed ocean and ice load functions'; 
+text51 = sprintf('   gmt makecpt -Cno_green -T%d/%d/%d > $dG/pale_i.cpt',min_i,max_i,ten_i); 
+text52 = sprintf('   gmt makecpt -Cno_green -T0/1/%f > $dG/pale_o.cpt',ten_o); 
+text53 = '   gmt psbasemap -Y19 -P -Bf180a90/f90a60WSEN -R0/360/-85/85 -JQ180/16 -K > $dG/ioFunction.ps'; 
+text54 = '   gmt pscontour -I -JQ -O -K -R0/360/-85/85 $dG/ocean.txt -C$dG/pale_o.cpt  >> $dG/ioFunction.ps'; 
+text55 = '   gmt pscoast -B -R -O -K -W0.1p -JQ -Df -A1000 >> $dG/ioFunction.ps'; 
+text56 = '   gmt psscale -U/.2/.2/"ISSM-PSL" -E -K -C$dG/pale_o.cpt -Ba:Reconstructed_ocean_function: -D8/-1/10/0.75h -O >> $dG/ioFunction.ps'; 
+% 
+text57 = '   gmt psbasemap -Y-13 -P -Bf180a90/f90a60WSEN -R0/360/-85/85 -JQ180/16 -K -O >> $dG/ioFunction.ps'; 
+text58 = '   gmt pscontour -I -JQ -O -K -R0/360/-85/85 $dG/iceLoad.txt -C$dG/pale_i.cpt  >> $dG/ioFunction.ps'; 
+text59 = '   gmt pscoast -B -R -O -K -W0.1p -JQ -Df -A1000 >> $dG/ioFunction.ps'; 
+text60 = '   gmt psscale -U/.2/.2/"ISSM-PSL" -E -K -C$dG/pale_i.cpt -Baf:Reconstructed_ice_load:/:m: -D8/-1/10/0.75h -O >> $dG/ioFunction.ps'; 
+text61 = '   open $dG/ioFunction.ps'; 
+text99 = 'esac'; 
+% 
+fid = fopen(sprintf('./data_gmt/plotResults.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',text0); 
+fprintf(fid,'%s\n',text20); 
+fprintf(fid,'%s\n',text0); 
+fprintf(fid,'%s\n',text21); 
+fprintf(fid,'%s\n',text0); 
+fprintf(fid,'%s\n',text22); 
+fprintf(fid,'%s\n',text23); 
+fprintf(fid,'%s\n',text24); 
+if (farrell==1)  % Farrell and Clark (1976) benchmark  
+   fprintf(fid,'%s\n',text25a); 
+   fprintf(fid,'%s\n',text26b); 
+else 
+   fprintf(fid,'%s\n',text25); 
+   fprintf(fid,'%s\n',text26); 
+   fprintf(fid,'%s\n',text27); 
+end
+% 
+fprintf(fid,'%s\n',text01); 
+fprintf(fid,'%s\n',text280); 
+fprintf(fid,'%s\n',text281); 
+fprintf(fid,'%s\n',text282); 
+fprintf(fid,'%s\n',text283); 
+fprintf(fid,'%s\n',text284); 
+fprintf(fid,'%s\n',text28); 
+fprintf(fid,'%s\n',text00); 
+%
+fprintf(fid,'%s\n',text0); 
+fprintf(fid,'%s\n',text50); 
+fprintf(fid,'%s\n',text51); 
+fprintf(fid,'%s\n',text52); 
+fprintf(fid,'%s\n',text53); 
+fprintf(fid,'%s\n',text54); 
+fprintf(fid,'%s\n',text55); 
+fprintf(fid,'%s\n',text56); 
+fprintf(fid,'%s\n',text01); 
+fprintf(fid,'%s\n',text57); 
+fprintf(fid,'%s\n',text58); 
+fprintf(fid,'%s\n',text59); 
+fprintf(fid,'%s\n',text60); 
+fprintf(fid,'%s\n',text61); 
+%
+fprintf(fid,'%s\n',text0); 
+fprintf(fid,'%s\n',text99); 
+% 
+fclose(fid); 
+
Index: /issm/trunk-jpl/src/m/psl/icelm.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/icelm.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/icelm.m	(revision 19296)
@@ -0,0 +1,24 @@
+function ice_lm = icelm(lMax,nPix,ice,sh) 
+
+%---------------------------------------------------------------------
+% icelm :: a function to compute SH coefficients of "ice load 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; 
+
+ice_lm = zeros(1,(lMax+1)^2); 
+
+for l=0:lMax
+   for m=-l:l
+      ice_lm(1+p1) = sum(sh(:,1+p1).*ice');
+      p1 = p1+1;
+   end
+end
+
+ice_lm = ice_lm*4*pi/nPix; 
+
Index: /issm/trunk-jpl/src/m/psl/intro.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/intro.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/intro.m	(revision 19296)
@@ -0,0 +1,24 @@
+function intro(lMax,nPix,zlr_pq,oce_pq) 
+
+%---------------------------------------------------------------------
+% intro :: a function to instruct what software should be installed prior to using this code. 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+disp([' === A brief instruction before you can run this code =============================== ']); 
+disp(['     1. Install MEALPix http://sourceforge.net/projects/mealpix/ ']); 
+disp(['        MEALPix is a Matlab implementation of HEALPix http://healpix.jpl.nasa.gov/ ']); 
+disp(['        We rely on MEALPix mainly for pixelization & spherical harmonics computation. ']); 
+disp(['        I found a bug in MEALPix. E-mail for help: surendra.adhikari@jpl.nasa.gov ']); 
+disp(['     2. Install (optional) GMT http://gmt.soest.hawaii.edu/ ']); 
+disp(['        GMT-5.1.1 is recommended for direct use of our gmt code. ']); 
+disp(['     3. If MEALPix and GMT are both installed, ']); 
+disp(['        turn both "pix_switch" and "gmt_switch" on in ./data_mat/para.m and rerun. ']); 
+disp(['     4. If MEALPix is installed and there is no intention of using GMT, ']); 
+disp(['        turn only "pix_switch" on in ./data_mat/para.m and rerun. ']); 
+disp([' === E-mail to surendra.adhikari@jpl.nasa.gov for any question ====================== ']); 
+
Index: /issm/trunk-jpl/src/m/psl/lonLat.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/lonLat.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/lonLat.m	(revision 19296)
@@ -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/psl/oceFun.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/oceFun.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/oceFun.m	(revision 19296)
@@ -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/psl/ocelm.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/ocelm.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/ocelm.m	(revision 19296)
@@ -0,0 +1,23 @@
+function oce_lm = ocelm(lMax,nPix,ocean,sh) 
+
+%---------------------------------------------------------------------
+% ocelm :: a function to compute SH coefficients of "ocean 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; 
+
+oce_lm = zeros(1,(lMax+1)^2); 
+
+for l=0:lMax
+   for m=-l:l
+      oce_lm(1+p1) = sum(sh(:,1+p1).*ocean');
+      p1 = p1+1;
+   end
+end
+
+oce_lm = oce_lm*4*pi/nPix;
+
Index: /issm/trunk-jpl/src/m/psl/results.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/results.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/results.m	(revision 19296)
@@ -0,0 +1,38 @@
+function results(nSide,nPix,data_dir,ice_func,oce_func,eustatic,seaLevel,displace) 
+
+%---------------------------------------------------------------------
+% results :: save results that could be plotted using GMT code (cf. gmtResults.m) 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+lon_lat = lonLat(nSide,nPix);  % (lon, lat) 
+lonLat_ice = zeros(3,nPix); 
+lonLat_ice(1,:) = lon_lat(1,:); 
+lonLat_ice(2,:) = lon_lat(2,:); 
+lonLat_ice(3,:) = ice_func; 
+dlmwrite(sprintf('./data_gmt/results_%s/iceLoad.txt',data_dir),lonLat_ice','delimiter','\t'); 
+lonLat_oce   = zeros(3,nPix); 
+lonLat_oce(1,:) = lon_lat(1,:); 
+lonLat_oce(2,:) = lon_lat(2,:); 
+lonLat_oce(3,:) = oce_func; 
+dlmwrite(sprintf('./data_gmt/results_%s/ocean.txt',data_dir),lonLat_oce','delimiter','\t'); 
+lonLat_eus = zeros(3,nPix); 
+lonLat_eus(1,:) = lon_lat(1,:); 
+lonLat_eus(2,:) = lon_lat(2,:); 
+lonLat_eus(3,:) = eustatic; 
+dlmwrite(sprintf('./data_gmt/results_%s/eustatic.txt',data_dir),lonLat_eus','delimiter','\t'); 
+lonLat_slr = zeros(3,nPix); 
+lonLat_slr(1,:) = lon_lat(1,:); 
+lonLat_slr(2,:) = lon_lat(2,:);  
+lonLat_slr(3,:) = seaLevel ; 
+dlmwrite(sprintf('./data_gmt/results_%s/seaLevel.txt',data_dir),lonLat_slr','delimiter','\t'); 
+lonLat_dis = zeros(3,nPix); 
+lonLat_dis(1,:) = lon_lat(1,:); 
+lonLat_dis(2,:) = lon_lat(2,:);  
+lonLat_dis(3,:) = displace ; 
+dlmwrite(sprintf('./data_gmt/results_%s/displace.txt',data_dir),lonLat_dis','delimiter','\t'); 
+
Index: /issm/trunk-jpl/src/m/psl/selm.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/selm.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/selm.m	(revision 19296)
@@ -0,0 +1,16 @@
+function se_lm = selm(rho_i,rho_o,ice_lm,oce_lm) 
+
+%---------------------------------------------------------------------
+% euslm :: a function to compute SH coefficients for eustatic terms S^E 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+se_lm = zeros(1,length(oce_lm)); 
+
+rho_r = rho_i/rho_o; 
+se_lm(1) = -sqrt(4*pi)*rho_r*ice_lm(1)/oce_lm(1);
+
Index: /issm/trunk-jpl/src/m/psl/xplm.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/xplm.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/xplm.m	(revision 19296)
@@ -0,0 +1,26 @@
+function xp_pq = xplm(lMax,nPix,rho_i,rho_e,ice_pq,oce_pq,loveH,loveK) 
+
+%---------------------------------------------------------------------
+% xplm :: a function to compute SH coefficient for X'  
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+q1 = 0; 
+
+for p=0:lMax
+   for q=-p:p 
+      if (p<1) 
+         xp_pq(1+q1) = 3*(rho_i/rho_e)*ice_pq(1+q1)/(2*p+1);
+      else 
+         xp_pq(1+q1) = 3*(rho_i/rho_e)*ice_pq(1+q1)*(1+loveK(p)-loveH(p))/(2*p+1);
+      end
+      q1 = q1+1; 
+   end 
+end 
+
+xp_pq(1) = xp_pq(1) - sum(xp_pq.*oce_pq)/oce_pq(1);
+
Index: /issm/trunk-jpl/src/m/psl/xxlm.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/xxlm.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/xxlm.m	(revision 19296)
@@ -0,0 +1,36 @@
+function x_lm = xxlm(lMax,nPix,rho_i,rho_e,ice_pq,oce_pq,ocean,sh,loveH,loveK) 
+
+%---------------------------------------------------------------------
+% xxlm :: a function to compute SH coefficient for X 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+q1 = 0; 
+q2 = 0; 
+
+% obtain xp_pq 
+xp_pq = xplm(lMax,nPix,rho_i,rho_e,ice_pq,oce_pq,loveH,loveK); 
+
+% compute xp, i.e. x' 
+xp = zeros(1,nPix);
+for p=0:lMax
+   for q = -p:p
+      xp = xp + xp_pq(1+q1).*sh(:,1+q1)';
+      q1 = q1 + 1;
+   end
+end
+
+x_lm = zeros(1,(lMax+1)^2);
+for l=0:lMax
+   for m=-l:l
+      x_lm(1+q2) = sum((xp'.*sh(:,1+q2)).*ocean');
+      q2 = q2+1;
+   end
+end
+% 
+x_lm = x_lm*4*pi/nPix;
+
Index: /issm/trunk-jpl/src/m/psl/yplm.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/yplm.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/yplm.m	(revision 19296)
@@ -0,0 +1,26 @@
+function yp_pq = yplm(lMax,nPix,rho_o,rho_e,zlr_pq,oce_pq,loveH,loveK) 
+
+%---------------------------------------------------------------------
+% yplm :: a function to compute SH coefficients for Y' 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+q1 = 0; 
+
+for p=0:lMax
+   for q=-p:p 
+      if (p<1) 
+         yp_pq(1+q1) = 3*(rho_o/rho_e)*zlr_pq(1+q1)/(2*p+1);
+      else 
+         yp_pq(1+q1) = 3*(rho_o/rho_e)*zlr_pq(1+q1)*(1+loveK(p)-loveH(p))/(2*p+1);
+      end
+      q1 = q1+1; 
+   end 
+end 
+
+yp_pq(1) = yp_pq(1) - sum(yp_pq.*oce_pq)/oce_pq(1);
+
Index: /issm/trunk-jpl/src/m/psl/yylm.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/yylm.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/yylm.m	(revision 19296)
@@ -0,0 +1,37 @@
+function y_lm = yylm(lMax,nPix,rho_o,rho_e,zlr_pq,oce_pq,ocean,sh,loveH,loveK) 
+
+%---------------------------------------------------------------------
+% yylm :: a function to compute SH coefficient for Y  
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+q1 = 0; 
+q2 = 0; 
+
+% obtain yp_pq
+yp_pq = yplm(lMax,nPix,rho_o,rho_e,zlr_pq,oce_pq,loveH,loveK); 
+
+% compute yp, i.e. y'
+yp = zeros(1,nPix);
+for p=0:lMax
+   for q = -p:p
+      yp = yp + yp_pq(1+q1).*sh(:,1+q1)';
+      q1 = q1 + 1;
+   end
+end
+
+y_lm = zeros(1,(lMax+1)^2);
+% 
+for l=0:lMax
+   for m=-l:l
+      y_lm(1+q2) = sum((yp'.*sh(:,1+q2)).*ocean');
+      q2 = q2+1;
+   end
+end
+% 
+y_lm = y_lm*4*pi/nPix;
+
Index: /issm/trunk-jpl/src/m/psl/zelm.m
===================================================================
--- /issm/trunk-jpl/src/m/psl/zelm.m	(revision 19296)
+++ /issm/trunk-jpl/src/m/psl/zelm.m	(revision 19296)
@@ -0,0 +1,14 @@
+function ze_lm = zelm(rho_i,rho_o,ice_lm,oce_lm) 
+
+%---------------------------------------------------------------------
+% zelm :: a function to compute SH coefficients for eustatic terms Z^E 
+%---------------------------------------------------------------------
+% This code is written as a part of the ISSM-PSL project
+% (c) S. Adhikari 
+%     Jet Propulsion Laboratory, Caltech 
+%     November 3, 2014
+%---------------------------------------------------------------------
+
+rho_r = rho_i/rho_o; 
+ze_lm = - (rho_r*ice_lm(1)/oce_lm(1))*oce_lm;
+