source: issm/trunk-jpl/test/SandBox/test2004.m@ 24763

%Test Name: Earth_Antarctica_GIA

%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.shp'];
%}}}

skipmesh=1;

%create sealevel model to hold our information: 
if ~skipmesh,
        sl=sealevelmodel();
%Create basins using boundaries from shapefile: %{{{
%HemisphereWest: {{{
        sl.addbasin(basin('continent','hemispherewest','name','hemispherewest','epsg',3587,'boundaries',{... %Peru projection 3587
        boundary('shppath',shppath,'shpfilename','HemisphereSplit','epsg',4326,'orientation','reverse'),...
        boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031),...
        boundary('shppath',shppath,'shpfilename','RonneBrunt','epsg',3031,'orientation','reverse'),...
        boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031),...
        boundary('shppath',shppath,'shpfilename','RonneFront','epsg',3031,'orientation','reverse'),...
        boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031),...
        boundary('shppath',shppath,'shpfilename','WestAntarctica2','epsg',3031,'orientation','reverse'),...
        boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031)...
        }));
        %}}}
        %HemisphereEast: {{{
        sl.addbasin(basin('continent','hemisphereeast','name','hemisphereeast','epsg',4462,'boundaries',{... %Australian projection lat,long
        boundary('shppath',shppath,'shpfilename','HemisphereSplit','epsg',4326),...
        boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031),...
        boundary('shppath',shppath,'shpfilename','EastAntarctica2','epsg',3031,'orientation','reverse'),...
        boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031)...
        }));
        %}}}
        %Antarctica excluding Ronne: {{{
        sl.addbasin(basin('continent','antarctica','name','antarctica-grounded','epsg',3031,'boundaries',{...
        boundary('shppath',shppath,'shpfilename','NorthAntarctica','epsg',3031),...
        boundary('shppath',shppath,'shpfilename','EastAntarctica2','epsg',3031),...
        boundary('shppath',shppath,'shpfilename','SouthAntarctica','epsg',3031),...
        boundary('shppath',shppath,'shpfilename','WestAntarctica2','epsg',3031)...
        boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031)...
        boundary('shppath',shppath,'shpfilename','RonneIceShelf','epsg',3031)...
        boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031)...
        boundary('shppath',shppath,'shpfilename','RonneBrunt','epsg',3031)...
        }));
        %}}}
        %Ronne: {{{
        sl.addbasin(basin('continent','antarctica','name','ronne','epsg',3031,'boundaries',{...
        boundary('shppath',shppath,'shpfilename','RonneWestSummit','epsg',3031),...
        boundary('shppath',shppath,'shpfilename','RonneIceShelf','epsg',3031),...
        boundary('shppath',shppath,'shpfilename','RonneEastSummit','epsg',3031),...
        boundary('shppath',shppath,'shpfilename','RonneFront','epsg',3031,'orientation','reverse')...
        }));
        %}}}
%}}}
%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);

        %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
%}}}
end
%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 :  %{{{
        %new type of mask: 
        md.mask=maskpsl; 

        %land and ocean: 
        md.mask.land_levelset=ones(md.mesh.numberofvertices,1);
        md.mask.ocean_levelset=zeros(md.mesh.numberofvertices,1);

        %ice levelset from domain outlines: 
        md.mask.ice_levelset=-ones(md.mesh.numberofvertices,1);
        pos=find(md.mesh.vertexonboundary); md.mask.ice_levelset(pos)=0;

        %ice front: 
        pos=find(md.mask.ice_levelset==0); md.mask.ocean_levelset(pos)=-1; md.mask.land_levelset(pos)=1;

        %now, glaciers: 
        md.mask.glacier_levelset=zeros(md.mesh.numberofvertices,1);
        %}}} 
        %initial grounding line position: {{{
        if bas.isnameany('antarctica-grounded'), 

                md.mask.groundedice_levelset=ones(md.mesh.numberofvertices,1);
        end
        if bas.isnameany('ronne'), 
                
                md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1);

                %correction to land and ocean levelset: ice shelf is not on land! 
                pos=find(md.mask.ice_levelset<=0 & md.mask.groundedice_levelset<=0);
                md.mask.ocean_levelset(pos)=1;
                md.mask.land_levelset(pos)=-1;
        end
        %}}}
        %latlong:  % {{{
        [md.mesh.long,md.mesh.lat]=gdaltransform(md.mesh.x,md.mesh.y,sprintf('EPSG:%i',md.mesh.epsg),'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'),

                delH=textread([modeldatapath '/AdhikariSciAdvTrends/AIS_delH_trend.txt']);
                longAIS=delH(:,1);
                latAIS=delH(:,2);
                delHAIS=delH(:,3);
                index=delaunay(latAIS,longAIS);

                lat=md.mesh.lat; 
                long=md.mesh.long+360; 
                pos=find(long>360);long(pos)=long(pos)-360;

                delHAIS=InterpFromMesh2d(index,longAIS,latAIS,delHAIS,long,lat);
                northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0;

                md.slr.deltathickness=delHAIS(md.mesh.elements)*[1;1;1]/3/100;

                md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
                md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
                md.slr.Ngia=zeros(md.mesh.numberofvertices,1);
                md.slr.Ugia=zeros(md.mesh.numberofvertices,1);

                md.dsl.global_average_thermosteric_sea_level_change=[0;0];
                md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
                md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
                md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1);

        end %}}}
        % material properties: {{{
        md.materials=materials('hydro');
        %}}}
        %diverse: {{{
        md.miscellaneous.name=bas.name;
        % }}}

        sl.icecaps{ind}=md;
end
%}}}
% ParameterizeContinents {{{

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,sprintf('EPSG:%i',md.mesh.epsg),'EPSG:4326'); 

        %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); 

        %grounded ice: 
        md.mask.groundedice_levelset=-ones(md.mesh.numberofvertices,1); 

        % }}}
        %}}}
        %slr loading/calibration:  {{{

        md.slr.sealevel=zeros(md.mesh.numberofvertices,1);
        md.slr.spcthickness=NaN*ones(md.mesh.numberofvertices,1);
        md.slr.Ngia=zeros(md.mesh.numberofvertices,1);
        md.slr.Ugia=zeros(md.mesh.numberofvertices,1);

        md.dsl.global_average_thermosteric_sea_level_change=[0;0];
        %md.slr.steric_rate=(1.1+.38)*ones(md.mesh.numberofvertices,1); %steric + water storage.
        md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
        md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
        md.slr.hydro_rate = zeros(md.mesh.numberofvertices,1);

        %love numbers: 
        md.slr.ocean_area_scaling=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);

error;

%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
%}}}}

% }}}
error;
%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};

%}}}