Index: ../trunk-jpl/test/SandBox/test2004.m
===================================================================
--- ../trunk-jpl/test/SandBox/test2004.m	(revision 24808)
+++ ../trunk-jpl/test/SandBox/test2004.m	(revision 24809)
@@ -1,9 +1,11 @@
 %Test Name: Earth_Antarctica_GIA
+		
+testagainst2002=1;
 
 %Data paths: {{{
 modeldatapath='/Users/larour/ModelData';
 shppath='/Users/larour/issm-jpl/proj-group/qgis/NightlyRun';
-gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1.shp'];
+gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1-NightlyRun.shp'];
 %}}}
 
 %create sealevel model to hold our information: 
@@ -145,22 +147,30 @@
 	end % }}}
 	%Slr: {{{
 	if bas.iscontinentany('antarctica'),
+		if testagainst2002,
+			md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
+			%antarctica
+			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;
+		else
+			delH=textread([modeldatapath '/AdhikariSciAdvTrends/AIS_delH_trend.txt']);
+			longAIS=delH(:,1);
+			latAIS=delH(:,2);
+			delHAIS=delH(:,3);
+			index=delaunay(latAIS,longAIS);
 
-		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;
 
-		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;
 
-		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
 
-		md.slr.deltathickness=delHAIS(md.mesh.elements)*[1;1;1]/3/100;
-
 		md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
 		md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
 		md.slr.Ngia=zeros(md.mesh.numberofvertices,1);
@@ -252,8 +262,21 @@
 	%}}}
 	%slr loading/calibration:  {{{
 
+	if testagainst2002,
+		md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
+		%greenland
+		late=sum(md.mesh.lat(md.mesh.elements),2)/3;
+		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;
+
+		%correct mask: 
+		md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
+	else
+		md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
+	end
+
 	md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
-	md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
 	md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
 	md.slr.Ngia=zeros(md.mesh.numberofvertices,1);
 	md.slr.Ugia=zeros(md.mesh.numberofvertices,1);
@@ -357,28 +380,27 @@
 md.slr.abstol=1e-3;
 md.slr.geodetic=1;
 
-%eustatic: 
-md.slr.rigid=0; md.slr.elastic=0; md.slr.rotation=0; 
-md.cluster=generic('name',oshostname(),'np',10);
+% max number of iteration reverted back to 10 (i.e., the original default value)
+md.slr.maxiter=10;
+
+%eustatic run:
+md.slr.rigid=0; md.slr.elastic=0;md.slr.rotation=0;
 md=solve(md,'Sealevelrise');
+Seustatic=md.results.SealevelriseSolution.Sealevel;
 
-%eustatic + rigid + elastic run: 
-md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=0; 
-md.cluster=generic('name',oshostname(),'np',16);
-md.verbose=verbose('111111111');
+%eustatic + rigid run:
+md.slr.rigid=1; md.slr.elastic=0;md.slr.rotation=0;
 md=solve(md,'Sealevelrise');
-SnoRotation=md.results.SealevelriseSolution.Sealevel;
+Srigid=md.results.SealevelriseSolution.Sealevel;
 
-%eustatic + rigid + elastic + rotation run: 
+%eustatic + rigid + elastic run:
+md.slr.rigid=1; md.slr.elastic=1;md.slr.rotation=0;
+md=solve(md,'Sealevelrise');
+Selastic=md.results.SealevelriseSolution.Sealevel;
+
+%eustatic + rigid + elastic + rotation run:
 md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=1;
-md.cluster=generic('name',oshostname(),'np',3);
-%md.verbose=verbose('111111111');
 md=solve(md,'Sealevelrise');
-SRotation=md.results.SealevelriseSolution.Sealevel;
+Srotation=md.results.SealevelriseSolution.Sealevel;
 
-%Fields and tolerances to track changes
-field_names     ={'noRotation','Rotation'};
-field_tolerances={1e-13,1e-13};
-field_values={SnoRotation,SRotation};
-
 %}}}