Index: issm/trunk-jpl/test/SandBox/test2004.m
===================================================================
--- issm/trunk-jpl/test/SandBox/test2004.m	(revision 24863)
+++ issm/trunk-jpl/test/SandBox/test2004.m	(revision 24905)
@@ -1,5 +1,5 @@
 %Test Name: Earth_Antarctica_GIA
 		
-testagainst2002=1;
+testagainst2002=0;
 
 %Data paths: {{{
@@ -29,8 +29,18 @@
 	}));
 	%}}}
+	%Ross: {{{
+	sl.addbasin(basin('continent','antarctica','name','ross','proj',proj3031,'boundaries',{...
+		boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
+		boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031),...
+		boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),...
+		boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031,'orientation','reverse')...
+		}));
+%}}}
 	%HemisphereEast: {{{
 	sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','proj',laea(0,+90),'boundaries',{... %Australian projection lat,long
 		boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326),...
 		boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
+		boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031),...
+		boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),...
 		boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031,'orientation','reverse'),...
 		boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031)...
@@ -41,4 +51,6 @@
 	boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031),...
 	boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031),...
+	boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),...
+	boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031,'orientation','reverse'),...
 	boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),...
 	boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031)...
@@ -61,8 +73,8 @@
 %Go through basins and mesh:  %{{{
 %meshing parameters:  {{{
-hmin=2000; hmax=10000; hmin=hmin*1000; hmax=hmax*1000; 
+hmin=500; hmax=700; hmin=hmin*1000; hmax=hmax*1000; 
 tolerance=100; %tolerance of 100m on Earth position when mergin 3d meshes
 threshold=5;
-defaultoptions={'KeepVertices',0,'MaxCornerAngle',0.0000000001,'NoBoundaryRefinment',1}; 
+defaultoptions={'KeepVertices',0,'MaxCornerAngle',0.0000000001,'NoBoundaryRefinement',1}; 
 alreadyloaded=0;
 %}}}
@@ -104,34 +116,12 @@
 	md=sl.icecaps{ind}; bas=sl.basins{ind}; 
 	%masks :  %{{{
-	%new type of mask: 
-	md.mask=maskpsl; 
-
-	%land and ocean: 
-	md.mask.land_levelset=ones(md.mesh.numberofvertices,1);
-	md.mask.ocean_levelset=zeros(md.mesh.numberofvertices,1);
-
 	%ice levelset from domain outlines: 
 	md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1);
-	pos=find(md.mesh.vertexonboundary); md.mask.ice_levelset(pos)=0;
-
-	%ice front: 
-	pos=find(md.mask.ice_levelset==0); md.mask.ocean_levelset(pos)=-1; md.mask.land_levelset(pos)=1;
-
-	%now, glaciers: 
-	md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1);
-	%}}} 
-	%initial grounding line position: {{{
+	
 	if bas.isnameany('antarctica-grounded'), 
-
 		md.mask.ocean_levelset=ones(md.mesh.numberofvertices,1);
 	end
-	if bas.isnameany('ronne'), 
-
+	if bas.isnameany('ronne','ross'), 
 		md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
-
-		%correction to land and ocean levelset: ice shelf is not on land! 
-		pos=find(md.mask.ice_levelset<=0 & md.mask.ocean_levelset<=0);
-		md.mask.ocean_levelset(pos)=1;
-		md.mask.land_levelset(pos)=-1;
 	end
 	%}}}
@@ -153,20 +143,15 @@
 			late=sum(md.mesh.lat(md.mesh.elements),2)/3;
 			longe=sum(md.mesh.long(md.mesh.elements),2)/3;
-			pos=find(late <-80);
-			md.slr.deltathickness(pos)=-100;
+			pos=find(late <-85);
+			ratio=0.225314032985172/0.193045366574523;
+			%ratio=   1.276564103522540/.869956;
+			md.slr.deltathickness(pos)=-100*ratio;
 		else
-			delH=textread([modeldatapath '/AdhikariSciAdvTrends/AIS_delH_trend.txt']);
-			longAIS=delH(:,1);
-			latAIS=delH(:,2);
-			delHAIS=delH(:,3);
-			index=delaunay(latAIS,longAIS);
-
-			lat=md.mesh.lat; 
-			long=md.mesh.long+360; 
-			pos=find(long>360);long(pos)=long(pos)-360;
-
+			delH=textread('../Data/AIS_delH_trend.txt');
+			longAIS=delH(:,1); latAIS=delH(:,2); delHAIS=delH(:,3); index=delaunay(longAIS,latAIS);
+			lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360;
 			delHAIS=InterpFromMesh2d(index,longAIS,latAIS,delHAIS,long,lat);
 			northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0;
-
+			
 			md.slr.deltathickness=delHAIS(md.mesh.elements)*[1;1;1]/3/100;
 		end
@@ -205,6 +190,4 @@
 
 	%mask:  %{{{
-	%new type of mask: 
-	md.mask=maskpsl;
 	%Figure out mask from initial mesh: deal with land and ocean masks (land include grounded ice).  %{{{
 	%first, transform land element  mask into vertex driven one.  
@@ -226,12 +209,7 @@
 	ind2=[md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3); md.mesh.elements(landelsconocean,1)];
 
-	h=waitbar(0,'Starting mask processing');
 
 	%edge ind1 and ind2: 
 	for i=1:length(ind1),
-		if ~mod(i,100),
-			waitbar(i/length(ind1),h,sprintf('%.2i%% along...',i/length(ind1)*100));
-		end
-
 		els1=md.mesh.vertexconnectivity(ind1(i),1: md.mesh.vertexconnectivity(ind1(i),end));
 		els2=md.mesh.vertexconnectivity(ind2(i),1: md.mesh.vertexconnectivity(ind2(i),end));
@@ -244,18 +222,15 @@
 		end
 	end
-	md.mask.land_levelset=land_mask;
-	close(h);
-	%}}}
-	%work on ocean, glaciers and ice:  {{{
-	%ocean: opposite of land: 
-	md.mask.ocean_levelset=-md.mask.land_levelset;
-
-	%initliaze: 
-	md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); 
-
-	%grounded ice: 
-	md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); 
-
-	md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1);
+
+	md.mask.ocean_levelset=land_mask;
+	md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);   %if there are glaciers, we'll adjust 
+
+	if testagainst2002,
+		% {{{
+		%greenland
+		pos=find(md.mesh.lat > 70 &  md.mesh.lat < 80 & md.mesh.long>-60 & md.mesh.long<-30);
+		md.mask.ice_levelset(pos)=-1;
+		% }}}
+	end
 
 	% }}}
@@ -263,5 +238,6 @@
 	%slr loading/calibration:  {{{
 
-	if testagainst2002,
+	if testagainst2002, 
+		% {{{
 		md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
 		%greenland
@@ -269,10 +245,38 @@
 		longe=sum(md.mesh.long(md.mesh.elements),2)/3;
 		pos=find(late > 70 &  late < 80 & longe>-60 & longe<-30);
-		md.slr.deltathickness(pos)=-100;
+		ratio=.3823/.262344;
+		%md.slr.deltathickness(pos)=-100*ratio;
 
 		%correct mask: 
-		md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
+		md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; 
+		% }}}
 	else
+
 		md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
+
+		delH=textread('../Data/GIS_delH_trend.txt');
+		longGIS=delH(:,1); latGIS=delH(:,2); delHGIS=delH(:,3); index=delaunay(longGIS,latGIS);
+		lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360;
+		delHGIS=InterpFromMeshToMesh2d(index,longGIS,latGIS,delHGIS,long,lat);
+		delHGISe=delHGIS(md.mesh.elements)*[1;1;1]/3;
+	
+		delH=textread('../Data/GLA_delH_trend_15regions.txt');
+		longGLA=delH(:,1); latGLA=delH(:,2); delHGLA=sum(delH(:,3:end),2); index=delaunay(longGLA,latGLA);
+		lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360;
+		delHGLA=InterpFromMeshToMesh2d(index,longGLA,latGLA,delHGLA,long,lat);
+		delHGLAe=delHGLA(md.mesh.elements)*[1;1;1]/3;
+
+		pos=find(delHGISe);
+		md.slr.deltathickness(pos)=delHGISe(pos)/100;
+		pos=find(delHGLAe);
+		md.slr.deltathickness(pos)=delHGLAe(pos)/100;
+
+		%adjust mask accordingly: 
+		pos=find(md.slr.deltathickness); 
+		flags=zeros(md.mesh.numberofvertices,1); 
+		flags(md.mesh.elements(pos,:))=1;
+		pos=find(flags);
+		md.mask.ice_levelset(pos)=-1;
+		md.mask.ocean_levelset(pos)=1;
 	end
 
@@ -288,6 +292,4 @@
 	md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1);
 
-	%love numbers: 
-	md.slr.ocean_area_scaling=1;
 	%}}}
 	%geometry:  {{{
@@ -301,5 +303,4 @@
 end
 % }}}
-
 %%Assemble Earth in 3D {{{
 
@@ -326,6 +327,4 @@
 sl.transfer('mask.ice_levelset');
 sl.transfer('mask.ocean_levelset');
-sl.transfer('mask.land_levelset');
-sl.transfer('mask.glacier_levelset');
 sl.transfer('geometry.bed');
 sl.transfer('mesh.lat');
@@ -355,5 +354,4 @@
 
 % }}}
-
 %Solve Sea-level equation on Earth only:  {{{
 md=sl.earth; %we don't do computations on ice sheets or land.
@@ -363,9 +361,10 @@
 
 %elastic loading from love numbers: 
-nlov=1001;
-md.slr.love_h = love_numbers('h'); md.slr.love_h(nlov+1:end)=[];
-md.slr.love_k = love_numbers('k'); md.slr.love_k(nlov+1:end)=[];
-md.slr.love_l = love_numbers('l'); md.slr.love_l(nlov+1:end)=[];
+nlov=101;
+md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[];
+md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[];
+md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[];
 md.slr.ocean_area_scaling = 0; 
+md.slr.loop_increment=200;
 
 %Miscellaneous
@@ -373,5 +372,5 @@
 
 %New stuff
-md.dsl.global_average_thermosteric_sea_level_change=[0;0];
+md.dsl.global_average_thermosteric_sea_level_change=[1.1+.38;0]; %steric + water storage AR5.
 
 %Solution parameters
@@ -379,4 +378,5 @@
 md.slr.abstol=1e-3;
 md.slr.geodetic=1;
+md.timestepping.time_step=1;
 
 % max number of iteration reverted back to 10 (i.e., the original default value)
@@ -385,4 +385,7 @@
 %eustatic run:
 md.slr.rigid=0; md.slr.elastic=0;md.slr.rotation=0;
+md.slr.requested_outputs= {'default',...
+	    'SealevelriseDeltathickness','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',...
+		    'SealevelNEsaRate', 'SealevelUEsaRate', 'SealevelNGiaRate', 'SealevelUGiaRate','SealevelEustaticMask','SealevelEustaticOceanMask'};
 md=solve(md,'Sealevelrise');
 Seustatic=md.results.SealevelriseSolution.Sealevel;
