Changeset 24905

Timestamp:

05/27/20 09:54:32 (5 years ago)

Author:

Eric.Larour

Message:

CHG: stable, trying to commit into nightly runs.

File:

• issm/trunk-jpl/test/SandBox/test2004.m (modified) (19 diffs)

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

@@ -1 +1 @@
 %Test Name: Earth_Antarctica_GIA
                 
-testagainst2002=1;
+testagainst2002=0;
 
 %Data paths: {{{
@@ -29 +29 @@
         }));
         %}}}
+        %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 +51 @@
         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 +73 @@
 %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 +116 @@
         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 +143 @@
                         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 +190 @@
 
         %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 +209 @@
         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 +222 @@
                 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 +238 @@
         %slr loading/calibration:  {{{
 
-        if testagainst2002,
+        if testagainst2002, 
+                % {{{
                 md.slr.deltathickness=zeros(md.mesh.numberofelements,1);
                 %greenland
@@ -269 +245 @@
                 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 +292 @@
         md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1);
 
-        %love numbers: 
-        md.slr.ocean_area_scaling=1;
         %}}}
         %geometry:  {{{
@@ -301 +303 @@
 end
 % }}}
-
 %%Assemble Earth in 3D {{{
 
@@ -326 +327 @@
 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 +354 @@
 
 % }}}
-
 %Solve Sea-level equation on Earth only:  {{{
 md=sl.earth; %we don't do computations on ice sheets or land.
@@ -363 +361 @@
 
 %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 +372 @@
 
 %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 +378 @@
 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 +385 @@
 %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;