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

../trunk-jpl/test/SandBox/test2004.m

@@ -1 +1 @@
 %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 +147 @@
         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 +262 @@
         %}}}
         %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 +380 @@
 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};
-
 %}}}