Index: /issm/trunk-jpl/examples/EsaGRACE/grace.m
===================================================================
--- /issm/trunk-jpl/examples/EsaGRACE/grace.m	(revision 22796)
+++ /issm/trunk-jpl/examples/EsaGRACE/grace.m	(revision 22796)
@@ -0,0 +1,121 @@
+% map grace loads in meters [m] of water equivalent height 
+function water_load = grace(index,lat,lon,tmin,tmax); 
+
+	%compute centroids using (lat,lon) data {{{
+	ne = length(index); % number of elements 
+	% lat -> [0,180]; long -> [0,360] to compute centroids 
+	lat=90-lat; 
+	lon(lon<0)=180+(180+lon(lon<0)); 
+	
+	ax_0=lat(index(:,1)); ay_0=lon(index(:,1)); 
+	bx_0=lat(index(:,2)); by_0=lon(index(:,2)); 
+	cx_0=lat(index(:,3)); cy_0=lon(index(:,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;
+
+	% [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;
+	% }}}
+
+	% Monthly GRACE data 
+	filename=['GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc']; 
+	time_0=ncread(filename,'time'); % days since 2002-01-01 00:00:00 UTC 
+	long_0=ncread(filename,'lon'); % longitudes: 0.27-359.75
+	lati_0=ncread(filename,'lat'); % latitudes: -89.75:89.75
+	rec=ncread(filename,'lwe_thickness');% rec_ensemble_mean [cm]
+
+	time_yr = 2002+time_0/365; % [yr] 
+
+	[nn_0,mm_0] = size(squeeze(rec(:,:,1))); 
+	for jj=1:mm_0     % chose a latitude
+		for kk=1:nn_0
+			ii=nn_0*(jj-1)+kk;
+			lat_0(ii)=lati_0(jj); 
+			lon_0(ii)=long_0(kk); 
+			tws_monthly(:,ii) = rec(kk,jj,:);
+		end 
+	end 
+	% 
+	%% translate variables as grace-related variables -- so I do not need to do too much editing 
+	grace_monthly=tws_monthly; 
+	grace_monthly(grace_monthly<-1344.6016 | grace_monthly>858.5046)=0; 
+
+	% detrend over the entire time period 
+	grace_monthly = detrend(grace_monthly); % 159X64800 => remove trends from each column! 
+
+	% fill out the blanks {{{ 
+	
+	lat_grace=lat_0; 
+	lon_grace=lon_0; 
+	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)=grace_monthly(:,mm); 
+			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)=grace_monthly(:,jj); 
+			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(length(time_0),num_plus); 
+	
+	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; 
+	load_grace_plus(:,1+num_org:num_plus)=load_new; %(:,:);
+	% }}}
+
+	% choose selected months ONLY 
+	[diff1,pos1] = min(abs(tmin-time_yr));
+	[diff2,pos2] = min(abs(tmax-time_yr)); 
+
+	time_yr=time_yr(pos1:pos2); 
+	load_grace_plus=load_grace_plus(pos1:pos2,:); 
+	num_yr=length(time_yr); 
+	water_load_0=zeros(ne,num_yr);
+
+	for jj=1:num_yr
+		water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element);
+		disp([num2str(jj),' of ',num2str(num_yr),' months done!']); 
+	end 
+
+	water_load=water_load_0/100;		% cm -> meters of water 
+	water_load(isnan(water_load))=0; 
+
Index: /issm/trunk-jpl/examples/EsaGRACE/runme.m
===================================================================
--- /issm/trunk-jpl/examples/EsaGRACE/runme.m	(revision 22796)
+++ /issm/trunk-jpl/examples/EsaGRACE/runme.m	(revision 22796)
@@ -0,0 +1,171 @@
+
+clear all;
+steps=[5]; % [1:6]; 
+
+if any(steps==0) % Download GRACE land_mass data {{{
+	disp('   Step 0: Download GRACE land_mass data');
+
+	% Connect to ftp server and download 
+	f = ftp('podaac-ftp.jpl.nasa.gov');
+	%dir(f)
+	cd(f,'allData/tellus/L3/land_mass/RL05/netcdf/');
+   mget(f,'GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc');
+
+	% display the content of the netcdf file. 
+	%ncdisp('GRCTellus.JPL.200204_201701.LND.RL05_1.DSTvSCS1411.nc')
+
+end % }}}
+
+if any(steps==1) % Global mesh creation {{{ 
+	disp('   Step 1: Global mesh creation');
+
+	resolution=300;	% km. uniform meshing... 
+	radius = 6.371012*10^3;	% mean radius of Earth, km
+
+	%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);
+
+	%figure out mask: 
+	md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 
+
+	save ./Models/EsaGRACE.Mesh md;
+	
+	plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
+	%export_fig('Fig1.pdf'); 
+
+end % }}} 
+
+if any(steps==2) % Define loads {{{ 
+	disp('   Step 2: Define loads in meters of ice height equivalent');
+	md = loadmodel('./Models/EsaGRACE.Mesh');
+
+	year_month = 2007+15/365;
+	time_min=year_month; 
+	time_max=year_month; 
+
+	% map GRACE water load on to the mesh for the seleted month(s) 
+	water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_min,time_max); 
+	
+	md.esa.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent 
+
+	save ./Models/EsaGRACE.Loads md; 
+	
+	plotmodel (md,'data',md.esa.deltathickness,...
+		'view',[90 -90],'caxis',[-.1 .1],...
+		'title','Ice height equivalent [m]');
+	%export_fig('Fig2.pdf'); 
+
+end % }}} 
+
+if any(steps==3) % Parameterization {{{ 
+	disp('   Step 3: Parameterization');
+	md = loadmodel('./Models/EsaGRACE.Loads');
+
+	% Love numbers and reference frame: CF or CM (choose one!)  
+	nlove=10001;	% up to 10,000 degree 
+	md.esa.love_h = love_numbers('h','CF'); md.esa.love_h(nlove+1:end)=[];
+	md.esa.love_l = love_numbers('l','CF'); md.esa.love_l(nlove+1:end)=[];
+
+	% Mask: for computational efficiency only those elements that have loads are convolved! 
+	md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded 
+	md.mask.ice_levelset = ones(md.mesh.numberofvertices,1); 
+	pos=find(md.esa.deltathickness~=0);
+	md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads 
+	md.mask.land_levelset = 1-md.mask.ocean_levelset; 
+
+	%% IGNORE BUT DO NOT DELETE %% {{{
+	% Geometry: Important only when you want to couple with Ice Flow Model 
+	di=md.materials.rho_ice/md.materials.rho_water;
+	md.geometry.thickness=ones(md.mesh.numberofvertices,1);
+	md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
+	md.geometry.base=md.geometry.surface-md.geometry.thickness;
+	md.geometry.bed=md.geometry.base;
+	% Materials: 
+	md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
+	md.materials.rheology_B=paterson(md.initialization.temperature);
+	md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+	% Miscellaneous: 
+	md.miscellaneous.name='EsaGRACE';
+	%% IGNORE BUT DO NOT DELETE %% }}}  
+	
+	save ./Models/EsaGRACE.Parameterization md; 
+
+end % }}} 
+
+if any(steps==4) % Solve {{{ 
+	disp('   Step 4: Solve Esa solver');
+	md = loadmodel('./Models/EsaGRACE.Parameterization');
+
+	% Request outputs 
+	md.esa.requested_outputs = {'EsaUmotion','EsaNmotion','EsaEmotion'};
+	
+	% Cluster info 
+	md.cluster=generic('name',oshostname(),'np',3); 
+	md.verbose=verbose('111111111');
+
+	% Solve 
+	md=solve(md,'Esa');
+
+	save ./Models/EsaGRACE.Solution md; 
+
+end % }}} 
+
+if any(steps==5) % Plot solutions {{{ 
+	disp('   Step 5: Plot solutions');
+	md = loadmodel('./Models/EsaGRACE.Solution');
+
+	% solutions. 
+	ur = md.results.EsaSolution.EsaUmotion*1000; % [mm] 
+	un = md.results.EsaSolution.EsaNmotion*1000; % [mm] 
+	ue = md.results.EsaSolution.EsaEmotion*1000; % [mm] 
+
+	res = 1.0; % degree 
+
+	% Make a grid of lats and lons, based on the min and max of the original vectors
+	[lat_grid, lon_grid] = meshgrid(linspace(-90,90,180/res), linspace(-180,180,360/res));
+	sol_grid = zeros(size(lat_grid)); 
+
+	sol = ue; 
+
+	% Make a interpolation object
+	F = scatteredInterpolant(md.mesh.lat,md.mesh.long,sol); 
+	F.Method = 'linear';
+	F.ExtrapolationMethod = 'linear'; 
+
+	% Do the interpolation to get gridded solutions... 
+	sol_grid = F(lat_grid, lon_grid);
+	sol_grid(isnan(sol_grid))=0; 
+
+	% set polar unphysical 0s to Nan 
+	sol_grid(lat_grid>85 & sol_grid==0) =NaN; 
+
+	set(0,'DefaultAxesFontSize',18,'DefaultAxesLineWidth',1,'DefaultTextFontSize',18,'DefaultLineMarkerSize',8)
+	figure1=figure('Position', [100, 100, 1000, 500]); 
+	gcf; 
+	load coast; 
+	cla; 
+	pcolor(lon_grid,lat_grid,sol_grid); shading flat; 
+	%caxis([-0.3 0.3])
+	hold on
+	plot(long,lat,'k');
+	%geoshow(lat,long,'DisplayType','polygon','FaceColor','white'); 	
+	hold off; 
+	c1=colorbar;
+	colormap(jet); 
+	xlim([-180 180]); 
+	ylim([-90 90]); 
+	grid on; 
+	title(['Average change in relative sea-level [mm/yr]']);
+	set(gcf,'color','w');
+
+	%export_fig('Fig5.pdf'); 
+
+end % }}} 
+
+if any(steps==6) % {{{ Transient 
+	disp('   Step 6: Transient run');
+
+end % }}} 
+