Index: /issm/trunk-jpl/examples/SlrGRACE/runme.m
===================================================================
--- /issm/trunk-jpl/examples/SlrGRACE/runme.m	(revision 22804)
+++ /issm/trunk-jpl/examples/SlrGRACE/runme.m	(revision 22804)
@@ -0,0 +1,274 @@
+
+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');
+
+	numrefine=1;
+	resolution=150*1e3;		% inital resolution [m]. It determines, e.g., whether we capture small islands. 
+	radius = 6.371012*10^6;	% mean radius of Earth, m
+	mindistance_coast=150*1e3;		% coastal resolution [m] 
+	mindistance_land=300*1e3; % resolution on the continents [m]
+	maxdistance=600*1e3;	 % max element size (on mid-oceans) [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*1e-3,'resolution',resolution*1e-3); % attributes should be in km. 
+
+	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*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*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_coast); distance(pos)=mindistance_coast;
+		
+		% refine on the continents
+		pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land); 
+		distance(pos2)=mindistance_land; 
+
+		dist=min(maxdistance,distance); % max size 1000 km 
+		%use distance to the coastline to refine mesh: 
+		md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist);
+	end
+
+	%figure out mask: 
+	md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long); 
+
+	save ./Models/SlrGRACE.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/SlrGRACE.Mesh');
+
+	% define time interval for analysis 
+	year_month = 2007+15/365;
+	time_range = [year_month 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_range(1),time_range(2)); 
+	
+	md.slr.deltathickness = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent 
+
+	save ./Models/SlrGRACE.Loads md; 
+	
+	plotmodel (md,'data',md.slr.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/SlrGRACE.Loads');
+
+	% Love numbers and reference frame: CF or CM (choose one!)  
+	nlove=10001;	% up to 10,000 degree 
+	md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlove+1:end)=[];
+	md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlove+1:end)=[];
+	md.slr.love_l = love_numbers('l','CM'); md.slr.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.slr.deltathickness~=0);
+	md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; % -1 = ice loads 
+	md.mask.land_levelset = 1-md.mask.ocean_levelset; 
+
+	% initialize 
+	md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
+	% steric rate
+	md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
+	% solution is scaled according to "exact" area of the oceans 
+	md.slr.ocean_area_scaling = 1;
+
+	% convergence criteria 
+	%md.slr.reltol=NaN;
+	%md.slr.abstol=1e-4;
+
+	%% 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='SlrGRACE';
+	%% IGNORE BUT DO NOT DELETE %% }}}  
+	
+	save ./Models/SlrGRACE.Parameterization md; 
+
+end % }}} 
+
+if any(steps==4) % Solve {{{ 
+	disp('   Step 4: Solve Slr solver');
+	md = loadmodel('./Models/SlrGRACE.Parameterization');
+
+	% What physics to capture? 
+	md.slr.rigid=1; 
+	md.slr.elastic=1; 
+	md.slr.rotation=1;
+
+	% Request outputs 
+%	md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'};
+	
+	% Cluster info 
+	md.cluster=generic('name',oshostname(),'np',3); 
+	md.verbose=verbose('111111111');
+
+	% Solve 
+	md=solve(md,'Slr');
+
+	save ./Models/SlrGRACE.Solution md; 
+
+end % }}} 
+
+if any(steps==5) % Plot solutions {{{ 
+	disp('   Step 5: Plot solutions');
+	md = loadmodel('./Models/SlrGRACE.Solution');
+
+	% loads and solutions. 
+	sol1 = md.slr.deltathickness*100; % WEH cm 
+	sol2 = md.results.SealevelriseSolution.Sealevel*1000; % [mm] 
+	sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [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)); 
+
+	for kk=1:2
+		sol=eval(sprintf('sol%d',kk));
+	
+		% if data are on elements, map those on to the vertices {{{
+		if length(sol)==md.mesh.numberofelements 
+			% map on to the vertices 
+			for jj=1:md.mesh.numberofelements
+				ii=(jj-1)*3;
+				pp(ii+1:ii+3)=md.mesh.elements(jj,:);
+			end
+			for jj=1:md.mesh.numberofvertices
+				pos=ceil(find(pp==jj)/3); 
+				temp(jj)=mean(sol(pos)); 
+			end
+			sol=temp'; 
+		end % }}}
+
+		% 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; 
+		sol_grid(lat_grid>85 & sol_grid==0) =NaN; % set polar unphysical 0s to 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; hold on; 
+		plot(long,lat,'k'); hold off; 
+		c1=colorbar;
+		colormap(jet); 
+		xlim([-180 180]); 
+		ylim([-90 90]); 
+		grid on; 
+		title(sol_name(kk)); 
+		set(gcf,'color','w');
+
+		%export_fig('Fig5.pdf'); 
+	end
+
+end % }}} 
+
+if any(steps==6) % Transient {{{ 
+	disp('   Step 6: Transient run');
+	md = loadmodel('./Models/SlrGRACE.Parameterization');
+
+	% read GRACE data during the chosen time period << steps=2 >>
+	disp('Projecting grace loads onto the mesh...');
+	time_range = 2007 + [15 350]/365; 
+	water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2)); 
+
+	% define ice load
+	[ne,nt]=size(water_load); 
+   md.slr.deltathickness = zeros(ne+1,nt);
+	md.slr.deltathickness(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice; % ice height equivalent 
+	md.slr.deltathickness(ne+1,:)=[1:nt]; % times 
+	
+	% define transient run 
+	md.transient.issmb=0;
+	md.transient.ismasstransport=0; 
+	md.transient.isstressbalance=0; 
+	md.transient.isthermal=0;
+	md.transient.isslr=1;
+	
+	% time stepping, output frequency etc. 
+	md.timestepping.start_time=-1;
+   md.timestepping.final_time=nt; 
+	md.timestepping.time_step=1; 
+   md.timestepping.interp_forcings=0;
+	md.settings.output_frequency=1;
+
+	% Request outputs 
+	md.slr.requested_outputs = {'SlrUmotion','SlrNmotion','SlrEmotion'};
+	
+	% Cluster info 
+	md.cluster=generic('name',oshostname(),'np',3); 
+	md.verbose=verbose('111111111');
+
+	% Solve 
+	md=solve(md,'tr');
+
+end % }}} 
+