clear all
addpath(genpath('/data/home/yxzhang/DarwinFuture/library'));
addpath(genpath('/data/home/yxzhang/NASA'));
addpath(genpath('/data/home/ppwu/draw'));
load("grid.mat");
load("mycmap.mat");
lon_144x90=lon_144x90-180;
landareas=shaperead('ne_10m_land.shp');

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];

Dir_Prefix = '/data/home/ympeng/MITgcm-Plastic/verification/global_hg_igsm/run_Covid19_optimum_microp_macrop/';

yy=[1,5,80];
ww=0.22; ll=0.185;
x1=0.11; x2=0.35; x3=0.59;
y1=0.19; y2=0.365; y3=0.54;
year=linspace(349920,520560,80);
year1=linspace(2021,2100,80);
% % -------------- data reading and processing -------------- % % (kg m-3)
for yr=yy
    %% Surface (kg m-3)
    data_s=zeros(144,90,22,40);
    for n=[1:12 19:30 37:40 43:46 49:52 55:58]
        dataf_s=rdmds(strcat(Dir_Prefix,'PTRtave',num2str(n,'%02d')),year(yr));
        dataf_s(dataf_s<0.0)=0.0;
        data_s(:,:,:,n) = dataf_s.*wvolume;
    end

    data2=zeros(144,90);
    for n=[1:12 19:30 37:40 43:46 49:52 55:58]
        if(density(n)<1048)
            data2(:,:)=data2(:,:)+data_s(:,:,1,n);
        else
            data2(:,:)=data2(:,:)+data_s(:,:,1,n)*0.05;
        end
    end
    
    dataf2=data2./ar*1000000*1000;
    data3=log10(dataf2);
    data3(data3==-Inf)=NaN;
    
    u=zeros(144,90);
    u(:,:)=data3(:,:);
    pdata=zeros(90,144);
    pdata=u(:,:)';
    pdata2(:,1:72)=pdata(:,73:144);
    pdata2(:,73:144)=pdata(:,1:72);
    
    % p1=figure(1);
    set(gcf,'Color',[1 1 1]);
    if yr==1
        subplot(3,3,1)
        pcolor(lon_144x90,lat_144x90,pdata2);
        set(gca,'position',[x1 y3 ww ll])
    elseif yr==5
        subplot(3,3,2)
        pcolor(lon_144x90,lat_144x90,pdata2);
        set(gca,'position',[x2 y3 ww ll])
    elseif yr==80
        subplot(3,3,3)
        pcolor(lon_144x90,lat_144x90,pdata2);
        set(gca,'position',[x3 y3 ww ll])
    end
    shading flat
    colormap(widl)
    axesm('mercator','Frame','off','Grid','off',...
        'origin',[0 180 0]);
    mapshow(landareas,'FaceColor',[0.45 0.45 0.45],...
        'EdgeColor',[0.45 0.45 0.45]);
    grid off
    
    caxis([-4,2]);
    hcb1=colorbar('position',[0.82 0.558 0.01 0.15],...
        'FontWeight','Bold','FontName',...
        'Liberation Serif','FontSize',8, ...
        'Ticks',-4:2:2,'YTickLabel',...
        {'10^{-4}','10^{-2}','10^{0}','10^{2}'});
    set(get(hcb1,'Title'),'string','g km^{-2}',...
        'FontWeight','Bold')
    
    %% Beach (kg m-3)
    data_b=zeros(144,90,40);
    for n=[1:12 19:30 37:40 43:46 49:52 55:58]
        dataf_bc=rdmds(strcat(Dir_Prefix,'diags/Beach',num2str(n,'%02d')),year(yr));   % kg m-3
        dataf_bc(dataf_bc<0.0)=0.0;
        data_b(:,:,n) = dataf_bc.*delz(1);
    end

    data2=squeeze(nansum(data_b*1000000*1000,3));
    data3=log10(data2);
    data3(data3==-Inf)=NaN;
    
    u=zeros(144,90);
    u(:,:)=data3(:,:);
    pdata=zeros(90,144);
    pdata=u(:,:)';
    pdata2(:,1:72)=pdata(:,73:144);
    pdata2(:,73:144)=pdata(:,1:72);
    
    % figure
    set(gcf,'Color',[1 1 1]);
    if yr==1
        subplot(3,2,4)
        pcolor(lon_144x90,lat_144x90,pdata2);
        set(gca,'position',[x1 y2 ww ll])
    elseif yr==5
        subplot(3,3,5)
        pcolor(lon_144x90,lat_144x90,pdata2);
        set(gca,'position',[x2 y2 ww ll])
    elseif yr==80
        subplot(3,3,6)
        pcolor(lon_144x90,lat_144x90,pdata2);
        set(gca,'position',[x3 y2 ww ll])
    end
    shading flat
    colormap(widl)
    axesm('mercator','Frame','off','Grid','off',...
        'origin',[0 180 0]);
    mapshow(landareas,'FaceColor',[0.45 0.45 0.45],...
        'EdgeColor',[0.45 0.45 0.45]);
    grid off
    
    caxis([-2,4]);
    hcb2=colorbar('position',[0.82 0.381 0.01 0.15],...
        'FontWeight','Bold','FontName',...
        'Liberation Serif','FontSize',8, ...
        'Ticks',-2:2:4,'YTickLabel',...
        {'10^{-2}','10^{0}','10^{2}','10^{4}'});
    
    %% Bottom (kg m-2)
    data_bt=zeros(144,90,40);
    for n=[1:12 19:30 37:40 43:46 49:52 55:58]
        dataf_bt=rdmds(strcat(Dir_Prefix,'diags/Bottom',num2str(n,'%02d')),year(yr));   % kg m-2
        dataf_bt(dataf_bt<0.0)=0.0;
        data_bt(:,:,n) = dataf_bt;
    end
    
    data2=squeeze(nansum(data_bt*1000000*1000,3));
    data3 = log10(data2);
    data3(data3==-Inf)=NaN;
    
    u=zeros(144,90);
    u(:,:)=data3(:,:);
    pdata=zeros(90,144);
    pdata=u(:,:)';
    pdata2(:,1:72)=pdata(:,73:144);
    pdata2(:,73:144)=pdata(:,1:72);
    
    % figure
    set(gcf,'Color',[1 1 1]);
    if yr==1
        subplot(3,3,7)
        pcolor(lon_144x90,lat_144x90,pdata2);
        set(gca,'position',[x1 y1 ww ll])
    elseif yr==5
        subplot(3,3,8)
        pcolor(lon_144x90,lat_144x90,pdata2);
        set(gca,'position',[x2 y1 ww ll])
    elseif yr==80
        subplot(3,3,9)
        pcolor(lon_144x90,lat_144x90,pdata2);
        set(gca,'position',[x3 y1 ww ll])
    end
    shading flat
    colormap(widl)
    axesm('mercator','Frame','off','Grid','off',...
        'origin',[0 180 0]);
    mapshow(landareas,'FaceColor',[0.45 0.45 0.45],...
        'EdgeColor',[0.45 0.45 0.45]);
    grid off
    
    caxis([-14,6]);
    hcb3=colorbar('position',[0.82 0.208 0.01 0.15],...
        'FontWeight','Bold','FontName',...
        'Liberation Serif','FontSize',8, ...
        'Ticks',-14:5:6,'YTickLabel',...
        {'10^{-14}','10^{-9}','10^{-4}','10^{1}','10^{6}'});
    
end
hAxis=axes('visible','off');
h1=text(-0.045,0.64,'Surface');
set(h1,'fontsize',10,'FontName','Liberation Serif',...
    'FontWeight','Bold','rotation',90,...
    'HorizontalAlignment','center')
h2=text(-0.045,0.43,'Beach');
set(h2,'fontsize',10,'FontName','Liberation Serif',...
    'FontWeight','Bold','rotation',90,...
    'HorizontalAlignment','center')
h3=text(-0.045,0.21,'Sediment');
set(h3,'fontsize',10,'FontName','Liberation Serif',...
    'FontWeight','Bold','rotation',90,...
    'HorizontalAlignment','center')

t1=text(-0.012,0.58,'a');
set(t1,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');
t2=text(0.295,0.58,'b');
set(t2,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');
t3=text(0.605,0.58,'c');
set(t3,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');
t4=text(-0.012,0.36,'d');
set(t4,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');
t5=text(0.295,0.36,'e');
set(t5,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');
t6=text(0.605,0.36,'f');
set(t6,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');
t7=text(-0.012,0.155,'g');
set(t7,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');
t8=text(0.295,0.155,'h');
set(t8,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');
t9=text(0.605,0.155,'i');
set(t9,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');

t10=text(0.19,0.15,'2021');
set(t10,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');
t11=text(0.495,0.15,'2025');
set(t11,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');
t12=text(0.807,0.15,'2100');
set(t12,'fontsize',9.5,'fontname','Liberation Serif',...
    'FontWeight','Bold');

tite('microplastic','fontsize',12,'fontname','Liberation Serif',...
     'FontWeight','Bold');
print(gcf,'-dpng','-r300',...
    strcat('/data/home/ympeng/MITgcm-Plastic/draw/Fig5.png'));