Index: /issm/trunk-jpl/examples/EsaWahr/runme.m
===================================================================
--- /issm/trunk-jpl/examples/EsaWahr/runme.m	(revision 22791)
+++ /issm/trunk-jpl/examples/EsaWahr/runme.m	(revision 22791)
@@ -0,0 +1,159 @@
+
+clear all;
+steps=[4];
+
+if any(steps==1)
+	disp('   Step 1: Mesh creation');
+
+	%Generate initial uniform mesh 
+	md=roundmesh(model,100000,10000);  % Domain radius and element size (meters) 
+
+	save ./Models/EsaWahr.Mesh md;
+	
+	plotmodel (md,'data','mesh');
+
+end
+
+if any(steps==2)
+	disp('   Step 2: Define loads');
+	md = loadmodel('./Models/EsaWahr.Mesh');
+
+	% primer 
+	nv=md.mesh.numberofvertices; 
+	ne=md.mesh.numberofelements; 
+	rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice; 
+
+	% elemental centroids 
+	index=md.mesh.elements;		
+	x_cent=mean(md.mesh.x(index),2); 
+	y_cent=mean(md.mesh.y(index),2); 
+
+	% Define a 20-km radius disc and unload it by 1 meter water height equivalent uniformly  
+	md.esa.deltathickness = zeros(ne,1); 
+	disc_radius=20; % km 
+	rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000;	% radial distance in km 
+	md.esa.deltathickness(rad_dist<=disc_radius) = -1.0*rho_w_i;	% 1 m water withdrawl 
+
+	save ./Models/EsaWahr.Loads md; 
+	
+	plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]');
+
+end
+
+if any(steps==3)
+	disp('   Step 3: Parameterization');
+	md = loadmodel('./Models/EsaWahr.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.ice_levelset = -ones(md.mesh.numberofvertices,1); % -1 = ice loads 
+	md.mask.groundedice_levelset = ones(md.mesh.numberofvertices,1); % 1 = ice is grounnded 
+
+	%% 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='EsaWahr';
+	
+	save ./Models/EsaWahr.Parameterization md; 
+
+end
+
+if any(steps==4)
+	disp('   Step 4: Solve Esa solver');
+	md = loadmodel('./Models/EsaWahr.Parameterization');
+
+	% Request outputs 
+	md.esa.requested_outputs = {'EsaUmotion','EsaXmotion','EsaYmotion'};
+	
+	% Cluster info 
+	md.cluster=generic('name',oshostname(),'np',3); 
+	md.verbose=verbose('111111111');
+
+	% Solve 
+	md=solve(md,'Esa');
+
+	save ./Models/EsaWahr.EsaSolutions md; 
+
+end
+
+return; 
+
+% high-res mesh...
+% plot step. 
+% analytical step. 
+
+	
+	vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] 
+	horz_n = md.results.EsaSolution.EsaNmotion*1000; % [mm] 
+	horz_e = md.results.EsaSolution.EsaEmotion*1000; % [mm] 
+	horz = sqrt(horz_n.^2+horz_e.^2); 
+
+	% grid data from the center 
+	xi=[0:100:60000]; % 100 m resolution 
+	yi=zeros(1,length(xi)); 
+	vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); 
+	horz_track=griddata(md.mesh.x,md.mesh.y,horz,xi,yi,'linear'); 
+
+	return; 
+
+	% figure S1a {{{ 
+	set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); 
+	figure('Position', [100, 100, 800, 600]);
+	plotmodel(md,'data',vert,...
+		'xTickLabel#all',[],'yTickLabel#all',[],...
+		'colormap#all','jet','colorbarcornerposition#all','south',...
+		'expdisp#all','../Exp/RoundDomain.exp',...
+		'gridded#all',1,...
+		'axispos#1',[0.02 0.505 0.47 0.47],...
+		'colorbarpos#1',[0.045,0.55,0.18,0.02],'colorbartitle#1','Vertical [mm]',...
+		'caxis#1',[0 3.5],...
+		'data',horz,...
+		'axispos#2',[0.505 0.505 0.47 0.47],...
+		'colorbarpos',[0.53,0.55,0.18,0.02],'colorbartitle#2','Horizontal [mm]',...
+		'caxis#2',[0 0.5],...
+		'data',horz_n,...
+		'axispos',[0.02 0.02 0.47 0.47],...
+		'colorbarpos',[0.045,0.065,0.18,0.02],'colorbartitle#3','North-south [mm]',...
+		'data',horz_e,...
+		'caxis#3-4',[-0.5 0.5],...
+		'axispos',[0.505 0.02 0.47 0.47],...
+		'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]'); 
+	% }}} 
+	%export_fig('./Fig/Fig_S1a.pdf'); 
+
+	% figure S1b {{{
+	set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); 
+	figure1=figure('Position', [100, 100, 700, 400]);
+	ylabel_1={'0',' ','1','','2','','3',''}; 
+	axes1 = axes('Layer','top','Position',[0.1 0.15 0.8 0.8],...
+		'XTick',[0:10:60],'xlim',[0 60],...
+		'ylim',[0 3.5],'Ytick',[0:0.5:3.5],'yticklabel',ylabel_1); 
+		box(axes1,'on'); hold(axes1,'all'); grid on; 
+		xlabel(axes1,'Radial distance [km]'); 
+		ylabel(axes1,'Displacement [mm]');
+		plot([20 20],[0 3.5],'-.k','linewidth',2,'parent',axes1); 
+		h1=plot(xi/1000,vert_track*917/1000,'color',[0.8500 0.3250 0.0980],'linewidth',3,'parent',axes1); 
+		h2=plot(xi/1000,horz_track*917/1000,'b','linewidth',3,'parent',axes1); 
+		% first box 
+		ag1 = gca;
+		leg1a = legend(ag1,[h1,h2],'Vertical (upward)','Horizontal (away from disc)',1);
+		set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',24); 
+		% 
+	set(gcf,'color','w');
+	% }}} 
+	%export_fig('./Fig/Fig_S1b.pdf'); 
+
+