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