Index: /issm/trunk-jpl/test/NightlyRun/test2004.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 24717) +++ /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 24717) @@ -0,0 +1,1154 @@ +%Test Name: Earth_Antarctica_GIA + +%Data paths: {{{ +modeldatapath='/Users/larour/ModelData'; +shppath='/Users/larour/issm-jpl/proj-group/qgis/Slr'; +gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1.shp']; +%}}} + +%create sealevel model to hold our information: +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)})); + %}}} +%}}} +%Go through basins and mesh: %{{{ +%meshing parameters: {{{ +hmin=500; hmax=2000; hmin=hmin*1000; hmax=hmax*1000; +tolerance=100; %tolerance of 100m on Earth position when mergin 3d meshes +threshold=1; +defaultoptions={'KeepVertices',0,'MaxCornerAngle',0.0000000001,'NoBoundaryRefinment',1}; +alreadyloaded=0; +%}}} +for i=sl.basinindx('basin','all'), + + bas=sl.basins{i}; + disp(sprintf('Meshing basin %s\n',bas.name)); + + %recover basin domain: + domain=bas.contour(); + + %recover coastline inside basin, using GSHHS_c_L1. It's a lat,long file, hence epsg 4326 + coastline=bas.shapefilecrop('shapefile',gshhsshapefile,'epsgshapefile',4326,'threshold',threshold); + + %mesh: + md=bamg(model,'domain',domain,'subdomains',coastline,'hmin',hmin,'hmax',hmax,defaultoptions{:}); + plotmodel(md,'data','mesh');pause(1); + error; + + %miscellaneous: + md.mesh.epsg=bas.epsg; md.miscellaneous.name=bas.name; + + %recover mask where we have land: + md.private.bamg.landmask=double(md.private.bamg.mesh.Triangles(:,4)>=1); + + %vertex connectivity: + md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); + + %add model to sl icecaps: + sl.addicecap(md); +end +%}}} +error; +%Parameterize ice sheets : {{{ + +for ind=sl.basinindx('continent',{'greenland','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; + + %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; + 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; + + %correction to land and ocean levelset: ice shelf is not on land! + pos=find(md.mask.ice_levelset<=0 & md.mask.groundedice_levelset<=0); + 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 + %}}} + %latlong: % {{{ + [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'),% {{{ + 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; + + end %}}} + % }}} + %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; +end +%}}} +% 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'}), + disp(sprintf('Masks for basin %s\n', sl.icecaps{ind}.miscellaneous.name)); + md=sl.icecaps{ind}; bas=sl.basins{ind}; + + %recover lat,long: + [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: + 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. + land=md.private.bamg.landmask; + land_mask=-ones(md.mesh.numberofvertices,1); + + landels=find(land); + land_mask(md.mesh.elements(landels,:))=1; + + %gothrough edges of each land element + connectedels=md.mesh.elementconnectivity(landels,:); + connectedisonocean=~land(connectedels); + sumconnectedisonocean=sum(connectedisonocean,2); + + %figure out which land elements are connected to the ocean: + landelsconocean=landels(find(sumconnectedisonocean)); + + ind1=[md.mesh.elements(landelsconocean,1); md.mesh.elements(landelsconocean,2); md.mesh.elements(landelsconocean,3)]; + 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)); + els=intersect(els1,els2); + + if length(find(land(els)))==1, + %this edge is on the beach, 0 the edge: + land_mask(ind1(i))=0; + land_mask(ind2(i))=0; + 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); + + %%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); + + %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; + % }}} + %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; + % }}} + sl.icecaps{ind}=md; +end +% }}} +%Assemble Earth in 3D {{{ + +%parameters: +plotting=1; +tolerance=100; +loneedgesdetect=0; + +%create earth model by concatenating all the icecaps in 3d: +sl.caticecaps('tolerance',tolerance,'loneedgesdetect',loneedgesdetect); + +%figure out how each icecap's mesh connects to the larger earth mesh: +sl.intersections(); + +%figure out connectivity: +disp('Mesh connectivity'); +sl.earth.mesh.vertexconnectivity=NodeConnectivity(sl.earth.mesh.elements,sl.earth.mesh.numberofvertices); + +%areas: +disp('Mesh nodal areas'); +sl.earth.mesh.area=averaging(sl.earth,GetAreas3DTria(sl.earth.mesh.elements,sl.earth.mesh.x,sl.earth.mesh.y,sl.earth.mesh.z),4); + +%transfer a list of fields from each icecap and continent back to Earth: +sl.transfer('mask.groundedice_levelset'); +sl.transfer('mask.ice_levelset'); +sl.transfer('mask.ocean_levelset'); +sl.transfer('mask.land_levelset'); +sl.transfer('mask.glacier_levelset'); +sl.transfer('geometry.thickness'); +sl.transfer('geometry.surface'); +sl.transfer('geometry.base'); +sl.transfer('geometry.bed'); +sl.transfer('initialization.temperature'); +sl.transfer('materials.rheology_B'); +sl.transfer('mesh.lat'); +sl.transfer('mesh.long'); +sl.transfer('slr.deltathickness'); +sl.transfer('slr.Ngia'); +sl.transfer('slr.Ugia'); + +%radius: +sl.earth.mesh.r=sqrt(sl.earth.mesh.x.^2+sl.earth.mesh.y.^2+sl.earth.mesh.z.^2); + +%check on the mesh transitions: {{{ +if plotting, + flags=ones(sl.earth.mesh.numberofvertices,1); + for i=1:length(sl.transitions), + flags(sl.transitions{i})=i; + end + plotmodel(sl.earth,'data',flags,'shading','faceted','coastline','on','coast_color','g') +end +%}}}} + +% }}} +%Solve Sea-level equation on Earth only: {{{ +md=sl.earth; %we don't do computations on ice sheets or land. + +%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)=[]; +md.slr.ocean_area_scaling = 0; + +%Miscellaneous +md.miscellaneous.name='test2004'; + +%New stuff +md.slr.spcthickness = NaN(md.mesh.numberofvertices,1); +md.slr.Ngia = zeros(md.mesh.numberofvertices,1); +md.slr.Ugia = zeros(md.mesh.numberofvertices,1); + +%Solution parameters +md.slr.reltol=NaN; +md.slr.abstol=1e-3; +md.slr.geodetic=1; + +%eustatic + rigid + elastic run: +md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=0; +md.cluster=generic('name',oshostname(),'np',3); +%md.verbose=verbose('111111111'); +md=solve(md,'Sealevelrise'); +SnoRotation=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; + +%Fields and tolerances to track changes +field_names ={'noRotation','Rotation'}; +field_tolerances={1e-13,1e-13}; +field_values={SnoRotation,SRotation}; + +%}}}