source: issm/oecreview/Archive/24684-25833/ISSM-24808-24809.diff@ 25834

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};
-
 %}}}