clear all
geodir = '/data/home/yxzhang/IGSM/regrid/';
ar =ncread(strcat(geodir,'grid_igsm.nc'),'rA');    % grid area
wet=ncread(strcat(geodir,'grid_igsm.nc'),'HFacC'); % fraction of water in grid (vertically)

delz=[10, 15, 20, 30, 40, ...
    55, 70, 90, 110, 135, ...
    160, 190, 225, 265, 310, ...
    360, 415, 475, 540, 610, ...
    685, 765];

tvolume=zeros(144,90,22);
for level=1:22
    tvolume(:,:,level)=ar*delz(level);
end
wvolume=tvolume.*wet;
 
depth={'10m','25m','45m','75m','115m', ...
    '170m','240m','330m','440m','575m', ...
    '735m','925m','1150m','1415m','1725m', ...
    '2085m','2500m','2975m','3515m','4125m', ...
    '4810m','5575m'};
density=[950,1048,1146,950,1048,1146,950,1048,1146, ...
    950,1048,1146,950,1048,1146,950,1048,1146, ...
    900,1048,1196,900,1048,1196,900,1048,1196, ...
    900,1048,1196,900,1048,1196,900,1048,1196, ...
    1410,1410,1410,1410,1410,1410, ...
    1050,1050,1050,1050,1050,1050, ...
    1050,1050,1050,1050,1050,1050, ...
    550,550,550,550,550,550]; % kg m-3

Dir_Prefix = '/data/home/ympeng/MITgcm-Plastic/verification/global_hg_igsm/run_Covid19_optimum_microp_macrop/';
year=linspace(349920,520560,80);
year1=linspace(2021,2100,80);
% PTRtave (kg m-3)
pdata_sf=zeros(80,1);
pdata_oc=zeros(80,1);
% % -------------- data reading and processing -------------- % %
for yr=1:80
    dataf_sf=zeros(144,90,22,60);
    for n=1:60
        data_sf=rdmds(strcat(Dir_Prefix,'PTRtave',num2str(n,'%02d')),year(yr));
        data_sf(data_sf<0.0)=0.0;
        dataf_sf(:,:,:,n)=data_sf.*wvolume;
    end
    
    dataf2_sf=zeros(144,90,22);
    for n=1:60
        if(density(n)<1048)
            dataf2_sf=dataf2_sf+dataf_sf(:,:,:,n);
        else
            dataf2_sf=dataf2_sf+dataf_sf(:,:,:,n)*0.05;
        end
    end
    data_oc=dataf2_sf./ar*1000000;
    
    data2_sf=data_oc(:,:,1);
    pdata_sf(yr,1)=sum(data2_sf(:));
    
    data2_oc=squeeze(sum(data_oc,3));
    pdata_oc(yr,1)=sum(data2_oc(:));
end
% Beach (kg m-3)
pdata_bc=zeros(80,1);
% % -------------- data reading and processing -------------- % %
for yr=1:80
    dataf_bc=zeros(144,90,60);
    for n=1:60
        data_bc=rdmds(strcat(Dir_Prefix,'diags/Beach',num2str(n,'%02d')),year(yr));
        data_bc(data_bc<0.0)=0.0;
        dataf_bc(:,:,n)=data_bc.*delz(1);   % kg km-2
    end
    data2_bc=dataf_bc*1000000;
    data3_bc=squeeze(sum(data2_bc,3));
    pdata_bc(yr,1)=sum(data3_bc(:));
end
% Bottom  (kg m-2)
pdata_bt=zeros(80,1);
% -------------- data reading and processing -------------- % %
for yr=1:80
    dataf_bt=zeros(144,90,60);
    for n=1:60
        data_bt=rdmds(strcat(Dir_Prefix,'diags/Bottom',num2str(n,'%02d')),year(yr));
        data_bt(data_bt<0.0)=0.0;
        dataf_bt(:,:,n)=data_bt;
    end
    data2_bt=dataf_bt*1000000;
    data3_bt=squeeze(sum(data2_bt,3));
    pdata_bt(yr,1)=sum(data3_bt(:));
end

% save fig4_ts.mat

%%
load fig4_ts.mat
p1=figure(1);
sum=pdata_oc+pdata_bc+pdata_bt;
perc_sf=(pdata_sf./sum)*100;
perc_oc=(pdata_oc./sum)*100;
perc_bc=(pdata_bc./sum)*100;
perc_bt=(pdata_bt./sum)*100;
perc=[perc_sf perc_oc perc_bc perc_bt];

subplot(2,1,1)
x=1:80;
[AX,H1,H2] = plotyy(x,pdata_sf,x,perc_sf,'plot');
set(gca,'position',[0.15 0.55 0.75 0.4]);
set(H1,'LineStyle','-','LineWidth',1.2)
set(H2,'LineStyle','--','LineWidth',1.2)
set(AX,'xlim',[-1 82]);
set(gca,'xtick',[]);

set(get(AX(1),'Ylabel'),'String','Concentration (kg km^{-2})',...
    'FontWeight','Bold','FontName',...
    'Liberation Serif','FontSize',13)
set(AX(1),'ylim',[0 10])
set(AX(1),'ytick',[0 5 10],...
    'FontWeight','Bold','FontName',...
    'Liberation Serif','FontSize',11)

set(get(AX(2),'Ylabel'),'String','%',...
    'FontWeight','Bold','FontName',...
    'Liberation Serif','FontSize',13)
set(AX(2),'ylim',[0 4])
set(AX(2),'ytick',[0 2 4],...
    'FontWeight','Bold','FontName',...
    'Liberation Serif','FontSize',11)

hh=legend('surface plastic concentration','% surface plastic');
set(hh,'box','off','FontName',...
    'Liberation Serif','FontSize',13)

subplot(2,1,2)
perc=[perc_bt,perc_bc,perc_oc];
hold on
c_matrix=[160,82,45
          210,180,140
          135,206,235]./255;   % color of top
color_matrix=zeros(243,3);
A=3:3:243;
B=2:3:242;
C=1:3:241;
for a=1:243
    if ismember(a,A)==1
        color_matrix(a,:)=c_matrix(3,:);
    elseif ismember(a,B)==1
        color_matrix(a,:)=c_matrix(2,:);
    elseif ismember(a,C)==1
        color_matrix(a,:)=c_matrix(1,:);
    end
end

for i=1:80
    b=bar(i:i+2,[perc(i,:,:);0,0,0;0,0,0],1,'stacked');
    set(b(1),'facecolor',color_matrix((i-1)*3+1,:))
    set(b(2),'facecolor',color_matrix((i-1)*3+2,:))
    set(b(3),'facecolor',color_matrix((i-1)*3+3,:))
end
set(gca,'position',[0.01 0.1 0.8 0.4]);
xlim([-1,82]);
ylim([0,121]);
set(gca,'position',[0.15 0.12 0.75 0.4],...
    'Xtick',[1 20:20:80],'Xticklabel',...
    {'2021' '2040' '2060' '2080' '2100'},...
    'ytick',[0 25 50 75 100],...
    'FontName','Liberation Serif','FontSize',11);
xlabel('Year','FontWeight','Bold','FontName',...
    'Liberation Serif','FontSize',13);
ylabel('%','FontWeight','Bold','FontName',...
    'Liberation Serif','FontSize',13);

h=legend('Sediment','Beach','Seawater');
set(h,'position',[0.3 0.46 0.5 0.05],'box','off',...
    'orientation','horizon')

print(gcf,'-dpng','-r300',...
    strcat('/data/home/ympeng/MITgcm-Plastic/draw/Fig4.png'));