Index: issm/trunk-jpl/examples/SlrGRACE_NIMS/runme.m
===================================================================
--- issm/trunk-jpl/examples/SlrGRACE_NIMS/runme.m	(revision 25938)
+++ issm/trunk-jpl/examples/SlrGRACE_NIMS/runme.m	(revision 25939)
@@ -2,8 +2,5 @@
 addpath('../Data','../Functions');
 
-steps=[1:8];
-
-try
- % [1:8]
+steps=[8];
 
 if any(steps==1) %{{{
@@ -16,5 +13,8 @@
 	md.mesh=gmshplanet('radius',radius,'resolution',resolution);
 
-	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/SlrGRACE_Mesh md;
@@ -39,14 +39,14 @@
 	for i=1:numrefine,
 
-		md.mask.ocean_levelset=gmtmask(md.mesh.lat,md.mesh.long);
+		ocean_mask=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);	coaste.lat=md.mesh.lat(pos);	coaste.long=md.mesh.long(pos); 
+		pos=find(ocean_mask);	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(j),
 				phi2=coaste.lat/180*pi; lambda2=coaste.long/180*pi;
 				deltaphi=abs(phi2-phi1); deltalambda=abs(lambda2-lambda1);
@@ -61,5 +61,5 @@
 		pos=find(distance<mindistance_coast); distance(pos)=mindistance_coast;
 		
-		pos2=find(md.mask.ocean_levelset~=1 & distance>mindistance_land);
+		pos2=find(ocean_mask~=1 & distance>mindistance_land);
 		distance(pos2)=mindistance_land;
 
@@ -72,5 +72,8 @@
 	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/SlrGRACE_Mesh_refined md;
@@ -127,9 +130,9 @@
 	mask = domain_mask(lat_element_0,lon_element_0,domain); 
 
-	md.slr.deltathickness = mask.*water_load*md.materials.rho_freshwater/md.materials.rho_ice;
+	md.solidearth.surfaceload.icethicknesschange = mask.*water_load*md.materials.rho_freshwater/md.materials.rho_ice;
 
 	save ./Models/SlrGRACE_Loads md;
 
-	plotmodel(md,'data',md.slr.deltathickness,'view',[90 90],'title','Ice height equivalent [m]');
+	plotmodel(md,'data',md.solidearth.surfaceload.icethicknesschange,'view',[90 90],'title','Ice height equivalent [m]');
 end %}}}
 
@@ -140,25 +143,18 @@
 	md.solidearth.lovenumbers = lovenumbers('maxdeg',10000);
 
-	md.mask.ocean_levelset = ones(md.mesh.numberofvertices,1);
-	md.mask.ice_levelset = ones(md.mesh.numberofvertices,1);
-	pos=find(md.slr.deltathickness~=0);
-	md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
-	md.mask.land_levelset = 1-md.mask.ocean_levelset;
-
-	md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
-	md.slr.steric_rate=zeros(md.mesh.numberofvertices,1);
-	md.slr.ocean_area_scaling = 1;
-
-	md.slr.spcthickness=NaN*zeros(md.mesh.numberofvertices,1);
-	md.slr.Ngia=zeros(md.mesh.numberofvertices,1);
-	md.slr.Ugia=zeros(md.mesh.numberofvertices,1);
-
-	md.slr.geodetic=1;
-
-	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;
+	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);
+
+	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);
@@ -175,7 +171,9 @@
 	md = loadmodel('./Models/SlrGRACE_Parameterization');
 
-	md.slr.rigid=1;
-	md.slr.elastic=1;
-	md.slr.rotation=1;
+	md.solidearth.settings.rigid=1;
+	md.solidearth.settings.elastic=1;
+	md.solidearth.settings.rotation=1;
+
+	md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'};
 
 	md.cluster=generic('name',oshostname(),'np',3);
@@ -191,6 +189,6 @@
 	md = loadmodel('./Models/SlrGRACE_Solution');
 
-	sol1 = md.slr.deltathickness*100;						% [cm]
-	sol2 = md.results.SealevelriseSolution.Sealevel*1000;	% [mm]
+	sol1 = md.solidearth.surfaceload.icethicknesschange*100;						% [cm]
+	sol2 = md.results.SealevelriseSolution.SealevelRSL*1000;	% [mm]
 
 	sol_name={'Change in water equivalent height [cm]', 'Relative sea-level [mm]'};
@@ -252,12 +250,12 @@
 
 	disp('Projecting  loads onto the mesh...');
-	time_range = [2005+15/365 2014+350/365];
-	%time_range = 2007 + [15 350]/365;
+	%time_range = [2005+15/365 2014+350/365];
+	time_range = 2007 + [15 350]/365;
 	water_load = grace(md.mesh.elements,md.mesh.lat,md.mesh.long,time_range(1),time_range(2));
 
 	[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;
-	md.slr.deltathickness(ne+1,:)=[1:nt]; % times
+	md.solidearth.surfaceload.icethicknesschange = zeros(ne+1,nt);
+	md.solidearth.surfaceload.icethicknesschange(1:ne,:) = water_load*md.materials.rho_freshwater/md.materials.rho_ice;
+	md.solidearth.surfaceload.icethicknesschange(ne+1,:)=[1:nt]; % times
 
 	md.transient.issmb=0;
@@ -266,4 +264,5 @@
 	md.transient.isthermal=0;
 	md.transient.isslr=1;
+	md.transient.isgia=1;
 
 	md.timestepping.start_time=-1;
@@ -273,4 +272,6 @@
 	md.settings.output_frequency=1;
 
+	md.solidearth.requested_outputs = {'Sealevel','SealevelRSL'};
+	
 	md.cluster=generic('name',oshostname(),'np',3);
 	md.verbose=verbose('111111111');
@@ -285,5 +286,5 @@
 	md = loadmodel('./Models/SlrGRACE_Transient');
 
-	time = md.slr.deltathickness(end,:);
+	time = md.solidearth.surfaceload.icethicknesschange(end,:);
 
 	tide_x = 37+29/60;	% degree N
@@ -291,6 +292,6 @@
 
 	for tt=1:length(time)
-		gmsl(tt) = md.results.TransientSolution(tt).SealevelRSLEustatic*1000;
-		sol = (md.results.TransientSolution(tt+1).Sealevel-md.results.TransientSolution(tt).Sealevel)*1000;
+		gmsl(tt) = md.results.TransientSolution(tt).Bslr*1000;
+		sol = md.results.TransientSolution(tt+1).SealevelRSL*1000;
 		slr_incheon(tt) = griddata(md.mesh.lat,md.mesh.long,sol,tide_x,tide_y);
 	end
