% Script written by Eliza Dawson to generate figures for Dawson et al.,
% 2022, "Ice mass loss sensitivity to the Antarctic ice sheet basal thermal state"

% To run this script, ISSM needs to be downloaded and installed
% (https://issm.jpl.nasa.gov/download/)

% Note that this script also uses ColorBrewer colormaps which can be downloaded
% here: https://www.mathworks.com/matlabcentral/fileexchange/45208-colorbrewer-attractive-and-distinctive-colormaps?s_tid=mwa_osa_a

% Settings: 
clear all; close all; clc
step=[1];
glacier='Antarctica';
rootdir='./'; 
org=organizer('prefix', [glacier '_'], 'steps', step);


% Step 1: Plot Antarctic Ice Sheet Basal Temperatures {{{
if perform(org, 'BasalTemps');
    md=loadmodel(org,'Initialization');
    T = md.results.ThermalSolution.Temperature;
    Tm = T + md.materials.beta.*md.initialization.pressure; % pressure adjusted temperature        
        
    plotmodel(md,'data',Tm,...
        'caxis',[273.5-19,273.5],'xlabel','m','ylabel','m','colorbartitle','K',...
        'title','Pressure Adjusted Basal Temperatures','axis','off','layer',1,'figure',1);
    c1 = colormap(brewermap(9,'*Blues'));
    c1 = c1(1:6,:);
    c2 = colormap(brewermap(22,'YlOrRd'));
    c3 = c2(end-3,:);
    c2 = c2(1:12,:);
    c = [c1 ; c2; c3];
    colorbar;
    colormap(c);
end


% Step 2: Spatial Plots of Experiment Zone Temperature Ranges {{{
if perform(org, 'TempMasks');
    md=loadmodel(org,'Initialization');
    T = md.results.ThermalSolution.Temperature;
    Tm = T + md.materials.beta.*md.initialization.pressure; % pressure adjusted temperature
    FricCoef = sqrt(md.friction.coefficient.^2 / md.constants.yts);
    FricCoef(md.mask.groundedice_levelset<0)=NaN;
    Tm(md.mask.groundedice_levelset<0)=NaN;
    
    load('Masks/BedVelMask_100_km_Of_100_m_per_yr.mat')

    Thw = 273.14; 
    Frz = 271.14; % Change the frozen threshold to change the frozen/thawable mask
    
    FricMasked = BedVelMask.*FricCoef;
    FricMasked(FricMasked==0)= NaN;
    TmMasked = BedVelMask.*Tm;
    TmMasked(TmMasked==0)= NaN;

    Tb = project2d(md,TmMasked,1);
    BFric = project2d(md,FricMasked,1);

    indicesT = find(Tb>Thw);
    indicesF = find(Tb<Frz);
    indicesFtoT = find(Tb<Thw & Tb>Frz);
    TbMask = NaN(length(md.mesh.x2d),1);
    TbMask(indicesF) = 1;
    TbMask(indicesFtoT) = 2;
    TbMask(indicesT) = 3;
    TbMask(md.mask.groundedice_levelset<0)=NaN;
    
    basin = {'DrainageBasins/Dp-E.exp','DrainageBasins/Cp-D.exp','DrainageBasins/A-Ap.exp','DrainageBasins/B-C.exp',...
    'DrainageBasins/Ap-B.exp','DrainageBasins/D-Dp.exp','DrainageBasins/G-H.exp','DrainageBasins/Ep-F.exp',...
    'DrainageBasins/Hp-I.exp','DrainageBasins/I-Ipp.exp','DrainageBasins/Ipp-J.exp','DrainageBasins/J-Jpp.exp',...
    'DrainageBasins/E-Ep.exp','DrainageBasins/F-G.exp','DrainageBasins/H-Hp.exp','DrainageBasins/K-A.exp',...
    'DrainageBasins/C-Cp.exp','DrainageBasins/Jpp-K.exp'};

    plotmodel(md,'data',TbMask,'layer',1,...
        'expdisp',{basin{1},basin{2},basin{3},basin{4},basin{5},basin{6},basin{7},basin{8},basin{9},basin{10},...
         basin{11},basin{12},basin{13},basin{14},basin{15},basin{16},basin{17},basin{18}},'linewidth',2,'expstyle',...
         {[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],...
         [0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],...
         [0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7]},...
         'axis','off','title', [num2str(Thw-Frz) 'deg thawable'])
    c1 = colormap(brewermap(12,'*Blues'));
    c1 = c1(3,:);
    c2 = colormap(brewermap(16,'YlOrRd'));
    c3 = c2(end-3,:);
    c2 = c2(5,:);
    c = [c1 ; c2; c3];
    colorbar;
    colormap(c);
    name = (['BTempMask' num2str(Thw-Frz) 'deg.eps']);   
end


% Step 3: Spatial Plots of Friction Coefficient Changes {{{
if perform(org, 'CompareFricReductions');

    md=loadmodel(org,'Transient_Exp4_100yrs'); % Can switch this to Exp1, Exp2, or Exp3 instead
    md2=loadmodel(org,'Transient_Ctrl_100yrs'); 
    
    FricExp = sqrt(md.friction.coefficient.^2 / md.constants.yts);
    FricCtrl = sqrt(md2.friction.coefficient.^2 / md2.constants.yts);
    
    basin = {'DrainageBasins/Dp-E.exp','DrainageBasins/Cp-D.exp','DrainageBasins/A-Ap.exp','DrainageBasins/B-C.exp',...
    'DrainageBasins/Ap-B.exp','DrainageBasins/D-Dp.exp','DrainageBasins/G-H.exp','DrainageBasins/Ep-F.exp',...
    'DrainageBasins/Hp-I.exp','DrainageBasins/I-Ipp.exp','DrainageBasins/Ipp-J.exp','DrainageBasins/J-Jpp.exp',...
    'DrainageBasins/E-Ep.exp','DrainageBasins/F-G.exp','DrainageBasins/H-Hp.exp','DrainageBasins/K-A.exp',...
    'DrainageBasins/C-Cp.exp','DrainageBasins/Jpp-K.exp'};

    plotmodel(md,'data',FricExp - FricCtrl,'layer',1,'figure',4,...
        'caxis',[-50 50],'colormap', colormap(brewermap([],'PuOr')),...
        'expdisp',{basin{1},basin{2},basin{3},basin{4},basin{5},basin{6},basin{7},basin{8},basin{9},basin{10},...
         basin{11},basin{12},basin{13},basin{14},basin{15},basin{16},basin{17},basin{18}},'linewidth',2,'expstyle',...
         {[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],...
         [0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],...
         [0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7],[0.7 0.7 0.7]},...
         'axis','off','title','Fric Coef Change (Exp4-Ctrl)')
end


% Step 4: Plots of CDFs for Friction Coefficient by Drainage Basin
if perform(org, 'FricCoefCDFs');
    md=loadmodel(org,'Initialization');
    T = md.results.ThermalSolution.Temperature;
    Tm = T + md.materials.beta.*md.initialization.pressure; % pressure adjusted temperature
    FricCoef = sqrt(md.friction.coefficient.^2 / md.constants.yts);
%     FricCoef = md.friction.coefficient;
    FricCoef(md.mask.groundedice_levelset<0)=NaN;
    Tm(md.mask.groundedice_levelset<0)=NaN;

    Thw = 273.14;
    Frz = 271.14;
    id ={'Dp-E','Cp-D','A-Ap','B-C','Ap-B','D-Dp','E-Ep','C-Cp','Jpp-K','K-A',... % E. Ant (10)
        'G-H','Ep-F','J-Jpp','F-G','H-Hp',... % W. Ant (5)
        'Hp-I','I-Ipp','Ipp-J'}; % Ant. Peninsula (3)

    load('Masks/BedVelMask_100_km_Of_100_m_per_yr.mat')

    for ii = 1:length(id)
        basin = ['DrainageBasins/' id{ii} '.exp'];
        
        distance=ExpToLevelSet(md.mesh.x2d,md.mesh.y2d,basin); %in meters
        distance(md.mask.groundedice_levelset<0)=NaN;
        pos = find(distance>0);
        basin_mask = zeros(md.mesh.numberofvertices,1);
        basin_mask(pos)=1;

        FricMasked = BedVelMask.*basin_mask.*FricCoef;
        FricMasked(FricMasked==0)= NaN;
        TmMasked = BedVelMask.*basin_mask.*Tm;
        TmMasked(TmMasked==0)= NaN;

        Tb = project2d(md,TmMasked,1);
        BFric = project2d(md,FricMasked,1);

        indicesT = find(Tb>Thw);
        indicesF = find(Tb<Frz);
        indicesFtoT = find(Tb<Thw & Tb>Frz);

        BFricCoefTh = BFric(indicesT);
        BFricCoefFz = BFric(indicesF);
        BFricCoefFtoT = BFric(indicesFtoT);
        
        figure(ii)
        hold on;
        [h,stats] = cdfplot(BFricCoefTh);
        plot(h.XData,h.YData,'r','LineWidth',2)
        [h,stats] = cdfplot(BFricCoefFtoT);
        plot(h.XData,h.YData,'color',[0.9 0.6 0],'LineWidth',2)
        [h,stats] = cdfplot(BFricCoefFz);
        plot(h.XData,h.YData,'color',[0.4 0.4 1],'LineWidth',2)
        xlabel('Friction Coefficient')
        title(['Empirical CDF Basin ' id{ii}]);
        set(gca,'FontSize',16,'LineWidth',1)
        hold off; 
    end
end


% Step 5: Plots of Regional Ice Sheet Mass Loss {{{
if perform(org, 'RegionalIceSheetMassLoss');
    
        model={'Transient_Ctrl_100yrs','Transient_Exp1_100yrs','Transient_Exp2_100yrs','Transient_Exp3_100yrs','Transient_Exp4_100yrs'};
        id ={'Dp-E','Cp-D','A-Ap','B-C','Ap-B','D-Dp','E-Ep','C-Cp','Jpp-K','K-A',... % E. Ant (10)
            'G-H','Ep-F','J-Jpp','F-G','H-Hp',... % W. Ant (5)
            'Hp-I','I-Ipp','Ipp-J'}; % Ant. Peninsula (3)
        yr=51;
        
        %%%%%%%%%%%%% calculate volume above flotation for each experiment
        
        clear meanbed; clear meanthickness; clear md
        for ii = 1
            md=loadmodel(org,model{ii});
            for jj = 1:yr
                for kk = 1:length(id)
                    Bed = md.results.TransientSolution(jj).Bed;
                    Thickness =  md.results.TransientSolution(jj).Thickness;
                    time(jj) =  md.results.TransientSolution(jj).time;

                    basin = ['DrainageBasins/' id{kk} '.exp'];
                    distance=ExpToLevelSet(md.mesh.x2d,md.mesh.y2d,basin); %in meters
                    distance(md.mask.groundedice_levelset<0)=NaN;
                    pos = find(distance>0);
                    basin_mask = zeros(md.mesh.numberofvertices,1);
                    basin_mask(pos)=1;
                    Thickness_basin = Thickness.*basin_mask;
                    Bed_basin = Bed.*basin_mask;

                    meanbed=mean(Bed_basin(md.mesh.elements2d),2);
                    meanthickness=mean(Thickness_basin(md.mesh.elements2d),2);

                    H_AF=meanthickness + min(md.materials.rho_water/md.materials.rho_ice.*meanbed,0);
                    H_AF(H_AF<0)=0;
                    areas=GetAreas(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d);
                    vol_AF_basin = areas.*H_AF;
                    vol_all_ctrl(kk,jj) = sum(vol_AF_basin);
                end
            end
        end

        clear meanbed; clear meanthickness; clear md
        for ii = 2
            md=loadmodel(org,model{ii});
            for jj = 1:yr
                for kk = 1:length(id)
                    Bed = md.results.TransientSolution(jj).Bed;
                    Thickness =  md.results.TransientSolution(jj).Thickness;

                    basin = ['DrainageBasins/' id{kk} '.exp'];
                    distance=ExpToLevelSet(md.mesh.x2d,md.mesh.y2d,basin); %in meters
                    distance(md.mask.groundedice_levelset<0)=NaN;
                    pos = find(distance>0);
                    basin_mask = zeros(md.mesh.numberofvertices,1);
                    basin_mask(pos)=1;
                    Thickness_basin = Thickness.*basin_mask;
                    Bed_basin = Bed.*basin_mask;

                    meanbed=mean(Bed_basin(md.mesh.elements2d),2);
                    meanthickness=mean(Thickness_basin(md.mesh.elements2d),2);

                    H_AF=meanthickness + min(md.materials.rho_water/md.materials.rho_ice.*meanbed,0);
                    H_AF(H_AF<0)=0;
                    areas=GetAreas(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d);
                    vol_AF_basin = areas.*H_AF;
                    vol_all_2K(kk,jj) = sum(vol_AF_basin);
                end
            end
        end
        clear meanbed; clear meanthickness; clear md
        for ii = 3
            md=loadmodel(org,model{ii});
            for jj = 1:yr
                for kk = 1:length(id)
                    Bed = md.results.TransientSolution(jj).Bed;
                    Thickness =  md.results.TransientSolution(jj).Thickness;

                    basin = ['DrainageBasins/' id{kk} '.exp'];
                    distance=ExpToLevelSet(md.mesh.x2d,md.mesh.y2d,basin); %in meters
                    distance(md.mask.groundedice_levelset<0)=NaN;
                    pos = find(distance>0);
                    basin_mask = zeros(md.mesh.numberofvertices,1);
                    basin_mask(pos)=1;
                    Thickness_basin = Thickness.*basin_mask;
                    Bed_basin = Bed.*basin_mask;

                    meanbed=mean(Bed_basin(md.mesh.elements2d),2);
                    meanthickness=mean(Thickness_basin(md.mesh.elements2d),2);

                    H_AF=meanthickness + min(md.materials.rho_water/md.materials.rho_ice.*meanbed,0);
                    H_AF(H_AF<0)=0;
                    areas=GetAreas(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d);
                    vol_AF_basin = areas.*H_AF;
                    vol_all_5K(kk,jj) = sum(vol_AF_basin);
                end
            end
        end
        
        clear meanbed; clear meanthickness; clear md
        for ii = 4
            md=loadmodel(org,model{ii});
            for jj = 1:yr
                for kk = 1:length(id)
                    Bed = md.results.TransientSolution(jj).Bed;
                    Thickness =  md.results.TransientSolution(jj).Thickness;

                    basin = ['DrainageBasins/' id{kk} '.exp'];
                    distance=ExpToLevelSet(md.mesh.x2d,md.mesh.y2d,basin); %in meters
                    distance(md.mask.groundedice_levelset<0)=NaN;
                    pos = find(distance>0);
                    basin_mask = zeros(md.mesh.numberofvertices,1);
                    basin_mask(pos)=1;
                    Thickness_basin = Thickness.*basin_mask;
                    Bed_basin = Bed.*basin_mask;

                    meanbed=mean(Bed_basin(md.mesh.elements2d),2);
                    meanthickness=mean(Thickness_basin(md.mesh.elements2d),2);

                    H_AF=meanthickness + min(md.materials.rho_water/md.materials.rho_ice.*meanbed,0);
                    H_AF(H_AF<0)=0;
                    areas=GetAreas(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d);
                    vol_AF_basin = areas.*H_AF;
                    vol_all_8K(kk,jj) = sum(vol_AF_basin);
                end
            end
        end
        
        clear meanbed; clear meanthickness; clear md
        for ii = 5
            md=loadmodel(org,model{ii});
            for jj = 1:yr
                for kk = 1:length(id)
                    Bed = md.results.TransientSolution(jj).Bed;
                    Thickness =  md.results.TransientSolution(jj).Thickness;

                    basin = ['DrainageBasins/' id{kk} '.exp'];
                    distance=ExpToLevelSet(md.mesh.x2d,md.mesh.y2d,basin); %in meters
                    distance(md.mask.groundedice_levelset<0)=NaN;
                    pos = find(distance>0);
                    basin_mask = zeros(md.mesh.numberofvertices,1);
                    basin_mask(pos)=1;
                    Thickness_basin = Thickness.*basin_mask;
                    Bed_basin = Bed.*basin_mask;

                    meanbed=mean(Bed_basin(md.mesh.elements2d),2);
                    meanthickness=mean(Thickness_basin(md.mesh.elements2d),2);

                    H_AF=meanthickness + min(md.materials.rho_water/md.materials.rho_ice.*meanbed,0);
                    H_AF(H_AF<0)=0;
                    areas=GetAreas(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d);
                    vol_AF_basin = areas.*H_AF;
                    vol_all_12K(kk,jj) = sum(vol_AF_basin);
                end
            end
        end
        
        
        mass_all_ctrl = 917.*vol_all_ctrl./(10^12); %total mass Gt
        mass_all_2K = 917.*vol_all_2K./(10^12); %total mass Gt
        mass_all_5K = 917.*vol_all_5K./(10^12); %total mass Gt
        mass_all_8K = 917.*vol_all_8K./(10^12); %total mass Gt
        mass_all_12K = 917.*vol_all_12K./(10^12); %total mass Gt

        c = colormap(brewermap(10,'Reds'));
        for b=1:18;

            figure(b)
            hold on;
            plot(time,mass_all_2K(b,:)-mass_all_ctrl(b,:),'color',[c(4,:)],'LineWidth',2);
            plot(time,mass_all_5K(b,:)-mass_all_ctrl(b,:),'color',[c(6,:)],'LineWidth',2);
            plot(time,mass_all_8K(b,:)-mass_all_ctrl(b,:),'color',[c(8,:)],'LineWidth',2);
            plot(time,mass_all_12K(b,:)-mass_all_ctrl(b,:),'color',[c(10,:)],'LineWidth',2);
            title(['Region ' id{b} ' Mass Loss']);
            legend('Exp 1','Exp2','Exp 3','Exp 4')
            xlabel('Years')
            ylabel('Change in Mass (Gt)')
            ylim([-2300 20])
            grid on
            set(gca,'FontSize',16,'LineWidth',1)
            hold off
        end

end



% Step 6: Plots of Thawable Area Per Basin {{{
if perform(org, 'ThawableArea');
    md=loadmodel(org,'Initialization');
    T = md.results.ThermalSolution.Temperature;
    Tm = T + md.materials.beta.*md.initialization.pressure; % pressure adjusted temperature
    
    load('Masks/BedVelMask_100_km_Of_100_m_per_yr.mat')

    Thw = 273.14;
    for Frz = 261.15:1:273.15
    
        id ={'D-Dp','C-Cp','Ap-B','Cp-D','A-Ap','B-C','K-A','Jpp-K','E-Ep','Dp-E',... % E. Ant (10)
        'G-H','H-Hp','Ep-F','F-G','J-Jpp',... % W. Ant (5)
        'Hp-I','I-Ipp','Ipp-J'}; % Ant. Peninsula (3)

        for ii = 1:length(id)

            basin = ['DrainageBasins/' id{ii} '.exp'];

            distance=ExpToLevelSet(md.mesh.x2d,md.mesh.y2d,basin); %in meters
            distance(md.mask.groundedice_levelset<0)=NaN;
            pos = find(distance>0);
            basin_mask = NaN(md.mesh.numberofvertices,1);
            basin_mask(pos)=1;

            TmMasked = BedVelMask.*basin_mask.*Tm; % Keep only coastal regions
            TmMasked(TmMasked==0)= NaN;
            Tb = project2d(md,TmMasked,1);
            indicesRegion = find(~isnan(Tb));
            RegionMask = zeros(length(md.mesh.x2d),1);
            RegionMask(indicesRegion)=1;
            Region = RegionMask(md.mesh.elements2d);
            areas = GetAreas(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d);
            RegionAreas = areas.*Region;
            totRegionArea(ii,Frz-260.15) = sum(mean(RegionAreas,2));

            indicesFtoT = find(Tb<Thw & Tb>Frz);
            ThawableMask = zeros(length(md.mesh.x2d),1);
            ThawableMask(indicesFtoT)=1;
            Thawable = ThawableMask(md.mesh.elements2d);
            ThawableAreas = areas.*Thawable;
            totThawableArea(ii,Frz-260.15) = sum(mean(ThawableAreas,2));
        end
    
    end
        
    figure(1) %Just West Ant.
    c = colormap(brewermap(7,'*PuRd'));
    p = plot([0:1:12],totThawableArea(11:15,:)'./totRegionArea(11:15,:)'*100,'LineWidth',2);
    set(p, {'color'}, num2cell(c(1:end-2,:),2));grid on
    set(gca,'XDir','reverse','FontSize',16,'LineWidth',2)
    xlabel('Thwable Range (degrees below PMP)')
    ylabel('% of basin area thawable')
    ylim([0 60])
    title('West Ant. Basins Thawable Extents')
    legend([id(11:15)],'Location','EastOutside')

    figure(2) % Just East Ant.
    c = colormap(brewermap(12,'*YlOrBr'));
    p = plot([0:1:12],totThawableArea(1:10,:)'./totRegionArea(1:10,:)'*100,'LineWidth',2);
    set(p, {'color'}, num2cell(c(1:end-2,:),2));grid on
    set(gca,'XDir','reverse','FontSize',16,'LineWidth',2)
    xlabel('Thwable Range (degrees below PMP)')
    ylabel('% of basin area thawable')
    ylim([0 60])
    title('East Ant. Basins Thawable Extents')
    legend([id(1:10)],'Location','EastOutside')
    
end

