Index: /issm/trunk-jpl/test/NightlyRun/test2004.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 24906) +++ /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 24906) @@ -0,0 +1,408 @@ +%Test Name: Earth_Antarctica_GIA + +testagainst2002=0; + +%Data paths: {{{ +modeldatapath='/Users/larour/ModelData'; +shppath='/Users/larour/issm-jpl/proj-group/qgis/NightlyRun'; +gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1-NightlyRun.shp']; +%}}} + +%create sealevel model to hold our information: + sl=sealevelmodel(); + +%Create basins using boundaries from shapefile: +%some projections we'll rely on: %{{{ +proj4326=epsg2proj(4326); +proj3031=epsg2proj(3031); +%}}} +%HemisphereWest: {{{ +sl.addbasin(basin('continent','hemispherewest','name','hemispherewest','proj',laea(0,-90),'boundaries',{... %Peru projection 3587 + boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneBrunt','proj',proj3031,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031)... + })); + %}}} + %Ross: {{{ + sl.addbasin(basin('continent','antarctica','name','ross','proj',proj3031,'boundaries',{... + boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031,'orientation','reverse')... + })); +%}}} + %HemisphereEast: {{{ + sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','proj',laea(0,+90),'boundaries',{... %Australian projection lat,long + boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326),... + boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031)... + })); +%}}} +%Antarctica excluding Ronne: {{{ +sl.addbasin(basin('continent','antarctica','name','antarctica-grounded','proj',proj3031,'boundaries',{... + boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031,'orientation','reverse'),... + boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031)... + boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031)... + boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031)... + boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031)... + boundary('shppath',shppath,'shpfilename','RonneBrunt','proj',proj3031)... + })); + %}}} + %Ronne: {{{ + sl.addbasin(basin('continent','antarctica','name','ronne','proj',proj3031,'boundaries',{... + boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),... + boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse')... + })); +%}}} + +%Meshing +%Go through basins and mesh: %{{{ +%meshing parameters: {{{ +hmin=500; hmax=700; hmin=hmin*1000; hmax=hmax*1000; +tolerance=100; %tolerance of 100m on Earth position when mergin 3d meshes +threshold=5; +defaultoptions={'KeepVertices',0,'MaxCornerAngle',0.0000000001,'NoBoundaryRefinement',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); + + %miscellaneous: + md.mesh.proj=bas.proj; 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 +%}}} + +%Parameterization: +%Parameterize ice sheets : {{{ + +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}; + %masks : %{{{ + %ice levelset from domain outlines: + md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1); + + if bas.isnameany('antarctica-grounded'), + md.mask.ocean_levelset=ones(md.mesh.numberofvertices,1); + end + if bas.isnameany('ronne','ross'), + md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); + end + %}}} + %latlong: % {{{ + [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326'); + %}}} + %geometry: {{{ + if bas.iscontinentany('antarctica'), + di=md.materials.rho_ice/md.materials.rho_water; + + disp(' reading bedrock'); + md.geometry.bed=interpBedmap2(md.mesh.x,md.mesh.y,'bed'); + end % }}} + %Slr: {{{ + if bas.iscontinentany('antarctica'), + if testagainst2002, + md.slr.deltathickness=zeros(md.mesh.numberofelements,1); + %antarctica + late=sum(md.mesh.lat(md.mesh.elements),2)/3; + longe=sum(md.mesh.long(md.mesh.elements),2)/3; + pos=find(late <-85); + ratio=0.225314032985172/0.193045366574523; + %ratio= 1.276564103522540/.869956; + md.slr.deltathickness(pos)=-100*ratio; + else + delH=textread('../Data/AIS_delH_trend.txt'); + longAIS=delH(:,1); latAIS=delH(:,2); delHAIS=delH(:,3); index=delaunay(longAIS,latAIS); + lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360; + delHAIS=InterpFromMesh2d(index,longAIS,latAIS,delHAIS,long,lat); + northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0; + + md.slr.deltathickness=delHAIS(md.mesh.elements)*[1;1;1]/3/100; + end + + 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; + % }}} + + sl.icecaps{ind}=md; +end +%}}} +% ParameterizeContinents {{{ + +sl.basinindx('continent',{'hemisphereeast','hemispherewest'}) + +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}; + + %recover lat,long: + [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326'); + + %mask: %{{{ + %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)]; + + + %edge ind1 and ind2: + for i=1:length(ind1), + 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.ocean_levelset=land_mask; + md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); %if there are glaciers, we'll adjust + + if testagainst2002, + % {{{ + %greenland + pos=find(md.mesh.lat > 70 & md.mesh.lat < 80 & md.mesh.long>-60 & md.mesh.long<-30); + md.mask.ice_levelset(pos)=-1; + % }}} + end + + % }}} + %}}} + %slr loading/calibration: {{{ + + if testagainst2002, + % {{{ + md.slr.deltathickness=zeros(md.mesh.numberofelements,1); + %greenland + late=sum(md.mesh.lat(md.mesh.elements),2)/3; + longe=sum(md.mesh.long(md.mesh.elements),2)/3; + pos=find(late > 70 & late < 80 & longe>-60 & longe<-30); + ratio=.3823/.262344; + %md.slr.deltathickness(pos)=-100*ratio; + + %correct mask: + md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; + % }}} + else + + md.slr.deltathickness=zeros(md.mesh.numberofelements,1); + + delH=textread('../Data/GIS_delH_trend.txt'); + longGIS=delH(:,1); latGIS=delH(:,2); delHGIS=delH(:,3); index=delaunay(longGIS,latGIS); + lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360; + delHGIS=InterpFromMeshToMesh2d(index,longGIS,latGIS,delHGIS,long,lat); + delHGISe=delHGIS(md.mesh.elements)*[1;1;1]/3; + + delH=textread('../Data/GLA_delH_trend_15regions.txt'); + longGLA=delH(:,1); latGLA=delH(:,2); delHGLA=sum(delH(:,3:end),2); index=delaunay(longGLA,latGLA); + lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360; + delHGLA=InterpFromMeshToMesh2d(index,longGLA,latGLA,delHGLA,long,lat); + delHGLAe=delHGLA(md.mesh.elements)*[1;1;1]/3; + + pos=find(delHGISe); + md.slr.deltathickness(pos)=delHGISe(pos)/100; + pos=find(delHGLAe); + md.slr.deltathickness(pos)=delHGLAe(pos)/100; + + %adjust mask accordingly: + pos=find(md.slr.deltathickness); + flags=zeros(md.mesh.numberofvertices,1); + flags(md.mesh.elements(pos,:))=1; + pos=find(flags); + md.mask.ice_levelset(pos)=-1; + md.mask.ocean_levelset(pos)=1; + end + + 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); + + %}}} + %geometry: {{{ + di=md.materials.rho_ice/md.materials.rho_water; + md.geometry.bed=-ones(md.mesh.numberofvertices,1); + % }}} + %materials: {{{ + md.materials=materials('hydro'); + % }}} + 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('force',1); + +%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.ice_levelset'); +sl.transfer('mask.ocean_levelset'); +sl.transfer('geometry.bed'); +sl.transfer('mesh.lat'); +sl.transfer('mesh.long'); +sl.transfer('slr.deltathickness'); +sl.transfer('slr.spcthickness'); +sl.transfer('slr.Ngia'); +sl.transfer('slr.Ugia'); +sl.transfer('slr.hydro_rate'); +sl.transfer('slr.sealevel'); +sl.transfer('dsl.sea_surface_height_change_above_geoid'); +sl.transfer('dsl.sea_water_pressure_change_at_sea_floor'); + +%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: {{{ +plotting=1; +if plotting, + flags=ones(sl.earth.mesh.numberofelements,1); + for i=1:length(sl.eltransitions), + flags(sl.eltransitions{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. + +%Materials: +md.materials=materials('hydro'); + +%elastic loading from love numbers: +nlov=101; +md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; +md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; +md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; +md.slr.ocean_area_scaling = 0; +md.slr.loop_increment=200; + +%Miscellaneous +md.miscellaneous.name='test2004'; + +%New stuff +md.dsl.global_average_thermosteric_sea_level_change=[1.1+.38;0]; %steric + water storage AR5. + +%Solution parameters +md.slr.reltol=NaN; +md.slr.abstol=1e-3; +md.slr.geodetic=1; +md.timestepping.time_step=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.slr.maxiter=10; + +%eustatic run: +md.slr.rigid=0; md.slr.elastic=0;md.slr.rotation=0; +md.slr.requested_outputs= {'default',... + 'SealevelriseDeltathickness','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... + 'SealevelNEsaRate', 'SealevelUEsaRate', 'SealevelNGiaRate', 'SealevelUGiaRate','SealevelEustaticMask','SealevelEustaticOceanMask'}; +md=solve(md,'Sealevelrise'); +Seustatic=md.results.SealevelriseSolution.Sealevel; + +%eustatic + rigid run: +md.slr.rigid=1; md.slr.elastic=0;md.slr.rotation=0; +md=solve(md,'Sealevelrise'); +Srigid=md.results.SealevelriseSolution.Sealevel; + +%eustatic + rigid + elastic run: +md.slr.rigid=1; md.slr.elastic=1;md.slr.rotation=0; +md=solve(md,'Sealevelrise'); +Selastic=md.results.SealevelriseSolution.Sealevel; + +%eustatic + rigid + elastic + rotation run: +md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=1; +md=solve(md,'Sealevelrise'); +Srotation=md.results.SealevelriseSolution.Sealevel; + +%}}} Index: sm/trunk-jpl/test/SandBox/test2004.m =================================================================== --- /issm/trunk-jpl/test/SandBox/test2004.m (revision 24905) +++ (revision ) @@ -1,408 +1,0 @@ -%Test Name: Earth_Antarctica_GIA - -testagainst2002=0; - -%Data paths: {{{ -modeldatapath='/Users/larour/ModelData'; -shppath='/Users/larour/issm-jpl/proj-group/qgis/NightlyRun'; -gshhsshapefile=[modeldatapath '/Gshhg/Shp/GSHHS_shp/c/GSHHS_c_L1-NightlyRun.shp']; -%}}} - -%create sealevel model to hold our information: - sl=sealevelmodel(); - -%Create basins using boundaries from shapefile: -%some projections we'll rely on: %{{{ -proj4326=epsg2proj(4326); -proj3031=epsg2proj(3031); -%}}} -%HemisphereWest: {{{ -sl.addbasin(basin('continent','hemispherewest','name','hemispherewest','proj',laea(0,-90),'boundaries',{... %Peru projection 3587 - boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RonneBrunt','proj',proj3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031)... - })); - %}}} - %Ross: {{{ - sl.addbasin(basin('continent','antarctica','name','ross','proj',proj3031,'boundaries',{... - boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031,'orientation','reverse')... - })); -%}}} - %HemisphereEast: {{{ - sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','proj',laea(0,+90),'boundaries',{... %Australian projection lat,long - boundary('shppath',shppath,'shpfilename','HemisphereSplit','proj',proj4326),... - boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossFront','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031)... - })); -%}}} -%Antarctica excluding Ronne: {{{ -sl.addbasin(basin('continent','antarctica','name','antarctica-grounded','proj',proj3031,'boundaries',{... - boundary('shppath',shppath,'shpfilename','NorthAntarctica','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','EastAntarctica2','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossWestFront','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RossIceShelf','proj',proj3031,'orientation','reverse'),... - boundary('shppath',shppath,'shpfilename','SouthAntarctica','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','WestAntarctica2','proj',proj3031)... - boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031)... - boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031)... - boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031)... - boundary('shppath',shppath,'shpfilename','RonneBrunt','proj',proj3031)... - })); - %}}} - %Ronne: {{{ - sl.addbasin(basin('continent','antarctica','name','ronne','proj',proj3031,'boundaries',{... - boundary('shppath',shppath,'shpfilename','RonneWestSummit','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RonneIceShelf','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RonneEastSummit','proj',proj3031),... - boundary('shppath',shppath,'shpfilename','RonneFront','proj',proj3031,'orientation','reverse')... - })); -%}}} - -%Meshing -%Go through basins and mesh: %{{{ -%meshing parameters: {{{ -hmin=500; hmax=700; hmin=hmin*1000; hmax=hmax*1000; -tolerance=100; %tolerance of 100m on Earth position when mergin 3d meshes -threshold=5; -defaultoptions={'KeepVertices',0,'MaxCornerAngle',0.0000000001,'NoBoundaryRefinement',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); - - %miscellaneous: - md.mesh.proj=bas.proj; 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 -%}}} - -%Parameterization: -%Parameterize ice sheets : {{{ - -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}; - %masks : %{{{ - %ice levelset from domain outlines: - md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1); - - if bas.isnameany('antarctica-grounded'), - md.mask.ocean_levelset=ones(md.mesh.numberofvertices,1); - end - if bas.isnameany('ronne','ross'), - md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1); - end - %}}} - %latlong: % {{{ - [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326'); - %}}} - %geometry: {{{ - if bas.iscontinentany('antarctica'), - di=md.materials.rho_ice/md.materials.rho_water; - - disp(' reading bedrock'); - md.geometry.bed=interpBedmap2(md.mesh.x,md.mesh.y,'bed'); - end % }}} - %Slr: {{{ - if bas.iscontinentany('antarctica'), - if testagainst2002, - md.slr.deltathickness=zeros(md.mesh.numberofelements,1); - %antarctica - late=sum(md.mesh.lat(md.mesh.elements),2)/3; - longe=sum(md.mesh.long(md.mesh.elements),2)/3; - pos=find(late <-85); - ratio=0.225314032985172/0.193045366574523; - %ratio= 1.276564103522540/.869956; - md.slr.deltathickness(pos)=-100*ratio; - else - delH=textread('../Data/AIS_delH_trend.txt'); - longAIS=delH(:,1); latAIS=delH(:,2); delHAIS=delH(:,3); index=delaunay(longAIS,latAIS); - lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360; - delHAIS=InterpFromMesh2d(index,longAIS,latAIS,delHAIS,long,lat); - northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0; - - md.slr.deltathickness=delHAIS(md.mesh.elements)*[1;1;1]/3/100; - end - - 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; - % }}} - - sl.icecaps{ind}=md; -end -%}}} -% ParameterizeContinents {{{ - -sl.basinindx('continent',{'hemisphereeast','hemispherewest'}) - -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}; - - %recover lat,long: - [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326'); - - %mask: %{{{ - %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)]; - - - %edge ind1 and ind2: - for i=1:length(ind1), - 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.ocean_levelset=land_mask; - md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); %if there are glaciers, we'll adjust - - if testagainst2002, - % {{{ - %greenland - pos=find(md.mesh.lat > 70 & md.mesh.lat < 80 & md.mesh.long>-60 & md.mesh.long<-30); - md.mask.ice_levelset(pos)=-1; - % }}} - end - - % }}} - %}}} - %slr loading/calibration: {{{ - - if testagainst2002, - % {{{ - md.slr.deltathickness=zeros(md.mesh.numberofelements,1); - %greenland - late=sum(md.mesh.lat(md.mesh.elements),2)/3; - longe=sum(md.mesh.long(md.mesh.elements),2)/3; - pos=find(late > 70 & late < 80 & longe>-60 & longe<-30); - ratio=.3823/.262344; - %md.slr.deltathickness(pos)=-100*ratio; - - %correct mask: - md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; - % }}} - else - - md.slr.deltathickness=zeros(md.mesh.numberofelements,1); - - delH=textread('../Data/GIS_delH_trend.txt'); - longGIS=delH(:,1); latGIS=delH(:,2); delHGIS=delH(:,3); index=delaunay(longGIS,latGIS); - lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360; - delHGIS=InterpFromMeshToMesh2d(index,longGIS,latGIS,delHGIS,long,lat); - delHGISe=delHGIS(md.mesh.elements)*[1;1;1]/3; - - delH=textread('../Data/GLA_delH_trend_15regions.txt'); - longGLA=delH(:,1); latGLA=delH(:,2); delHGLA=sum(delH(:,3:end),2); index=delaunay(longGLA,latGLA); - lat=md.mesh.lat; long=md.mesh.long+360; pos=find(long>360);long(pos)=long(pos)-360; - delHGLA=InterpFromMeshToMesh2d(index,longGLA,latGLA,delHGLA,long,lat); - delHGLAe=delHGLA(md.mesh.elements)*[1;1;1]/3; - - pos=find(delHGISe); - md.slr.deltathickness(pos)=delHGISe(pos)/100; - pos=find(delHGLAe); - md.slr.deltathickness(pos)=delHGLAe(pos)/100; - - %adjust mask accordingly: - pos=find(md.slr.deltathickness); - flags=zeros(md.mesh.numberofvertices,1); - flags(md.mesh.elements(pos,:))=1; - pos=find(flags); - md.mask.ice_levelset(pos)=-1; - md.mask.ocean_levelset(pos)=1; - end - - 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); - - %}}} - %geometry: {{{ - di=md.materials.rho_ice/md.materials.rho_water; - md.geometry.bed=-ones(md.mesh.numberofvertices,1); - % }}} - %materials: {{{ - md.materials=materials('hydro'); - % }}} - 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('force',1); - -%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.ice_levelset'); -sl.transfer('mask.ocean_levelset'); -sl.transfer('geometry.bed'); -sl.transfer('mesh.lat'); -sl.transfer('mesh.long'); -sl.transfer('slr.deltathickness'); -sl.transfer('slr.spcthickness'); -sl.transfer('slr.Ngia'); -sl.transfer('slr.Ugia'); -sl.transfer('slr.hydro_rate'); -sl.transfer('slr.sealevel'); -sl.transfer('dsl.sea_surface_height_change_above_geoid'); -sl.transfer('dsl.sea_water_pressure_change_at_sea_floor'); - -%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: {{{ -plotting=1; -if plotting, - flags=ones(sl.earth.mesh.numberofelements,1); - for i=1:length(sl.eltransitions), - flags(sl.eltransitions{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. - -%Materials: -md.materials=materials('hydro'); - -%elastic loading from love numbers: -nlov=101; -md.slr.love_h = love_numbers('h','CM'); md.slr.love_h(nlov+1:end)=[]; -md.slr.love_k = love_numbers('k','CM'); md.slr.love_k(nlov+1:end)=[]; -md.slr.love_l = love_numbers('l','CM'); md.slr.love_l(nlov+1:end)=[]; -md.slr.ocean_area_scaling = 0; -md.slr.loop_increment=200; - -%Miscellaneous -md.miscellaneous.name='test2004'; - -%New stuff -md.dsl.global_average_thermosteric_sea_level_change=[1.1+.38;0]; %steric + water storage AR5. - -%Solution parameters -md.slr.reltol=NaN; -md.slr.abstol=1e-3; -md.slr.geodetic=1; -md.timestepping.time_step=1; - -% max number of iteration reverted back to 10 (i.e., the original default value) -md.slr.maxiter=10; - -%eustatic run: -md.slr.rigid=0; md.slr.elastic=0;md.slr.rotation=0; -md.slr.requested_outputs= {'default',... - 'SealevelriseDeltathickness','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... - 'SealevelNEsaRate', 'SealevelUEsaRate', 'SealevelNGiaRate', 'SealevelUGiaRate','SealevelEustaticMask','SealevelEustaticOceanMask'}; -md=solve(md,'Sealevelrise'); -Seustatic=md.results.SealevelriseSolution.Sealevel; - -%eustatic + rigid run: -md.slr.rigid=1; md.slr.elastic=0;md.slr.rotation=0; -md=solve(md,'Sealevelrise'); -Srigid=md.results.SealevelriseSolution.Sealevel; - -%eustatic + rigid + elastic run: -md.slr.rigid=1; md.slr.elastic=1;md.slr.rotation=0; -md=solve(md,'Sealevelrise'); -Selastic=md.results.SealevelriseSolution.Sealevel; - -%eustatic + rigid + elastic + rotation run: -md.slr.rigid=1; md.slr.elastic=1; md.slr.rotation=1; -md=solve(md,'Sealevelrise'); -Srotation=md.results.SealevelriseSolution.Sealevel; - -%}}}