Index: ../trunk-jpl/examples/Functions/grace.m
===================================================================
--- ../trunk-jpl/examples/Functions/grace.m	(revision 26257)
+++ ../trunk-jpl/examples/Functions/grace.m	(revision 26258)
@@ -1,8 +1,11 @@
 % map grace loads in meters [m] of water equivalent height 
-function water_load = grace(index,lat,lon,tmin,tmax); 
+function water_load = grace(index,lat,lon,tmin,tmax,onvertex); 
 
 	%compute centroids using (lat,lon) data {{{
 	ne = length(index); % number of elements 
+	nv = length(lat); % number of vertices
+	% [lat,lon] \in [-90:90,-180,180]; 
+	lat_vertex_0=lat; long_vertex_0=lon; 
 	% lat -> [0,180]; long -> [0,360] to compute centroids 
 	lat=90-lat; 
 	lon(lon<0)=180+(180+lon(lon<0)); 
@@ -115,10 +118,18 @@
 	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);
+	if onvertex
+		water_load_0=zeros(nv,num_yr);
+	else
+		water_load_0=zeros(ne,num_yr);
+	end
 
 	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);
+		if onvertex
+			water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_vertex_0,lon);
+		else 
+			water_load_0(:,jj) = griddata(lat_grace_plus,lon_grace_plus,load_grace_plus(jj,:),lat_element_0,lon_element);
+		end
 		disp([num2str(jj),' of ',num2str(num_yr),' months done!']); 
 	end 
 
Index: ../trunk-jpl/examples/SlrFarrell/runme.m
===================================================================
--- ../trunk-jpl/examples/SlrFarrell/runme.m	(revision 26257)
+++ ../trunk-jpl/examples/SlrFarrell/runme.m	(revision 26258)
@@ -1,10 +1,7 @@
 clear all;
 
-steps=[1:5];
+steps=[4];
 
-try
- % [1:5]
-
 if any(steps==1) 
 	disp('   Step 1: Global mesh creation');
 
@@ -19,16 +16,16 @@
 	md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3); % attributes in [km]
 
 	for i=1:numrefine,
-		md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
+		ocean_mask_levelset=gmtmask(md.mesh.lat,md.mesh.long);
 
 		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);
+		pos=find(~ocean_mask_levelset);	coaste.lat=md.mesh.lat(pos);	coaste.long=md.mesh.long(pos);
+		pos=find(ocean_mask_levelset);	coasto.lat=md.mesh.lat(pos);	coasto.long=md.mesh.long(pos);
 
 		for j=1:md.mesh.numberofvertices
 			phi1=md.mesh.lat(j)/180*pi; lambda1=md.mesh.long(j)/180*pi;
-			if md.mask.ocean_levelset(j),
+			if ocean_mask_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));
@@ -41,7 +38,7 @@
 		end
 		pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
 
-		pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
+		pos2=find(ocean_mask_levelset~=1 & distance>mindistance_land);
 		distance(pos2)=mindistance_land;
 
 		dist=min(maxdistance,distance);
@@ -48,7 +45,10 @@
 		md.mesh=gmshplanet('radius',radius*1e-3,'resolution',resolution*1e-3,'refine',md.mesh,'refinemetric',dist);
 	end
 
-	md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
+	ocean_mask=gmtmask(md.mesh.lat,md.mesh.long);
+	pos = find(ocean_mask==0); 
+	md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); 
+	md.mask.ocean_levelset(pos)=1; 
 
 	save ./Models/SlrFarrell_Mesh md;
 
@@ -59,13 +59,19 @@
 	disp('   Step 2: Define source as in Farrell, 1972, Figure 1');
 	md = loadmodel('./Models/SlrFarrell_Mesh');
 
-	md.slr.sealevel=md.mask.ocean_levelset; % 1 m SLR everywhere
-	md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
-	md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
+	md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofvertices,1);
+	md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); 
 
+	pos = find(md.mask.ocean_levelset==-1); 
+	md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
+	md.solidearth.initialsealevel(pos)=1; % 1 m SLR everywhere
+	md.dsl.global_average_thermosteric_sea_level_change=[0;0];
+	md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
+	md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
+
 	save ./Models/SlrFarrell_Loads md;
 
-	plotmodel (md,'data',md.slr.sealevel,'view',[90 90],'title#all','Initial sea-level [m]');
+	plotmodel (md,'data',md.solidearth.initialsealevel,'view',[90 90],'title#all','Initial sea-level [m]');
 end 
 
 if any(steps==3) 
@@ -74,18 +80,17 @@
 
 	md.solidearth.lovenumbers=lovenumbers('maxdeg',10000);
 
-	md.mask.land_levelset = 1-md.mask.ocean_levelset;
-	md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
-	md.mask.ocean_levelset = -ones(md.mesh.numberofvertices,1);
-	pos=find(md.mesh.lat <-80);
-	md.mask.ice_levelset(pos(1))=-1; % ice yes!
-	md.mask.ocean_levelset(pos(1))=1; % ice grounded!
+	%md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
+	%md.mask.ocean_levelset = -ones(md.mesh.numberofvertices,1);
+	%pos=find(md.mesh.lat <-80);
+	%md.mask.ice_levelset(pos(1))=-1; % ice yes!
+	%md.mask.ocean_levelset(pos(1))=1; % ice grounded!
 
-	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;
+	% arbitary to pass consistency check. 
+	md.geometry.bed=-ones(md.mesh.numberofvertices,1);
+	md.geometry.surface=ones(md.mesh.numberofvertices,1);
+	md.geometry.base=md.geometry.bed; 
+	md.geometry.thickness=md.geometry.surface-md.geometry.base; 
 
 	md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
 	md.materials.rheology_B=paterson(md.initialization.temperature);
@@ -103,7 +108,7 @@
 	md.cluster=generic('name',oshostname(),'np',3);
 	md.verbose=verbose('111111111');
 
-	md.slr.reltol = 0.1/100; % percent change in solution
+	md.solidearth.settings.reltol = 0.1/100; % percent change in solution
 
 	md=solve(md,'Slr');
 
Index: ../trunk-jpl/examples/EsaGRACE/runme.m
===================================================================
--- ../trunk-jpl/examples/EsaGRACE/runme.m	(revision 26257)
+++ ../trunk-jpl/examples/EsaGRACE/runme.m	(revision 26258)
@@ -1,7 +1,7 @@
 clear all;
 addpath('../Data','../Functions');
 
-steps=[1]; % [1:5]
+steps=[5]; 
 
 if any(steps==1) 
 	disp('   Step 1: Global mesh creation');
Index: ../trunk-jpl/examples/EsaWahr/runme.m
===================================================================
--- ../trunk-jpl/examples/EsaWahr/runme.m	(revision 26257)
+++ ../trunk-jpl/examples/EsaWahr/runme.m	(revision 26258)
@@ -1,7 +1,7 @@
 clear all;
 addpath('../Functions');
 
-steps=[1]; % [1:7]
+steps=[2]; 
 
 if any(steps==1) 
 	disp('   Step 1: Mesh creation');
Index: ../trunk-jpl/examples/SlrGRACE/runme.m
===================================================================
--- ../trunk-jpl/examples/SlrGRACE/runme.m	(revision 26257)
+++ ../trunk-jpl/examples/SlrGRACE/runme.m	(revision 26258)
@@ -1,9 +1,9 @@
 clear all;
 addpath('../Data','../Functions');
 
-steps=[7];
+steps=[1:7];
 
-if any(steps==1)
+if any(steps==1) % {{{ 
 	disp('   Step 1: Global mesh creation');
 
 	numrefine=5;
@@ -50,8 +50,8 @@
 	save ./Models/SlrGRACE_Mesh md;
 
 	plotmodel (md,'data',md.mask.ocean_levelset,'edgecolor','k');
-end
-if any(steps==2)
+end % }}} 
+if any(steps==2) % {{{ 
 	disp('   Step 2: Define loads in meters of ice height equivalent');
 	md = loadmodel('./Models/SlrGRACE_Mesh');
 
@@ -58,40 +58,59 @@
 	year_month = 2007+15/365;
 	time_range = [year_month year_month];
 
-	water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
+	onvertex = 1; % map data on vertex. If 0, it maps on the elemental centroid. 
+	water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2),onvertex);
+	rho_water2ice = md.materials.rho_freshwater/md.materials.rho_ice; 
+	ice_load = water_load*rho_water2ice; % ice height equivalent. 
 
-	md.solidearth.surfaceload.icethicknesschange = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
+	%Geometry for the bed, arbitrary thickness of 100:
+	md.geometry.bed=zeros(md.mesh.numberofvertices,1);
+	md.geometry.base=md.geometry.bed;
+	md.geometry.thickness = 100*ones(md.mesh.numberofvertices,1);
+	md.geometry.surface=md.geometry.bed+md.geometry.thickness;
 
+	md.masstransport.spcthickness = [md.geometry.thickness + ice_load; 0];
+
+	md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); 
+
 	save ./Models/SlrGRACE_Loads md;
 
-	plotmodel (md,'data',md.solidearth.surfaceload.icethicknesschange,'view',[90 -90],'caxis',[-.1 .1],'title','Ice height equivalent [m]');
-end
-if any(steps==3)
+	plotmodel (md,'data',ice_load,...
+		'edgecolor','k','view',[45 45],'caxis',[-.1 .1],...
+		'title','Ice height equivalent [m]');
+end % }}} 
+if any(steps==3) % {{{ 
 	disp('   Step 3: Parameterization');
 	md = loadmodel('./Models/SlrGRACE_Loads');
-
+	
+	md.mask.ice_levelset=-md.mask.ocean_levelset;
+	
 	md.solidearth.lovenumbers = lovenumbers('maxdeg',10000);
+	md.solidearth.settings.reltol=NaN;
+	md.solidearth.settings.abstol=1e-3;
+	md.solidearth.settings.computesealevelchange=1;
+	md.solidearth.settings.isgrd=1;
+	md.solidearth.settings.grdmodel=1;
+	md.solidearth.settings.maxiter=10; 
+	
+	%time stepping:
+	md.timestepping.start_time=0;
+	md.timestepping.time_step=1;
+	md.timestepping.final_time=1;
 
-	md.mask.ice_levelset=-md.mask.ocean_levelset;
-	md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
-	md.dsl.global_average_thermosteric_sea_level_change=[0;0];
-	md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
-	md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
+	%masstransport:
+	md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+	md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+	md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+	md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+	md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
+	md.initialization.str=0;
 
-	md.solidearth.settings.ocean_area_scaling=1;
-
-	% arbitary to pass consistency check.	md.geometry.bed=-ones(md.mesh.numberofvertices,1);
-	md.geometry.surface=ones(md.mesh.numberofvertices,1);
-	md.geometry.base=md.geometry.bed;	md.geometry.thickness=md.geometry.surface-md.geometry.base;
-	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);
-
 	md.miscellaneous.name='SlrGRACE';
 
 	save ./Models/SlrGRACE_Parameterization md;
-end
-if any(steps==4)
+end % }}} 
+if any(steps==4) % {{{ 
 	disp('   Step 4: Solve Slr solver');
 	md = loadmodel('./Models/SlrGRACE_Parameterization');
 
@@ -99,21 +118,28 @@
 	md.solidearth.settings.elastic=1;
 	md.solidearth.settings.rotation=1;
 
-	md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'};
+	%Physics:
+	md.transient.issmb=0;
+	md.transient.isstressbalance=0;
+	md.transient.isthermal=0;
+	md.transient.ismasstransport=1;
+	md.transient.isslc=1;
+	
+	md.solidearth.requested_outputs={'Sealevel','Bed'}; 
 
 	md.cluster=generic('name',oshostname(),'np',3);
 	md.verbose=verbose('111111111');
 
-	md=solve(md,'Slr');
+	md=solve(md,'Transient');
 
 	save ./Models/SlrGRACE_Solution md;
-end
-if any(steps==5)
+end % }}}
+if any(steps==5) % {{{ 
 	disp('   Step 5: Plot solutions');
 	md = loadmodel('./Models/SlrGRACE_Solution');
 
-	sol1 = md.solidearth.surfaceload.icethicknesschange*100; % [cm]
-	sol2 = md.results.SealevelriseSolution.SealevelRSL*1000;	% [mm]
+	sol1 = (md.masstransport.spcthickness(1:end-1)-md.geometry.thickness)*100; % [cm]
+	sol2 = (md.results.TransientSolution.Sealevel-md.results.TransientSolution.Bed)*1000;	% [mm]
 
 	sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'};
 	fig_name={'Fig_dH.pdf','Fig_slr.pdf'};
@@ -153,7 +179,7 @@
 		if (kk==1)
 			geoshow(flipud(coastlat),flipud(coastlon),'DisplayType','polygon','FaceColor','white');
 		else
-			geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor',[.5 1 .5]);
+			geoshow(coastlat,coastlon,'DisplayType','polygon','FaceColor','none');
 		end
 		plot(coastlon, coastlat,'k'); hold off;
 		c1=colorbar;
@@ -166,8 +192,8 @@
 		set(gcf,'color','w');
 		%export_fig(fig_name{kk});
 	end
-end
-if any(steps==6)
+end % }}} 
+if any(steps==6) % {{{
 	disp('   Step 6: Transient run');
 	md = loadmodel('./Models/SlrGRACE_Parameterization');
 
@@ -200,8 +226,8 @@
 	md=solve(md,'tr');
 
 	save ./Models/SlrGRACE_Transient md;
-end
-if any(steps==7)
+end % }}} 
+if any(steps==7) % {{{ 
 	disp('   Step 7: Plot transient');
 	md = loadmodel('./Models/SlrGRACE_Transient');
 
@@ -283,4 +309,5 @@
 	ylabel('GMSL [mm/yr]');
 	set(gcf,'color','w');
 	%export_fig('Fig7.pdf');
-end
+end % }}} 
+