Index: /issm/trunk-jpl/test/SandBox/test2004.m =================================================================== --- /issm/trunk-jpl/test/SandBox/test2004.m (revision 24762) +++ /issm/trunk-jpl/test/SandBox/test2004.m (revision 24763) @@ -3,221 +3,53 @@ %Data paths: {{{ modeldatapath='/Users/larour/ModelData'; -shppath='/Users/larour/issm-jpl/proj-group/qgis/Slr'; +shppath='/Users/larour/issm-jpl/proj-group/qgis/NightlyRun'; gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1.shp']; %}}} +skipmesh=1; + %create sealevel model to hold our information: -sl=sealevelmodel(); - +if ~skipmesh, + sl=sealevelmodel(); %Create basins using boundaries from shapefile: %{{{ -sl.addbasin(basin('continent','southamerica','name','southamerica','epsg',5530,'boundaries',{... %{{{ - boundary('shppath',shppath,'shpfilename','PanamaWest','epsg',4326),... - boundary('shppath',shppath,'shpfilename','SouthAmericaNorthWest','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','NorthAmericaSouthWest','epsg',4326),... - boundary('shppath',shppath,'shpfilename','SouthAmericaWest','epsg',4326),... - boundary('shppath',shppath,'shpfilename','MarieByrd1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Ellsworth2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Ellsworth1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Bellingshausen2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Bellingshausen1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Wilkins1','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Wilkins2Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Peninsula','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Palmer1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Palmer1','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Ronne2Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Ronne2b','epsg',3031),... - boundary('shppath',shppath,'shpfilename','SlessorSummit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Slessor1','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','BruntSummit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Brunt1','epsg',3031),... - boundary('shppath',shppath,'shpfilename','DronningMaud1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','SouthAmericaEast','epsg',4326),... - boundary('shppath',shppath,'shpfilename','MidAtlantic','epsg',4326),... - boundary('shppath',shppath,'shpfilename','SouthAmericaNorthEast','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','PanamaEast','epsg',4326),... - boundary('shppath',shppath,'shpfilename','Panama','epsg',4326)})); -%}}} -sl.addbasin(basin('continent','northamerica','name','northamerica','epsg',3628,'boundaries',{... %{{{ - boundary('shppath',shppath,'shpfilename','NorthAmericaSouthWest','epsg',4326),... - boundary('shppath',shppath,'shpfilename','SouthAmericaNorthWest','epsg',4326),... - boundary('shppath',shppath,'shpfilename','PanamaWest','epsg',4326),... - boundary('shppath',shppath,'shpfilename','Panama','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','PanamaEast','epsg',4326),... - boundary('shppath',shppath,'shpfilename','SouthAmericaNorthEast','epsg',4326),... - boundary('shppath',shppath,'shpfilename','MidAtlantic','epsg',4326),... - boundary('shppath',shppath,'shpfilename','Atlantic','epsg',4326),... - boundary('shppath',shppath,'shpfilename','South2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','South2','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','South1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','SouthWest1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','SouthWest1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Russell1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Russell1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Jakobshavn1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Jakobshavn1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Store1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Store1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','NorthWest1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','NorthWest1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Heilprin1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Petermann2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Petermann1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Petermann1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Ostenfeld1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Ostenfeld1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Academy1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Academy1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','ArcticCanada','epsg',4326),... - boundary('shppath',shppath,'shpfilename','NorthAmericaNorthWest','epsg',4326),... - boundary('shppath',shppath,'shpfilename','NorthAmericaWest','epsg',4326)})); -%}}} -sl.addbasin(basin('continent','australia','name','australia','epsg',4462,'boundaries',{... %{{{ - boundary('shppath',shppath,'shpfilename','Ross3Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Australia','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','SouthAsiaSummit','epsg',4326),... - boundary('shppath',shppath,'shpfilename','SouthAsia','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','IndianOceanSummit','epsg',4326),... - boundary('shppath',shppath,'shpfilename','IndianOcean','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','AmerySummit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','QueenMary2','epsg',3031),... - boundary('shppath',shppath,'shpfilename','QueenMary1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Wilkes3','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Wilkes2Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Adelie2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Adelie1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Victoria2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Victoria1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Ross4','epsg',3031,'orientation','reverse')})); -%}}} -sl.addbasin(basin('continent','eurasia','name','eurasia','epsg',3416,'boundaries',{... %{{{ - boundary('shppath',shppath,'shpfilename','South2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Atlantic','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','MidAtlantic','epsg',4326),... - boundary('shppath',shppath,'shpfilename','SouthAmericaEast','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','DronningMaud1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','DronningMaud1','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','QueenMaud2Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','QueenMaud2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Kemp2Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Kemp2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','KempSummit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Amery5','epsg',3031),... - boundary('shppath',shppath,'shpfilename','AmerySummit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','IndianOcean','epsg',4326),... - boundary('shppath',shppath,'shpfilename','IndianOceanSummit','epsg',4326),... - boundary('shppath',shppath,'shpfilename','SouthAsia','epsg',4326),... - boundary('shppath',shppath,'shpfilename','SouthAsiaSummit','epsg',4326),... - boundary('shppath',shppath,'shpfilename','ArcticRussia','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','NorthAmericaNorthWest','epsg',4326),... - boundary('shppath',shppath,'shpfilename','ArcticCanada','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Academy1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','NEGIS1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','NEGIS1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Bistrup2','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Bistrup1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Daugaard-Jensen5','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Daugaard-Jensen-Summit4','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Geikie2','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Geikie1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Kangerlussuaq3','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Kangerlussuaq2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Helheim3','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Helheim2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','KogeBugt3','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','KogeBugt2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','SouthEast3','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','SouthEast2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','South3','epsg',3413,'orientation','reverse')})); -%}}} -sl.addbasin(basin('continent','pacific','name','pacific','epsg',3031,'boundaries',{... %{{{ - boundary('shppath',shppath,'shpfilename','Ross3Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','MarieByrd2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','MarieByrd1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','SouthAmericaWest','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','NorthAmericaSouthWest','epsg',4326),... - boundary('shppath',shppath,'shpfilename','NorthAmericaWest','epsg',4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','NorthAmericaNorthWest','epsg',4326),... - boundary('shppath',shppath,'shpfilename','ArcticRussia','epsg',4326),... - boundary('shppath',shppath,'shpfilename','SouthAsiaSummit','epsg',4326),... - boundary('shppath',shppath,'shpfilename','Australia','epsg',4326)})); -%}}} -sl.addbasin(basin('continent','greenland','name','greenland','epsg',3413,'boundaries',{... %{{{ - boundary('shppath',shppath,'shpfilename','South2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','South2','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','South1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','SouthWest1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','SouthWest1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Russell1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Russell1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Jakobshavn1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Jakobshavn1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Store1','epsg',3413,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Store1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','NorthWest1','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','NorthWest1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Heilprin1','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','Petermann2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Petermann1','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','Petermann1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Ostenfeld1','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','Ostenfeld1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Academy1','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','Academy1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','NEGIS1','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','NEGIS1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Bistrup2','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','Bistrup1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Daugaard-Jensen5','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','Daugaard-Jensen-Summit4','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Geikie2','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','Geikie1Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Kangerlussuaq3','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','Kangerlussuaq2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','Helheim3','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','Helheim2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','KogeBugt3','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','KogeBugt2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','SouthEast3','epsg',3413,'orientation','reverse'),..., - boundary('shppath',shppath,'shpfilename','SouthEast2Summit','epsg',3413),... - boundary('shppath',shppath,'shpfilename','South3','epsg',3413,'orientation','reverse')})); -%}}} -sl.addbasin(basin('continent','antarctica','name','antarctica','epsg',3031,'boundaries',{... %{{{ - boundary('shppath',shppath,'shpfilename','DronningMaud1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','DronningMaud1','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','QueenMaud2Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','QueenMaud2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Kemp2Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Kemp2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','KempSummit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Amery5','epsg',3031),... - boundary('shppath',shppath,'shpfilename','AmerySummit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','QueenMary2','epsg',3031),... - boundary('shppath',shppath,'shpfilename','QueenMary1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Wilkes3','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Wilkes2Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Adelie2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Adelie1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Victoria2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Victoria1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Ross4','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Ross3Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','MarieByrd2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','MarieByrd1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Ellsworth2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Ellsworth1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Bellingshausen2','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','Bellingshausen1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Wilkins1','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Wilkins2Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Peninsula','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Palmer1Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Palmer1','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Ronne2Summit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Ronne2b','epsg',3031),... - boundary('shppath',shppath,'shpfilename','SlessorSummit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Slessor1','epsg',3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','BruntSummit','epsg',3031),... - boundary('shppath',shppath,'shpfilename','Brunt1','epsg',3031)})); +%HemisphereWest: {{{ + sl.addbasin(basin('continent','hemispherewest','name','hemispherewest','epsg',3587,'boundaries',{... %Peru projection 3587 + boundary('shppath',shppath,'shpfilename','HemisphereSplit','epsg',4326,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031),... + boundary('shppath',shppath,'shpfilename','RonneBrunt','epsg',3031,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031),... + boundary('shppath',shppath,'shpfilename','RonneFront','epsg',3031,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031),... + boundary('shppath',shppath,'shpfilename','WestAntarctica2','epsg',3031,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031)... + })); + %}}} + %HemisphereEast: {{{ + sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','epsg',4462,'boundaries',{... %Australian projection lat,long + boundary('shppath',shppath,'shpfilename','HemisphereSplit','epsg',4326),... + boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031),... + boundary('shppath',shppath,'shpfilename','EastAntarctica2','epsg',3031,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031)... + })); + %}}} + %Antarctica excluding Ronne: {{{ + sl.addbasin(basin('continent','antarctica','name','antarctica-grounded','epsg',3031,'boundaries',{... + boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031),... + boundary('shppath',shppath,'shpfilename','EastAntarctica2','epsg',3031),... + boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031),... + boundary('shppath',shppath,'shpfilename','WestAntarctica2','epsg',3031)... + boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031)... + boundary('shppath',shppath,'shpfilename','RonneIceShelf','epsg',3031)... + boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031)... + boundary('shppath',shppath,'shpfilename','RonneBrunt','epsg',3031)... + })); + %}}} + %Ronne: {{{ + sl.addbasin(basin('continent','antarctica','name','ronne','epsg',3031,'boundaries',{... + boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031),... + boundary('shppath',shppath,'shpfilename','RonneIceShelf','epsg',3031),... + boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031),... + boundary('shppath',shppath,'shpfilename','RonneFront','epsg',3031,'orientation','reverse')... + })); %}}} %}}} @@ -244,5 +76,4 @@ md=bamg(model,'domain',domain,'subdomains',coastline,'hmin',hmin,'hmax',hmax,defaultoptions{:}); plotmodel(md,'data','mesh');pause(1); - error; %miscellaneous: @@ -259,73 +90,37 @@ end %}}} -error; +end %Parameterize ice sheets : {{{ -for ind=sl.basinindx('continent',{'greenland','antarctica'}), +for ind=sl.basinindx('continent',{'antarctica'}), disp(sprintf('Parameterizing basin %s\n', sl.icecaps{ind}.miscellaneous.name)); md=sl.icecaps{ind}; bas=sl.basins{ind}; - - %basin/continent specific code: {{{ - if bas.iscontinentany('antarctica'), - openstreet=openstreetantarctica; - domainoutline=domainoutlinea; - end; - if bas.iscontinentany('greenland'), - openstreet=openstreetgreenland; - domainoutline=domainoutlineg; - end - %}}} %masks : %{{{ %new type of mask: md.mask=maskpsl; - %land and ocean from open street: - flags=1-(~ContourToNodes(md.mesh.x,md.mesh.y,openstreet,0)); - pos=find(flags==0); flags(pos)=-1; - md.mask.land_levelset=flags; - md.mask.ocean_levelset=-flags; + %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: - flags= -(~ContourToNodes(md.mesh.x,md.mesh.y,[shppath '/' domainoutline],1)); - mds=extract(md,[shppath '/' domainoutline]); - flags(mds.mesh.extractedvertices)=~mds.mesh.vertexonboundary; - pos=find(flags==0 & md.mesh.vertexonboundary); flags(pos)=1; - md.mask.ice_levelset=-flags; - - %on ice front, we are not on ocean; + 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: - mesh_glacier_levelset=md.private.bamg.landmask; %see meshing for why md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1); - pos=find(mesh_glacier_levelset); - md.mask.glacier_levelset(md.mesh.elements(pos,:))=1; - - pos=find(md.mask.glacier_levelset); - if ~isempty(pos), - mds=extract(md,md.mask.glacier_levelset); - md.mask.ice_levelset(mds.mesh.extractedvertices)=-(~mds.mesh.vertexonboundary); - end %}}} %initial grounding line position: {{{ - if bas.iscontinentany('antarctica'), % {{{ - - %figure out initial grounding line position from the bedmap2 dataset: - land_type=interpBedmap2(md.mesh.x,md.mesh.y,'icemask_grounded_and_shelves'); - - pos=find(isnan(land_type)); - pos2=find(~isnan(land_type)); - for i=1:length(pos), - in=find_point(md.mesh.x(pos2),md.mesh.y(pos2),md.mesh.x(pos(i)),md.mesh.y(pos(i))); - land_type(pos(i))=land_type(pos2(in)); - end - - gridoniceshelf=zeros(md.mesh.numberofvertices,1); - gridoniceshelf(land_type>0.9)=land_type(land_type>0.9); - gridoniceshelf(gridoniceshelf>0)=-1; - gridoniceshelf(gridoniceshelf~=-1)=1; - - md.mask.groundedice_levelset=gridoniceshelf; + if bas.isnameany('antarctica-grounded'), + + md.mask.groundedice_levelset=ones(md.mesh.numberofvertices,1); + end + if bas.isnameany('ronne'), + + md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1); %correction to land and ocean levelset: ice shelf is not on land! @@ -333,15 +128,4 @@ md.mask.ocean_levelset(pos)=1; md.mask.land_levelset(pos)=-1; - %}}} - elseif bas.iscontinentany('greenland'), % {{{ - - mask=interpBedmachine(md.mesh.x,md.mesh.y,'mask','greenland'); - - gridoniceshelf=ones(md.mesh.numberofvertices,1); - pos=find(md.mask.ice_levelset>0); groundedice_levelset(pos)=-1; - pos=find(mask<=0); gridoniceshelf(pos)=0; - - md.mask.groundedice_levelset=gridoniceshelf; - %}}} end %}}} @@ -349,570 +133,46 @@ [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326'); %}}} - %flowequation: % {{{ - md=setflowequation(md,'SSA','all'); - %}}} %geometry: {{{ - if bas.iscontinentany('antarctica'),% {{{ + if bas.iscontinentany('antarctica'), di=md.materials.rho_ice/md.materials.rho_water; - disp(' reading dem'); - surface=interpBedmachine(md.mesh.x,md.mesh.y,'surface','antarctica'); - surface=griddata(md.mesh.x(~isnan(surface)),md.mesh.y(~isnan(surface)),surface(~isnan(surface)),md.mesh.x,md.mesh.y,'nearest'); - surfaceold=surface; - - disp(' reading firn layer'); - load(firnpath) - firn_layer=InterpFromGridToMesh(x1,y1,firn,md.mesh.x,md.mesh.y,0); - disp(' reading bedrock'); - base=interpBedmachine(md.mesh.x,md.mesh.y,'bed','antarctica'); - bedmap=interpBedmap2(md.mesh.x,md.mesh.y,'bed'); - base(base<-8000 | isnan(base))=bedmap(base<-8000 | isnan(base)); - - %thickness: - thickness=surface-base; - - %hydrostatic equilibrium on ice shelves: - pos=find(md.mask.groundedice_levelset<0); thickness(pos)=surface(pos)/(1-di)-firn_layer(pos); - surface(pos)=(1-di)*thickness(pos); - base(pos)=-di*thickness(pos); - - %firn layer over grounded ice: - pos=find(md.mask.groundedice_levelset>=0); thickness(pos)=thickness(pos)-firn_layer(pos); - - %check again: - pos=find(thickness<1); thickness(pos)=1; - - %reset: - surface=base+thickness; - - %make sure we are at hydrostatic at the grounding line: - bed=base;bed(md.mask.groundedice_levelset<0)=base(md.mask.groundedice_levelset<0)-2000; - - %set: - md.geometry.surface=surface; - md.geometry.bed=bed; - md.geometry.base=base; - md.geometry.thickness=thickness; - %}}} - elseif bas.iscontinentany('greenland'),% {{{ - mask=interpBedmachine(md.mesh.x,md.mesh.y,'mask','greenland'); - surface=interpBedmachine(md.mesh.x,md.mesh.y,'surface','greenland'); - thickness=interpBedmachine(md.mesh.x,md.mesh.y,'thickness','greenland'); - base=surface-thickness; - di=md.materials.rho_ice/md.materials.rho_water; - - pos=find(isnan(base) | isnan(surface)); - thickness(pos)=1; - surface(pos)=(1-di)*thickness(pos); - base(pos)=-di*thickness(pos); - - %places where thickness is 0: - pos=find(thickness==0); - thickness(pos)=.1; surface(pos)=base(pos)+thickness(pos); - - %adjust bed and base: - bed=base; - pos=find(md.mask.groundedice_levelset>=0); bed(pos)=base(pos); - pos=find(md.mask.groundedice_levelset<0); bed(pos)=base(pos)-100; - - md.geometry.surface=surface; - md.geometry.bed=bed; - md.geometry.base=bed; - md.geometry.thickness=thickness; + md.geometry.bed=interpBedmap2(md.mesh.x,md.mesh.y,'bed'); + end % }}} + %Slr: {{{ + if bas.iscontinentany('antarctica'), + + 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; + + 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; + + 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); + md.slr.Ugia=zeros(md.mesh.numberofvertices,1); + + md.dsl.global_average_thermosteric_sea_level_change=[0;0]; + md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); + md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); + md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1); end %}}} + % material properties: {{{ + md.materials=materials('hydro'); + %}}} + %diverse: {{{ + md.miscellaneous.name=bas.name; % }}} - %velocities: %{{{ - if bas.iscontinentany('antarctica'),% {{{ - - velnc =[jplsvn '/database/antarctica/velocity_ref_v28Apr2011bamber_900m_median.nc']; - - %read and process netcdf - xmin = double(ncreadatt(velnc,'/','xmin')); - ymax = double(ncreadatt(velnc,'/','ymax')); - nx = double(ncreadatt(velnc,'/','nx')); - ny = double(ncreadatt(velnc,'/','ny')); - spacing = double(ncreadatt(velnc,'/','spacing')); - vx = double(ncread(velnc,'vx')); - vy = double(ncread(velnc,'vy')); - x=xmin+(0:1:nx)'*spacing; - y=(ymax-ny*spacing)+(0:1:ny)'*spacing; - vx=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0); - vy=InterpFromGridToMesh(x,y,flipud(vy'),md.mesh.x,md.mesh.y,0); - vel=sqrt(vx.^2+vy.^2); - - %}}} - elseif bas.iscontinentany('greenland'),% {{{ - - velocities=refinementmetric(md.mesh,'greenland','greenland','nsidcvxvy'); - vx=velocities(:,1); vy=velocities(:,2); vel=sqrt(vx.^2+vy.^2); - - velocities=refinementmetric(md.mesh,'greenland','greenland','joughinvxvy'); - vxfar=velocities(:,1); vyfar=velocities(:,2); velfar=sqrt(vxfar.^2+vyfar.^2); - - %use velfar inland (from mosaic), and 2008 velocities nearcoastline: - pos=find(vel==0 | isnan(vel)); vel(pos)=velfar(pos); vx(pos)=vxfar(pos);vy(pos)=vyfar(pos); - - %set water velocity to 0: - pos=find(isnan(vel)); vel(pos)=0; vx(pos)=0; vy(pos)=0; - - end % }}} - - md.inversion.vx_obs=vx; - md.inversion.vy_obs=vy; - md.inversion.vel_obs=vel; - md.initialization.vx=md.inversion.vx_obs; - md.initialization.vy=md.inversion.vy_obs; - md.initialization.vz=zeros(md.mesh.numberofvertices,1); - md.initialization.vel=md.inversion.vel_obs; - - % }}} - %friction: {{{ - [sx,sy,s]=slope(md); sslope=averaging(md,s,0); - - pos=find(md.inversion.vel_obs==0); - vel=max(md.inversion.vel_obs,0.1); vel(pos)=1; - - md.friction.coefficient=sqrt(md.materials.rho_ice*md.geometry.thickness.*(sslope)./((md.materials.rho_ice*md.geometry.thickness+md.materials.rho_water*md.geometry.base).*vel/md.constants.yts)); - md.friction.coefficient=min(md.friction.coefficient,200); - - min_drag_coeff=35; - min_drag_coeff_outlets=5; - threshhold_drag_coeff=70; - background_drag_coeff=200; - - maxvel=30; - pos=find(md.friction.coefficient 1990 as average, and perturbation from 2005 -> 2010 - - ncdata = sprintf('%s/%s-%i%s',smbpath,'smb.KNMI',1979,'.ANT3K27.DRIFTALB.nc'); - lat = ncread(ncdata,'lat'); lon = ncread(ncdata,'lon'); - lat = double(lat); lon = double(lon); - lat = lat(:); lon = lon(:); - [x,y]=ll2xy(lat,lon,-1); - index=BamgTriangulate(x,y); - di=md.materials.rho_ice/md.materials.rho_water; - - smbs=zeros(length(x),numel(Years)); - for i=1:length(Years), - year=Years(i); - ncdata = sprintf('%s/%s-%i%s',smbpath,'smb.KNMI',year,'.ANT3K27.DRIFTALB.nc'); - smb=ncread(ncdata,'smb'); smb=sum(smb(:,:,:),3)/1000/di; smb=smb(:); - smbs(:,i)=smb; - end - smb=InterpFromMeshToMesh2d(index,x,y,smbs,md.mesh.x,md.mesh.y); - - %mean from 1979 to 1990: - pos=find(Years>=1979 & Years<=1990); - smbmean7990=mean(smb(:,pos),2); - smbmean7990=griddata(md.mesh.x(~isnan(smbmean7990)),md.mesh.y(~isnan(smbmean7990)),smbmean7990(~isnan(smbmean7990)),md.mesh.x,md.mesh.y,'nearest'); - - %mean from 2005-2010: - pos=find(Years>=2005 & Years<=2010); - smbmean0510=mean(smb(:,pos),2); - smbmean0510=griddata(md.mesh.x(~isnan(smbmean0510)),md.mesh.y(~isnan(smbmean0510)),smbmean0510(~isnan(smbmean0510)),md.mesh.x,md.mesh.y,'nearest'); - - md.smb.mass_balance = smbmean0510;; - - pos=find(md.mask.glacier_levelset); - md.smb.mass_balance(pos)=smbmean0510(pos)-smbmean7990(pos); - - %}}} - %temperature: {{{ - temperaturepath=['/Users/larour/ModelData/RACMO2Antarctica/annt33.mat']; - load(temperaturepath); - index=delaunay(x1,y1); - md.initialization.temperature=InterpFromMeshToMesh2d(index,x1(:),y1(:),annt33(:),md.mesh.x,md.mesh.y,'default',273.15); - md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base); - % }}} - %}}} - elseif bas.iscontinentany('greenland'),% {{{ - %% Annual mean for the period given by Years (smb0) - Years = 1960:2014; %we are going to need 1960 -> 1990 as average, and perturbation from 2010 -> 2014. - - smb2 = nan(md.mesh.numberofvertices,numel(Years)); - st2 = nan(md.mesh.numberofvertices,numel(Years)); - - %triangulate first: - ncdata = [MARname '/MARv3.5-15km-monthly-ERA-40-1971.nc']; - lat = ncread(ncdata,'LAT'); lon = ncread(ncdata,'LON'); - lat = double(lat); lon = double(lon); - lat = lat(:); lon = lon(:); - [x,y]=ll2xy(lat,lon,+1,45,70); - index=BamgTriangulate(x,y); - - smb1s=zeros(length(x),numel(Years)); - st1s=zeros(length(x),numel(Years)); - - for ii = 1:numel(Years) - if Years(ii)<=1978 - ncdata = [MARname '/MARv3.5-15km-monthly-ERA-40-' num2str(Years(ii)) '.nc']; - elseif Years(ii)>1978 - ncdata = [MARname '/MARv3.5-15km-monthly-ERA-Interim-' num2str(Years(ii)) '.nc']; - end - - st1 = ncread(ncdata,'STcorr'); % Surface temperature - st1(st1<-0.9000e30)=nan; - st1 = double(mean(st1,3)); % Annual mean surface temp - - st1s(:,ii)=st1(:); - - smb1 = ncread(ncdata,'SMBcorr'); % in mmWe/month - smb1(smb1<-0.900e30)=nan; - smb1 = double(sum(smb1,3)/1000); %mWe/year - - smb1s(:,ii)=smb1(:); - - end - - st2s=InterpFromMeshToMesh2d(index,x,y,st1s,md.mesh.x,md.mesh.y); - smb2s=InterpFromMeshToMesh2d(index,x,y,smb1s,md.mesh.x,md.mesh.y); - - smb19601990=mean(smb2s(:,1:31),2); - smb20102014=mean(smb2s(:,end-4:end),2); - - st20102014 = mean(st2s(:,end-4:end),2); - - md.initialization.temperature=273.25+st20102014; - md.initialization.pressure=md.materials.rho_ice*md.constants.g*(md.geometry.surface-md.geometry.base); - md.smb.mass_balance = smb20102014; - - pos=find(md.mask.glacier_levelset); - md.smb.mass_balance(pos)=smb20102014(pos)-smb19601990(pos); - - end % }}} - %}}} - %Mechanical boundary conditions: {{{ - %intialize spc boundary conditions: - md.stressbalance.spcvx=NaN*ones(md.mesh.numberofvertices,1); - md.stressbalance.spcvy=NaN*ones(md.mesh.numberofvertices,1); - md.stressbalance.spcvz=NaN*ones(md.mesh.numberofvertices,1); - md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6); - - %boundary conditions according to masks: - - %water (freeze velocity to 0) - pos=find(md.mask.ocean_levelset>=0); - md.stressbalance.spcvx(pos)=0; - md.stressbalance.spcvy(pos)=0; - md.stressbalance.spcvz(pos)=0; - - %glaciers (freeeze velocity to observed) - pos=find(md.mask.glacier_levelset); - md.stressbalance.spcvx(pos)=0; %flux divergence will be equated to smb perturbation. - md.stressbalance.spcvy(pos)=0; - - md.stressbalance.spcvz(pos)=0; - - %land with no ice (freeeze velocity to 0) - pos=find(md.mask.land_levelset>=0 & md.mask.ice_levelset>0 & ~md.mask.glacier_levelset); - md.stressbalance.spcvx(pos)=0; - md.stressbalance.spcvy(pos)=0; - md.stressbalance.spcvz(pos)=0; - - %if strictly ice (and not glacier), then no boundary condition: - pos=find(md.mask.ice_levelset<=0 & ~md.mask.glacier_levelset); - md.stressbalance.spcvx(pos)=NaN; - md.stressbalance.spcvy(pos)=NaN; - md.stressbalance.spcvz(pos)=NaN; - - %constrain boundaries of basin to observed velocities: - pos=find(md.mesh.vertexonboundary & md.mask.ice_levelset<=0); - md.stressbalance.spcvx(pos)=md.inversion.vx_obs(pos); - md.stressbalance.spcvy(pos)=md.inversion.vy_obs(pos); - md.stressbalance.spcvz(pos)=0; - - %fixing quirks here and there: - if strcmpi(bas.name,'Ross'), - flags = ContourToNodes(mdband.mesh.x,mdband.mesh.y,[antarcticashppath '/Murdo.exp'],1); - pos=find(flags); - - md.stressbalance.spcvx(pos)=md.inversion.vx_obs(pos); - md.stressbalance.spcvy(pos)=md.inversion.vy_obs(pos); - md.stressbalance.spcvz(pos)=0; - end - - %prognostic - md.masstransport.spcthickness=NaN*ones(md.mesh.numberofvertices,1); - % }}} - %Thermal boundary conditions: {{{ - if bas.iscontinentany('antarctica'),% {{{ - searisepath =[modeldatapath '/SeaRISE/Antarctica5km_shelves_v1.0']; - heatfluxpath =[searisepath '/bheatflx_fox.mat']; - - if strcmpi(bas.name,'ellsworth'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=79.3397; - md.basalforcings.deepwater_elevation=-1000; - md.basalforcings.upperwater_elevation=-141.59; - %}}} - elseif strcmpi(bas.name,'amery'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=8; - md.basalforcings.deepwater_elevation=-1500; - md.basalforcings.upperwater_elevation=-400; % }}} - elseif strcmpi(bas.name,'ross'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=5; - md.basalforcings.deepwater_elevation=-1100; - md.basalforcings.upperwater_elevation=-200; % }}} - elseif strcmpi(bas.name,'ronne'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=4; - md.basalforcings.deepwater_elevation=-1800; - md.basalforcings.upperwater_elevation=-400; % }}} - elseif strcmpi(bas.name,'brunt'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=.3; - md.basalforcings.deepwater_elevation=-900; - md.basalforcings.upperwater_elevation=-100; - % }}} - elseif strcmpi(bas.name,'kemp'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=.2; - md.basalforcings.deepwater_elevation=-800; - md.basalforcings.upperwater_elevation=-100; % }}} - elseif strcmpi(bas.name,'queenMaud'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=.5; - md.basalforcings.deepwater_elevation=-650; - md.basalforcings.upperwater_elevation=-180; % }}} - elseif strcmpi(bas.name,'dronningMaud'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=3; - md.basalforcings.deepwater_elevation=-800; - md.basalforcings.upperwater_elevation=-300; % }}} - elseif strcmpi(bas.name,'slessor'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=1.4; - md.basalforcings.deepwater_elevation=-1200; - md.basalforcings.upperwater_elevation=-500; % }}} - elseif strcmpi(bas.name,'victoria'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=10; - md.basalforcings.deepwater_elevation=-1000; - md.basalforcings.upperwater_elevation=-200; % }}} - elseif strcmpi(bas.name,'adelie'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=4; - md.basalforcings.deepwater_elevation=-1000; - md.basalforcings.upperwater_elevation=-200; % }}} - elseif strcmpi(bas.name,'wilkes'), % {{{ - md.basalforcings=linerbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=14; - md.basalforcings.deepwater_elevation=-1500; - md.basalforcings.upperwater_elevation=-50; % }}} - elseif strcmpi(bas.name,'wilkins'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=6; - md.basalforcings.deepwater_elevation=-600; - md.basalforcings.upperwater_elevation=-50; % }}} - elseif strcmpi(bas.name,'queenmary'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=10; - md.basalforcings.deepwater_elevation=-1200; - md.basalforcings.upperwater_elevation=-200; % }}} - elseif strcmpi(bas.name,'mariebyrd'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=8; - md.basalforcings.deepwater_elevation=-700; - md.basalforcings.upperwater_elevation=-200; % }}} - elseif strcmpi(bas.name,'bellingshausen'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=7; - md.basalforcings.deepwater_elevation=-550; - md.basalforcings.upperwater_elevation=-150; % }}} - elseif strcmpi(bas.name,'palmer'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=.1; - md.basalforcings.deepwater_elevation=-1000; - md.basalforcings.upperwater_elevation=-50; % }}} - elseif strcmpi(bas.name,'peninsula'), % {{{ - md.basalforcings=linearbasalforcings(md); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.deepwater_melting_rate=2; - md.basalforcings.deepwater_elevation=-550; - md.basalforcings.upperwater_elevation=-50; % }}} - else - %melt rates: {{{ - md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - load('./Data/monthly_mean_melt_adj09.mat'); fw2=fw_flux; - load('./Data/monthly_mean_melt_adj04.mat'); fw1=fw_flux; - fw_flux(:,:,end+1:end+size(fw2,3))=fw2; - load('./Data/grid_info_cp09.mat') - [xi,yi]= ll2xy(lat,lon,-1,0,71); index=delaunay(xi,yi); - mrate_mesh=InterpFromMeshToMesh2d(index,xi(:),yi(:),-1*reshape(squeeze(mean(fw_flux,3)),prod(size(xi)),1),md.mesh.x,md.mesh.y,'default',0); - pos=find(md.mask.groundedice_levelset<0 & mrate_mesh~=0); - md.basalforcings.floatingice_melting_rate(pos)=mrate_mesh(pos); - % }}} -end -md.thermal.spctemperature=[md.initialization.temperature;1]; %impose observed temperature on surface -load(heatfluxpath) -md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,bheatflx_fox,md.mesh.x,md.mesh.y,0); -% }}} -elseif bas.iscontinentany('greenland'),% {{{ - %Initialize melt rates: {{{ - md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); - md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); - - %Specificy melt rates and distances: {{{ - %Seroussi et al., 2011 - melt rates at grounding line of 79N are about 50 m/yr - %Rignot and Steffen, 2008 - Petermann high rates up to 20 km out - if strcmpi(bas.name,'negis'), - iceshelves={'gl_zac_1996','gl_79n_1996','gl_StorBistrup96'}; - distances=[10,10,10]; - meltrates=[7,2,4]; - meltratesgl=[25,22,4]; %35,35 - elseif strcmpi(bas.name,'ostenfeld'), - distances=[10,10,10,10,10]; - iceshelves={'gl_Ryder96','gl_Ostenfeld96','gl_Harder96','gl_Steensby96','gl_Brikkerne96'}; - meltrates=[7,11,4,5,4]; - meltratesgl=[25,26,4,5,4]; - elseif strcmpi(bas.name,'petermann'), - distances=[20,10]; - iceshelves={'gl_Peter96','gl_Humb96'}; - meltrates=[6,4]; %5 - meltratesgl=[22,4]; %26 - elseif strcmpi(bas.name,'academy'), - distances=[10]; - iceshelves={'gl_Hagen96'}; - meltrates=[4]; - meltratesgl=[4]; - elseif strcmpi(bas.name,'jakobshavn'), - distances=[10]; - iceshelves={'gl_Jak'}; - meltrates=[4]; - meltratesgl=[109]; %Prescott et al, 2003 - else - iceshelves={}; - end - %}}} - - for i=1:length(iceshelves), - iceshelf=iceshelves{i}; - distance=distances(i); - meltrate=meltrates(i); - meltrategl=meltratesgl(i); - - in= ContourToNodes(md.mesh.x,md.mesh.y,[shppath '/' iceshelf '_closed_polygon.shp'],1); - segs=MeshProfileIntersection(md.mesh.elements,md.mesh.x,md.mesh.y,[shppath '/' iceshelf '_closed_polygon.shp']); - pos=find(~isnan(segs(:,1))); - segs=[segs(pos,1:2); segs(pos,3:4)]; - - %figure out nodes within ~10 km of grounding line: - v=find(in); - glpt=zeros(md.mesh.numberofvertices,1); - a=zeros(md.mesh.numberofvertices,1); - for i=v' - x=md.mesh.x(i); y=md.mesh.y(i); - dist=sqrt((segs(:,1)-x).^2+(segs(:,2)-y).^2); - if min(dist)3, - md.basalforcings.floatingice_melting_rate(pos)=griddata(md.mesh.x(pos2),md.mesh.y(pos2),md.basalforcings.floatingice_melting_rate(pos2),md.mesh.x(pos),md.mesh.y(pos),'nearest'); - end - end - end - %}}} - %Thermal spcs and heatflux: %{{{ - md.thermal.spctemperature=md.initialization.temperature; %impose observed temperature on surface - - load(seariseg_heatflux) - [xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,39,71); - md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,bheatflx,xi,yi,0); - %}}} -end % }}} -%}}} -%Slr: {{{ -if bas.iscontinentany('antarctica'), - - 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; - - 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; -else - delH=textread([modeldatapath '/AdhikariSciAdvTrends/GIS_delH_trend.txt']); - longGIS=delH(:,1); - latGIS=delH(:,2); - delHGIS=delH(:,3); - - index=delaunay(latGIS,longGIS); - - lat=md.mesh.lat; - long=md.mesh.long+360; - pos=find(long>360);long(pos)=long(pos)-360; - - delHGIS=InterpFromMesh2d(index,longGIS,latGIS,delHGIS,long,lat); - md.slr.deltathickness=delHGIS(md.mesh.elements)*[1;1;1]/3/100; -end -md.slr.sealevel=zeros(md.mesh.numberofvertices,1); -md.slr.steric_rate=zeros(md.mesh.numberofvertices,1); -md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1); -md.slr.Ngia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ngia,md.mesh.long,md.mesh.lat); -md.slr.Ugia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ugia,md.mesh.long,md.mesh.lat); -%}}} -% material properties: {{{ -disp(' creating flow law paramter'); -md.materials.rheology_B=paterson(md.initialization.temperature); -pos=find(md.materials.rheology_B<=0); -md.materials.rheology_B(pos)=-md.materials.rheology_B(pos); -md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); -%}}} -%diverse: {{{ -md.miscellaneous.name=bas.name; -% }}} sl.icecaps{ind}=md; @@ -921,11 +181,5 @@ % ParameterizeContinents {{{ -%glacier load: {{{ -delH=textread([modeldatapath '/AdhikariSciAdvTrends/GLA_delH_trend_15regions.txt']); -longglaciers=delH(:,1); latglaciers=delH(:,2); delHglaciers=sum(delH(:,3:end),2); -indexglaciers=delaunay(latglaciers,longglaciers); -%}}} - -for ind=sl.basinindx('continent',{'southamerica','northamerica','australia','eurasia','pacific'}), +for ind=sl.basinindx('continent',{'hemisphereeast','hemispherewest'}), disp(sprintf('Masks for basin %s\n', sl.icecaps{ind}.miscellaneous.name)); md=sl.icecaps{ind}; bas=sl.basins{ind}; @@ -934,10 +188,4 @@ [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326'); - %interpolate glacier loads onto mesh: %{{{ - lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360; - mdelHglaciers=InterpFromMesh2d(indexglaciers,longglaciers,latglaciers,delHglaciers,long,lat); - mdelHglacierse=mdelHglaciers(md.mesh.elements)*[1;1;1]/3; - pos=find(~md.private.bamg.landmask); mdelHglacierse(pos)=0; - % }}} %mask: %{{{ %new type of mask: @@ -990,71 +238,36 @@ md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); - %%cross check that whereever we have an ice load, the mask is <0 on each vertex: if it's not th ecase, zero the ice load: - pos=find(mdelHglacierse); - vertices=md.mesh.elements(pos,:); %potentially having an ice load. - pos=find(md.mask.land_levelset(vertices)>=0); - icevertices=vertices(pos); - md.mask.ice_levelset(icevertices)=-md.mask.land_levelset(icevertices); - - %then take care of glaciers: don't! this is determined by the land mask. - %pos=find(md.slr.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; - %pos=find(mdelHglacierse); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; - - %grounded ice: md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1); - pos=find(mdelHglaciers); md.mask.groundedice_levelset(pos)=1; - - - %now, glaciers: - md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1); - pos=find(mdelHglacierse); md.mask.glacier_levelset(md.mesh.elements(pos,:))=1; + % }}} %}}} %slr loading/calibration: {{{ - %load glaciers and ice caps from GRACE: - md.slr.deltathickness=zeros(md.mesh.numberofelements,1); - pos=find(mdelHglacierse); md.slr.deltathickness(pos)=mdelHglacierse(pos)/100; - - %initialize: + 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); + md.slr.Ugia=zeros(md.mesh.numberofvertices,1); + + md.dsl.global_average_thermosteric_sea_level_change=[0;0]; + %md.slr.steric_rate=(1.1+.38)*ones(md.mesh.numberofvertices,1); %steric + water storage. + md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); + md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); + md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1); %love numbers: - md.slr.steric_rate=(1.1+.38)*ones(md.mesh.numberofvertices,1); %steric + water storage. md.slr.ocean_area_scaling=1; - md.slr.Ngia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ngia,md.mesh.long,md.mesh.lat); - md.slr.Ugia=InterpFromMesh2d(indexGIA,longGIA,latGIA,Ugia,md.mesh.long,md.mesh.lat); %}}} %geometry: {{{ di=md.materials.rho_ice/md.materials.rho_water; - - md.slr.spcthickness=ones(md.mesh.numberofvertices+1,2); - %time tags: make sure they are far in the past and future. - t1=-10000; t2=+10000; md.slr.spcthickness(end,1)=-10000; - md.slr.spcthickness(end,2)=+10000; - - %impart the glacier load: - pos=find(mdelHglaciers); - md.slr.spcthickness(pos,2)=md.slr.spcthickness(pos,1)+(t2-t1)*mdelHglaciers(pos)/100*di; - %mdelHglaciers is cm/yr water equivalent - - md.geometry.thickness=ones(md.mesh.numberofvertices,1); - md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1); - md.geometry.base=md.geometry.surface-md.geometry.thickness; - md.geometry.bed=md.geometry.base; + md.geometry.bed=-ones(md.mesh.numberofvertices,1); % }}} %materials: {{{ - md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1); - md.materials.rheology_B=paterson(md.initialization.temperature); - md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); - % }}} - %Solution: {{{ - md.transient.issmb=0; md.transient.isstressbalance=0; - md.transient.ismasstransport=0; md.transient.isthermal=0; - md.transient.isslr=0; md.slr.loop_increment=300; + md.materials=materials('hydro'); % }}} sl.icecaps{ind}=md; end % }}} + %Assemble Earth in 3D {{{ @@ -1066,4 +279,6 @@ %create earth model by concatenating all the icecaps in 3d: sl.caticecaps('tolerance',tolerance,'loneedgesdetect',loneedgesdetect); + +error; %figure out how each icecap's mesh connects to the larger earth mesh: @@ -1110,4 +325,5 @@ % }}} +error; %Solve Sea-level equation on Earth only: {{{ md=sl.earth; %we don't do computations on ice sheets or land.