
name='Earth.PSL.';

org=organizer('steps',steps,...
'download',0,'repository','./Models','prefix',name,...
'upload_server','astrid.jpl.nasa.gov','upload_login','larour','upload_path','/u/astrid-r1b/ModelData/Models');

clustername=oshostname();
yts=365.25*24*3600;

printresults=0;
loadresults=1;

% Data paths %{{{
addpath ./Data ./Par ./Functions/ ./Results 
para; % load parameters 
%}}}
% Cluster parameters{{{
cluster_pfe=pfe;
cluster_pfe.login='elarour';
cluster_pfe.numnodes=2;
cluster_pfe.cpuspernode=12;
cluster_pfe.time=120*60;
cluster_pfe.interactive=2; %2; %6;
cluster_pfe.port=1025;
cluster_pfe.processor='ivy';
cluster_pfe.codepath='/nobackupnfs2/elarour/trunk-jpl-adolc/bin';
cluster_pfe.executionpath='/nobackupp8/elarour/execution/pleiades';

cluster_generic=generic('name',oshostname,'np',12,'verbose',0);
if strcmpi(clustername,'pfe'),
	cluster=cluster_pfe;

elseif strcmpi(clustername,oshostname()),
	cluster=cluster_generic;
else
	error('cluster not supported yet');
end% }}}

% setup of Earth (spherical mesh using gmsh) 
if perform(org,'Mesh') % {{{

	numrefine=2;
	resolution=500; %km inital resolution of the mesh prior to refinement.
	%numrefine=1;
	%resolution=200; %km inital resolution of the mesh prior to refinement.
	radius=6000; %km
	mindistance=50000; %m

	%mesh earth: 
	md=model; 
	md.mask=maskpsl(); % use maskpsl class (instead of mask) to store the ocean function as a ocean_levelset 
	md.mesh=gmshplanet('radius',radius,'resolution',resolution);
	md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 

	for i=1:numrefine,
	
		%figure out mask: 
		md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 

		%figure out distance to the coastline, in lat,long (not x,y,z): 
		% {{{
		distance=zeros(md.mesh.numberofvertices,1);

		pos=find(~md.mask.ocean_levelset);	coaste.lat=md.mesh.lat(pos);	coaste.long=md.mesh.long(pos);  
		pos=find(md.mask.ocean_levelset);	coasto.lat=md.mesh.lat(pos);	coasto.long=md.mesh.long(pos);  

		for j=1:md.mesh.numberofvertices,
			%figure out nearest coastline (using the great circle distance)
			phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi; 
			if md.mask.ocean_levelset(j),
				phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi; 
				deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
				d=radius*1000*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2));
			else
				phi2=coasto.lat/180*pi; lambda2=coasto.long/180*pi; 
				deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
				d=radius*1000*2*asin(sqrt(sin(deltaphi/2).^2+cos(phi1).*cos(phi2).*sin(deltalambda/2).^2));
			end
			distance(j)=min(d);
		end
		pos=find(distance<mindistance); distance(pos)=mindistance;
		%}}}

		%use distance to the coastline to refine mesh: 
		md.mesh=gmshplanet('radius',6000,'resolution',resolution,'refine',md.mesh,'refinemetric',distance/2);
	end

	%figure out mask: 
	md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 
	%savemodel(org,md);
end
% }}}
% Gravitationally consistent sea-level and geodetic signatures on a rotating, elastic earth 
% Here we compute solutions for ANTARCTIC ICE SHEET, whose mass balance is infered from GRACE 
if perform(org,'PSL_rotating_AIS') %  {{{
	
	%load the mesh and define a few things {{{ 
	if (exist('md','var')~=1 || exist('love','var')~=1)
		md=loadmodel(org,'Mesh'); 
		love = dlmread('love_numbers_10k.txt'); 
		loveH = love(:,2);  % radial displacement (height) 
		loveK = love(:,3);  % gravitational potential (phi) 
	end

	ne=md.mesh.numberofelements;	% number of elements 
	nv=md.mesh.numberofvertices;	% number of vertices 
	rho_i=md.materials.rho_ice;				rho_o=md.materials.rho_water; 
	ice_rat=md.materials.rho_ice/rho_e;		oce_rat=md.materials.rho_water/rho_e; 
	index=md.mesh.elements; 
	sigma_r=omega*(moi_C-moi_A)/moi_A;				% sigma r  
	sigma_s=sigma_r*(1-tideLoveK/fluidLove);		% sigma r star 
	%}}} 

	%compute centroids using (lat,long) data {{{ 
	
	% 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:ne
		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; 
	% 
	lat_element=(ax_0+bx_0+cx_0)/3; 
	lon_element=(ay_0+by_0+cy_0)/3;

	% }}} 
	
	%compute area of elements {{{ 
	if (exist('areas','var')~=1)
		areas=GetAreas3DTria(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z); 
	end 
	% }}} 

	if (exist('oce','var')~=1)
		disp([' === Computing ocean function ======================================== ']);
	% define ocean function {{{ 
	
	% to find a ocean function, [lat,lon] \in [-90:90,-180,180]; 
	lat_element_0 = 90-lat_element;		lon_element_0 = lon_element;
	lon_element_0(lon_element>180) = (lon_element(lon_element>180)-180) - 180;

	oce=gmtmask(lat_element_0,lon_element_0); 

	% }}} 
		disp(['     done! ']);	
	end 

	if (exist('grace_monthly','var')~=1 || exist('ice_load','var')~=1) 
	% load the GRACE monthly (ice) LOAD (in meters) {{{  
	%load the grace monthly data 
	time=dlmread('grace_months.txt')./12;		% years 
	grace_monthly=dlmread('grace_monthly_load.txt'); 
	% fill out the blanks {{{ 
	
	lat_grace=90-grace_monthly(:,1); 
	lon_grace=grace_monthly(:,2); 
	num_org=length(lon_grace); 

	qq=1;			mm=1; 
	for jj=2:num_org-1
		if (lat_grace(jj)~=lat_grace(jj+1))
			lat_new(qq)=lat_grace(jj); 
			lon_new(qq)=lon_grace(jj)+(lon_grace(jj)-lon_grace(jj-1)); 
			load_new(qq,1:141)=grace_monthly(mm,3:143); 
			lat_new(qq+1)=lat_grace(jj); 
			lon_new(qq+1)=lon_grace(mm)-(lon_grace(jj)-lon_grace(jj-1)); 
			load_new(qq+1,1:141)=grace_monthly(jj,3:143); 
			qq=qq+2; 
			mm=jj+1; % to find out the value for monthly data 
		end
	end
	
	num_add=length(lat_new); 
	num_plus=num_org+num_add;

	lat_grace_plus=zeros(num_plus,1); 
	lon_grace_plus=zeros(num_plus,1); 
	load_grace_plus=zeros(num_plus,141); 
	
	lat_grace_plus(1:num_org)=lat_grace;
	lat_grace_plus(1+num_org:num_plus)=lat_new; 
	lon_grace_plus(1:num_org)=lon_grace;
	lon_grace_plus(1+num_org:num_plus)=lon_new; 
	load_grace_plus(1:num_org,:)=grace_monthly(:,3:143);
	load_grace_plus(1+num_org:num_plus,:)=load_new(:,:);
	% }}}
	
	ice_load=zeros(ne,141);
	for jj=1:141
		ice_load(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(:,jj),lat_element,lon_element);
	end 

	% antarctic mask 
	rad=25*pi/180; 
	phi0=-90*pi/180;		lambda0=0*pi/180; 
	mask_ant=zeros(ne,1); for j=1:ne 
		phi=lat_element_0(j)/180*pi;	lambda=lon_element_0(j)/180*pi; 
		delPhi=abs(phi-phi0);			delLambda=abs(lambda-lambda0);
		dist=2*asin(sqrt(sin(delPhi/2).^2+cos(phi).*cos(phi0).*sin(delLambda/2).^2));
		if (dist <= rad) 
			mask_ant(j)= (oce(j)~=1);
		end 
	end
	
	% ice load antarctica 
	ice_ant=bsxfun(@times,ice_load,mask_ant);
	ice_ant(isnan(ice_ant))=0; 
	
	% monthly ice volume/mass 
	%Gt_ice = sum(bsxfun(@times,ice_ant,areas),1).*1e-9;	% Gt of ice 
	%test=fit(time',Gt_ice','poly1'); 
	ice_ant=ice_ant.*90/127.3;		% Gt of ice (simple rescaling)!!! 
	% ice thickness distribution antarctica wrt April, 2002 
	ice_wrt0=bsxfun(@minus,ice_ant,ice_ant(:,1))*rho_o/rho_i; 
	% }}} 
	end

	if (exist('G_rigid','var')~=1 || exist('G_elast','var')~=1 || exist('U_elast','var')~=1)
		disp([' === Computing Green functions ======================================= ']);
	% green's function computation {{{ 
	love_num=10000; 
	love_K=loveK(1:love_num+1);	love_H=loveH(1:love_num+1); 
	G_rigid=zeros(nv,ne);			G_elast=zeros(nv,ne);			U_elast=zeros(nv,ne);  
	love_H_inf=love_H-ones(length(love_H),1).*loveH(love_num+1); 
	love_K_inf=love_K-ones(length(love_K),1).*loveK(love_num+1); 

	phi1=zeros(ne,1);					lambda1=zeros(ne,1); 
	phi2=lat_element_0./180.*pi;	lambda2=lon_element_0./180.*pi; 
	delPhi=zeros(nv,1);				delLambda=zeros(nv,1); 

	theta_rad=[0:0.01:180].*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,love_K_inf'),2); 
	elast_loveH=(loveH(love_num+1).*0.5./sin(0.5.*theta_rad))'...
		+ sum(bsxfun(@times,legendreP,love_H_inf'),2); 
	
	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)); 
		G_rigid(j,:)=0.5./sin(0.5.*alpha); % analytical soln 
		G_elast(j,:)=interp1(xx,elast_loveK,cos(alpha)); 
		U_elast(j,:)=interp1(xx,elast_loveH,cos(alpha)); 
		%display([num2str(j),' done']);
	end 
	% }}} 
		disp(['     done! ']);
	end

	% some additional initialization {{{ 
	lMax=2;
	% lat -> [0,180]; long -> [0,360] to compute spherical harmonics  
	sh=sharmonics(lat_element,lon_element,lMax);  
	sh_vertex=sharmonics(lat,lon,lMax); 
	% SH coefficients of ocean and ice function 
	delH_lm=zeros(9,length(time)); 
	for j=1:length(time)
		delH_lm(:,j) = shlm_element(md.mesh.elements,sh,lMax,ice_wrt0(:,j),areas); 
	end

	oce_area_frac=sum(oce.*areas)./sum(areas); 
	Galpha=G_rigid+G_elast-U_elast;		Ualpha=U_elast; 

	rot_lm=zeros(9,1); 
	delS=zeros(nv,length(time)); 
	delU=zeros(nv,length(time)); 
	delN=zeros(nv,length(time)); 
	eus_test=zeros(nv,length(time)); 
	c1=zeros(1,length(time-1));	c2=zeros(1,length(time-1));	c3=zeros(1,length(time-1));
	xi1=zeros(1,length(time));		xi2=zeros(1,length(time));		xi3=zeros(1,length(time));	
	xp=zeros(1,length(time));		yp=zeros(1,length(time));		
	% }}} 

	for kk=2:length(time)
		disp(['     time # ', num2str(kk), ' :::::::::::::::::::::']);
		%%% step 2- compute sea-level at other geophysical parameters at t=t_k {{{
		%%% in this method: at kk=0, (m1,m2) comes from the previous simulation
		term0=-rho_i*sum(ice_wrt0(:,kk).*areas)/(rho_o*sum(oce.*areas)); 
		term1=3*ice_rat.*sum(bsxfun(@times,Galpha,(ice_wrt0(:,kk).*areas)'),2)./sum(areas); 
		func3=mean(term1(index),2).*oce;
		term3=ones(nv,1).*sum(func3.*areas)./(oce_area_frac*sum(areas)); 
		% 
		delS(:,kk)=term0+term1-term3; % initialization 
		norm_diff = 1.0;		p = 0; 
		while norm_diff > convg_accuracy
			delS_temp=delS(:,kk); 
			norm_old = sqrt(sum(delS_temp.^2)); 
			%
			term2=3*oce_rat.*sum(bsxfun(@times,Galpha,(mean(delS_temp(index),2).*oce.*areas)'),2)./sum(areas); 
			func4=mean(term2(index),2).*oce;
			term4=ones(nv,1).*sum(func4.*areas)./(oce_area_frac*sum(areas)); 
			% 
			%inertia perturbation due to surface load and consequent elastic response of the Earth
			delS_oce=mean(delS_temp(index),2).*oce;
			delS_lm = shlm_element(md.mesh.elements,sh,lMax,delS_oce,areas); 
			% 
			dI_13 = -(1+loveK(3)) * (4*pi/sqrt(15)) * rad_e^4 * (rho_i*delH_lm(8,kk) + rho_o*delS_lm(8));
			dI_23 = -(1+loveK(3)) * (4*pi/sqrt(15)) * rad_e^4 * (rho_i*delH_lm(6,kk) + rho_o*delS_lm(6));
			dI_33 = -(1+loveK(3)) * (8*pi/(3*sqrt(5))) * rad_e^4 * (rho_i*delH_lm(7,kk) + rho_o*delS_lm(7));
			% 
			m1 = c1(kk)-dI_13*(sigma_r/(omega*(moi_C-moi_A)));		
			m2 = c2(kk)-dI_23*(sigma_r/(omega*(moi_C-moi_A)));	
			m3 = c3(kk)-dI_33/((1+4*tideLoveK*(moi_C-moi_A)/(3*fluidLove*moi_C))*moi_C);
			% SH coefficints of rotational potential
			rot_lm(5) = -(2/sqrt(60)) *rad_e^2*omega^2 *m1*m2;                     % l=2,m=-2
			rot_lm(6) = -(1/sqrt(15)) *rad_e^2*omega^2 *m2*(1+m3);                    % l=2,m=-1
			rot_lm(7) = (1/(6*sqrt(5))) *rad_e^2*omega^2 *(m1^2+m2^2-2*m3^2-4*m3); % l=2,m=0
			rot_lm(8) = -(1/sqrt(15)) *rad_e^2*omega^2 *m1*(1+m3);                    % l=2,m=1
			rot_lm(9) = -(1/sqrt(60)) *rad_e^2*omega^2 *(m1^2-m2^2);                 % l=2,m=2
			% 
			term5=(1+tideLoveK-tideLoveH).*sum(bsxfun(@times,sh_vertex,rot_lm'),2)./acc_g;
			func5=mean(term5(index),2).*oce;
			term6=ones(nv,1).*sum(func5.*areas)./(oce_area_frac*sum(areas)); 

			delS(:,kk)=term0+term1+term2-term3-term4+(term5-term6); 
			norm_new  = sqrt(sum(delS(:,kk).^2));
			norm_diff = abs(norm_new-norm_old);
			p = p+1;
			if norm_diff > convg_accuracy
				disp(['          iteration # ', num2str(p), ' :: difference in norm = ', num2str(norm_diff)]);
			else
				xp(kk)=m1;		yp(kk)=-m2;							% pole positions 
				xi1(kk)=dI_13/((1+loveK(3))*(moi_C-moi_A));	% excitation function 
				xi2(kk)=dI_23/((1+loveK(3))*(moi_C-moi_A));	% excitation function 
				xi3(kk)=-dI_33/((1+loveK(3))*moi_C);			% excitation function 
				disp(['          iteration # ', num2str(p), ' :: difference in norm = ', num2str(norm_diff)]);
				disp(['          solution converged! ']);
			end
		end
		%%% here you may compute delU and delN 
		% computation of vertical displacement 
		term1_U=3*ice_rat.*sum(bsxfun(@times,Ualpha,(ice_wrt0(:,kk).*areas)'),2)./sum(areas); 
		term2_U=3*oce_rat.*sum(bsxfun(@times,Ualpha,(mean(delS_temp(index),2).*oce.*areas)'),2)./sum(areas); 
		delU(:,kk)=term1_U+term2_U;
		% computation of geoid change 
		%delN(:,kk)=delS(:,kk)+delU(:,k);  %%% this can be computed anytime 
		% }}} 
		%%% step 2- compute c1,c2,c3 at t=t_k+1 {{{
		if (kk<length(time))
			disp(['          +++++ computing (c1, c2, c3) for the next loading point +++++ ']);
			t_e= yts*(time(kk+1)-time(kk)); % (relative evaluation time) in seconds 
			delS_m1m2=delS(:,kk); % initialization 
			norm_diff = 1.0;		p = 0; 
			while norm_diff > convg_accuracy
				norm_old = sqrt(sum(delS_m1m2.^2)); 
				%
				term2=3*oce_rat.*sum(bsxfun(@times,Galpha,(mean(delS_m1m2(index),2).*oce.*areas)'),2)./sum(areas); 
				func4=mean(term2(index),2).*oce;
				term4=ones(nv,1).*sum(func4.*areas)./(oce_area_frac*sum(areas)); 
				% 
				n1=(dI_13*(cos(sigma_s*t_e)-1) + dI_23*sin(sigma_s*t_e))/sigma_s; 
				n2=(dI_13*sin(sigma_s*t_e) - dI_23*(cos(sigma_s*t_e)-1))/sigma_s; 
				% 
				m1 = cos(sigma_s*t_e)*c1(kk) - sin(sigma_s*t_e)*c2(kk)...
					- (sigma_r/(omega*(moi_C-moi_A))) * (dI_13-(omega+sigma_s)*(cos(sigma_s*t_e)*n1 + sin(sigma_s*t_e)*n2)); 
				m2 = sin(sigma_s*t_e)*c1(kk) + cos(sigma_s*t_e)*c2(kk)...
					- (sigma_r/(omega*(moi_C-moi_A))) * (dI_23-(omega+sigma_s)*(sin(sigma_s*t_e)*n1 - cos(sigma_s*t_e)*n2)); 
				% m3 needs not be (re)computed as it does not depend on time 

				% SH coefficints of rotational potential
				rot_lm(5) = -(2/sqrt(60)) *rad_e^2*omega^2 *m1*m2;                     % l=2,m=-2
				rot_lm(6) = -(1/sqrt(15)) *rad_e^2*omega^2 *m2*(1+m3);                    % l=2,m=-1
				rot_lm(7) = (1/(6*sqrt(5))) *rad_e^2*omega^2 *(m1^2+m2^2-2*m3^2-4*m3); % l=2,m=0
				rot_lm(8) = -(1/sqrt(15)) *rad_e^2*omega^2 *m1*(1+m3);                    % l=2,m=1
				rot_lm(9) = -(1/sqrt(60)) *rad_e^2*omega^2 *(m1^2-m2^2);                 % l=2,m=2
				% 
				term5=(1+tideLoveK-tideLoveH).*sum(bsxfun(@times,sh_vertex,rot_lm'),2)./acc_g;
				func5=mean(term5(index),2).*oce;
				term6=ones(nv,1).*sum(func5.*areas)./(oce_area_frac*sum(areas)); 

				delS_m1m2=term0+term1+term2-term3-term4+(term5-term6); 
				norm_new  = sqrt(sum(delS_m1m2.^2));
				norm_diff = abs(norm_new-norm_old);
				p = p+1;
				if norm_diff > convg_accuracy
					disp(['          iteration # ', num2str(p), ' :: difference in norm = ', num2str(norm_diff)]);
				else
					c1(kk+1)=m1;	c2(kk+1)=m2;	c3(kk+1)=m3;	% constants for the next step 
					disp(['          iteration # ', num2str(p), ' :: difference in norm = ', num2str(norm_diff)]);
					disp(['          solution converged! ']);
				end
			end
			% }}} 
		end
	end

	% compute linear trends in solutions {{{ 
	% compute eustatic values 
	for kk=1:141
		eus_temp_delS=delS(:,kk); 
		eus_test_func=mean(eus_temp_delS(index),2).*oce;
		eust(kk)=sum(eus_test_func.*areas)./(oce_area_frac*sum(areas));
	end 
	% calculate the 'linear' trend:: delS  
	for jj=1:nv 
		delS_trend_0=fit(time',delS(jj,:)','poly1'); 
		delS_trend(jj)=delS_trend_0(1)-delS_trend_0(0); 
	end
	% eustatic of the 'linear' trend 
	eus_test_func=mean(delS_trend(index),2).*oce;
	eust_trend=sum(eus_test_func.*areas)./(oce_area_frac*sum(areas));
	% calculate the 'linear' trend:: delU   
	for jj=1:nv 
		delU_trend_0=fit(time',delU(jj,:)','poly1'); 
		delU_trend(jj)=delU_trend_0(1)-delU_trend_0(0); 
	end
	% monthly ice volume/mass 
	Gt_ice = sum(bsxfun(@times,ice_ant,areas),1).*1e-9;	% Gt of ice 
	% spatial trend of mass evolution (meter of ice) 
	for jj=1:ne
		ice_trend_0=fit(time',ice_wrt0(jj,:)','poly1'); 
		delH_trend(jj)=ice_trend_0(1)-ice_trend_0(0); 
	end
	% compute geocentric motion in mm 
	lMax=1;	% compute degree 1 terms 
	sh=sharmonics(lat_element,lon_element,lMax);  
	% SH coefficients of ocean and ice function 
	delH_lm=zeros(4,length(time)); 
	for j=1:length(time)
		delH_lm(:,j) = shlm_element(md.mesh.elements,sh,lMax,ice_wrt0(:,j),areas); 
	end
	delS_lm=zeros(4,length(time)); 
	for jj=1:length(time) 
		delS_temp=delS(:,jj); 
		delS_oce=mean(delS_temp(index),2).*oce;
		delS_lm(:,jj) = shlm_element(md.mesh.elements,sh,lMax,delS_oce,areas); 
	end 
	load_lm = rho_i*delH_lm + rho_o*delS_lm; % this is kg/m^2 
	coeff = sqrt(3)*(1-(-1.28740059+2*(-0.89858519))/3)/rho_e; 
	XX=1000.*coeff.*load_lm(4,:); % mm cosine term 
	YY=1000.*coeff.*load_lm(2,:); % mm sine term   
	ZZ=1000.*coeff.*load_lm(3,:); % mm 
	% }}} 

	% save results 
	%save('./Results/AIS_GRACE.mat','xi1','xi2','xi3','xp','yp','time','eust',...
	%	'delS','delU','ice_wrt0','delS_trend','delU_trend','delH_trend','eust_trend',...
	%	'Gt_ice','lat','lon','lat_element_0','lon_element_0','XX','YY','ZZ');	

end 
% }}} 
% 

