Index: /issm/trunk-jpl/test/SandBox/test2004.m =================================================================== --- /issm/trunk-jpl/test/SandBox/test2004.m (revision 24904) +++ /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;