Index: /issm/trunk-jpl/test/SandBox/test2004.m =================================================================== --- /issm/trunk-jpl/test/SandBox/test2004.m (revision 24774) +++ /issm/trunk-jpl/test/SandBox/test2004.m (revision 24775) @@ -7,57 +7,59 @@ %}}} -skipmesh=1; -skipparameterization=0; - %create sealevel model to hold our information: -if ~skipmesh, sl=sealevelmodel(); -%Create basins using boundaries from shapefile: %{{{ + +%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','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)... +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)... })); %}}} %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)... + 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','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','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)... })); %}}} - %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')... - })); - %}}} -%}}} + 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=2000; hmin=hmin*1000; hmax=hmax*1000; +hmin=2000; hmax=10000; hmin=hmin*1000; hmax=hmax*1000; tolerance=100; %tolerance of 100m on Earth position when mergin 3d meshes -threshold=1; +threshold=5; defaultoptions={'KeepVertices',0,'MaxCornerAngle',0.0000000001,'NoBoundaryRefinment',1}; alreadyloaded=0; @@ -79,5 +81,5 @@ %miscellaneous: - md.mesh.epsg=bas.epsg; md.miscellaneous.name=bas.name; + md.mesh.proj=bas.proj; md.miscellaneous.name=bas.name; %recover mask where we have land: @@ -91,6 +93,6 @@ end %}}} -end -if ~skipparameterization, + +%Parameterization: %Parameterize ice sheets : {{{ @@ -123,5 +125,5 @@ end if bas.isnameany('ronne'), - + md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1); @@ -133,5 +135,5 @@ %}}} %latlong: % {{{ - [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326'); + [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326'); %}}} %geometry: {{{ @@ -190,5 +192,5 @@ %recover lat,long: - [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326'); + [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,md.mesh.proj,'EPSG:4326'); %mask: %{{{ @@ -244,5 +246,5 @@ %grounded ice: md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1); - + md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1); @@ -252,4 +254,5 @@ md.slr.sealevel=zeros(md.mesh.numberofvertices,1); + md.slr.deltathickness=zeros(md.mesh.numberofelements,1); md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1); md.slr.Ngia=zeros(md.mesh.numberofvertices,1); @@ -275,6 +278,6 @@ end % }}} -end -%Assemble Earth in 3D {{{ + +%%Assemble Earth in 3D {{{ %parameters: @@ -287,5 +290,5 @@ %figure out how each icecap's mesh connects to the larger earth mesh: -sl.intersections(); +sl.intersections('force',1); %figure out connectivity: @@ -307,7 +310,11 @@ sl.transfer('mesh.long'); sl.transfer('slr.deltathickness'); +sl.transfer('slr.spcthickness'); sl.transfer('slr.Ngia'); sl.transfer('slr.Ugia'); -error; +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: @@ -315,8 +322,9 @@ %check on the mesh transitions: {{{ +plotting=1; if plotting, - flags=ones(sl.earth.mesh.numberofvertices,1); - for i=1:length(sl.transitions), - flags(sl.transitions{i})=i; + 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') @@ -325,7 +333,10 @@ % }}} -error; + %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: @@ -340,7 +351,5 @@ %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); +md.dsl.global_average_thermosteric_sea_level_change=[0;0]; %Solution parameters @@ -349,8 +358,13 @@ md.slr.geodetic=1; +%eustatic: +md.slr.rigid=0; md.slr.elastic=0; md.slr.rotation=0; +md.cluster=generic('name',oshostname(),'np',10); +md=solve(md,'Sealevelrise'); + %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.cluster=generic('name',oshostname(),'np',16); +md.verbose=verbose('111111111'); md=solve(md,'Sealevelrise'); SnoRotation=md.results.SealevelriseSolution.Sealevel;