clear all;close all;

TRC='CFC12';
% SAT='999';
% SAT='100';
% SAT='98';
ocnname1={'Arctic Ocean','Nordic Seas','North Atlantic','Tropical Atlantic',...
    'South Atlantic','North Indian','South India','North Pacific','Tropical Pacific','South Pacific','Southern Ocean'};
ocnname2={'AO','NS','NA','TA','SA','NI','SI','NP','TP','SP','SO'};
%%
load grid_tnx1v1.mat parea
parea0=parea';

%fpath='~/archive';
fpath='~/archive2';
dic01=ncread([fpath,'/NOIIANOC_mr03_ctr2/ocn/hist/NOIIANOC_mr03_ctr2.micom.hbgcy.2002.nc'],'dissiclvl');
dic02=ncread([fpath,'/NOIIANOC_mr03_hist2/ocn/hist/NOIIANOC_mr03_hist2.micom.hbgcy.2002.nc'],'dissiclvl');
cantmodel=(dic02-dic01)/1025*1e6;   %mol/m3->umol/kg
depth=ncread('~/tools/data/tnx1v1_depth.nc','depth');
kz=length(depth);

 %idlage=ncread('~/archive/NOIIANOC_mr03_hist2/ocn/hist/NOIIANOC_mr03_hist2.micom.hy.2002.nc','idlagelvl');
idlage=ncread('~/output/carbon/NOIIANOC_mr03_hist2/mage/mage_by_idlage_lvl.nc','mu');

% load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_CFC12_r10_s100_t2002_depth.mat
load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_CFC12_r10_s999_t2002_depth_nordic.mat
CantA=Cant_ttd;
mageA=mage;
% load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_CFC12_r10_s999_t2002_depth_nordic.mat
load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_SF6_r10_s86_t2002_depth_nordic.mat
CantB=Cant_ttd;
mageB=mage;

cantttd=CantA-CantB;
mage1=mageA-mageB;
mage1(isinf(mage1))=NaN;

% load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_CFC12_r8_s100_t2002_depth.mat
load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_CFC12_r8_s999_t2002_depth_nordic.mat
CantA=Cant_ttd;
mageA=mage;
% load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_CFC12_r8_s999_t2002_depth_nordic.mat
load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_SF6_r8_s86_t2002_depth_nordic.mat
CantB=Cant_ttd;
mageB=mage;
cantttd2=CantA-CantB;
mage2=mageA-mageB;
mage2(isinf(mage2))=NaN;

% load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_CFC12_r12_s100_t2002_depth.mat
load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_CFC12_r12_s999_t2002_depth_nordic.mat
CantA=Cant_ttd;
mageA=mage;
% load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_CFC12_r12_s999_t2002_depth_nordic.mat
load ~/output/carbon/NOIIANOC_mr03_hist2/Cant/NOIIANOC_mr03_hist2_Cant_SF6_r12_s86_t2002_depth_nordic.mat
CantB=Cant_ttd;
mageB=mage;
cantttd3=CantA-CantB;
mage3=mageA-mageB;
mage3(isinf(mage3))=NaN;

pmask=ncread('~/tools/data/mask_tnx1v1_basins.nc','pmask');
pmask=double(pmask);
pmask=repmat(pmask,[1 1 kz]);
%for n=2:12
for n=3		% only nordic seas

% unweighted
    %cantmodelp=cantmodel;
    %cantmodelp(find(pmask~=n))=NaN;
    %cantmodelp=squeeze(nanmean(nanmean(cantmodelp)));
%
    %cantttdp=cantttd;
    %cantttdp(find(pmask~=n))=NaN;
    %cantttdp=squeeze(nanmean(nanmean(cantttdp)));
%
    %cantttdp2=cantttd2;
    %cantttdp2(find(pmask~=n))=NaN;
    %cantttdp2=squeeze(nanmean(nanmean(cantttdp2)));
%
    %cantttdp3=cantttd3;
    %cantttdp3(find(pmask~=n))=NaN;
    %cantttdp3=squeeze(nanmean(nanmean(cantttdp3)));
   % 
    %magep1=mage1;
    %magep1(find(pmask~=n))=NaN;
    %magep1=squeeze(nanmean(nanmean(magep1)));
   % 
    %magep2=mage2;
    %magep2(find(pmask~=n))=NaN;
    %magep2=squeeze(nanmean(nanmean(magep2)));
%
    %magep3=mage3;
    %magep3(find(pmask~=n))=NaN;
    %magep3=squeeze(nanmean(nanmean(magep3)));
   % 
    %idlagep=idlage;
    %idlagep(find(pmask~=n))=NaN;
    %idlagep=squeeze(nanmean(nanmean(idlagep)));

% weighted
       cantmodel(find(pmask~=n))=NaN;
       for zn=1:kz
	   parea=parea0;
           cantmodelzn=squeeze(cantmodel(:,:,zn));
           parea(isnan(cantmodelzn))=NaN;
           cantmodelp(zn)=squeeze(nansum(nansum(cantmodelzn.*parea)))/squeeze(nansum(nansum(parea)));
       end

       cantttd(find(pmask~=n))=NaN;
       for zn=1:kz
	   parea=parea0;
           cantttdzn=squeeze(cantttd(:,:,zn));
           parea(isnan(cantttdzn))=NaN;
           cantttdp(zn)=squeeze(nansum(nansum(cantttdzn.*parea)))/squeeze(nansum(nansum(parea)));
       end

       cantttd2(find(pmask~=n))=NaN;
       for zn=1:kz
	   parea=parea0;
           cantttdzn2=squeeze(cantttd2(:,:,zn));
           parea(isnan(cantttdzn2))=NaN;
           cantttdp2(zn)=squeeze(nansum(nansum(cantttdzn2.*parea)))/squeeze(nansum(nansum(parea)));
       end
       
        cantttd3(find(pmask~=n))=NaN;
       for zn=1:kz
	   parea=parea0;
           cantttdzn3=squeeze(cantttd3(:,:,zn));
           parea(isnan(cantttdzn3))=NaN;
           cantttdp3(zn)=squeeze(nansum(nansum(cantttdzn3.*parea)))/squeeze(nansum(nansum(parea)));
       end
       
    mage1(find(pmask~=n))=NaN;
       for zn=1:kz
	   parea=parea0;
           magezn1=squeeze(mage1(:,:,zn));
           parea(isnan(magezn1))=NaN;
           magep1(zn)=squeeze(nansum(nansum(magezn1.*parea)))/squeeze(nansum(nansum(parea)));
       end

    mage2(find(pmask~=n))=NaN;
       for zn=1:kz
	   parea=parea0;
           magezn2=squeeze(mage2(:,:,zn));
           parea(isnan(magezn2))=NaN;
           magep2(zn)=squeeze(nansum(nansum(magezn2.*parea)))/squeeze(nansum(nansum(parea)));
       end
%        
    mage3(find(pmask~=n))=NaN;
       for zn=1:kz
	   parea=parea0;
           magezn3=squeeze(mage3(:,:,zn));
           parea(isnan(magezn3))=NaN;
           magep3(zn)=squeeze(nansum(nansum(magezn3.*parea)))/squeeze(nansum(nansum(parea)));
       end
    
    idlage(find(pmask~=n))=NaN;
       for zn=1:kz
	   parea=parea0;
           idlagezn=squeeze(idlage(:,:,zn));
           parea(isnan(idlagezn))=NaN;
           idlagep(zn)=squeeze(nansum(nansum(idlagezn.*parea)))/squeeze(nansum(nansum(parea)));
       end

   % 
    % plot
    cname='cbpaired';
    mycolors;
    close all;
    figure;
    %subaxis(1,2,1,'Spacing',0.02);
    hold on;
    %plot(cantmodelp,depth,'k-','linewidth',3)
    plot(cantttdp,depth,'color',red,'linewidth',3)
    plot(cantttdp2,depth,'linestyle','--','color',blue,'linewidth',2)
    plot(cantttdp3,depth,'linestyle','--','color',green,'linewidth',2)
    set(gca,'ydir','reverse','ylim',[0 3500],'fontsize',24,'fontweight','bold');
    set(gca,'xlim',[-3.5 0.5]);

    grid on;box on;
%     legend({'TTD','Model'},'location','southeast');
    legend({'TTD,  r=1.0','TTD, r=0.8','TTD, r=1.2'});
    xlabel('Cant (umol/kg)','fontsize',24,'fontweight','bold');
    ylabel('Depth (m)','fontsize',24,'fontweight','bold');
    title(['Cant'],'fontsize',24,'fontweight','bold');
    set(gcf,'color','white');figsize(11.69,8.26);

    Path=['/Users/heyc/Micom/output/carbon/NOIIANOC_mr03_hist2/profile/diff/'];
%     figname=['NOIIANOC_mr03_profile_',ocnname2{n-1},'_',TRC,'_r10_s100-999_T2002_r0812cant'];
    figname=['NOIIANOC_mr03_profile_',ocnname2{n-1},'_CFC12-SF6_r10_s999vs86_T2002_r0812cant'];
    %export_fig (gcf, [Path,figname], '-png','-m2','-r200', '-a2', '-q101', '-painters', '-rgb');
    save2pdf([Path,figname]);

    mfile=[mfilename('fullpath'),'.m'];
    fighist(mfile,Path,figname)
    
    %
    figure;
    %subaxis(1,2,2,'Spacing',0.02);
    hold on;
    %plot(idlagep,depth,'k-','linewidth',3)
    plot(magep1,depth,'color',red,'linewidth',3)
    plot(magep2,depth,'linestyle','--','color',blue,'linewidth',3)
    plot(magep3,depth,'linestyle','--','color',green,'linewidth',3)

    set(gca,'ydir','reverse','yaxislocation','right','ylim',[0 3500],'fontsize',24,'fontweight','bold');
    xlim=get(gca,'xlim');
    %set(gca,'xlim',[-abs(xlim(1)/5), xlim(2)]);
    set(gca,'xlim',[-5,30]);
    grid on;box on;
%     legend({'TTD','Model'},'location','southwest');
    legend({'TTD,  r=1.0','TTD, r=0.8','TTD, r=1.2'});
    xlabel('Age (yr)','fontsize',24,'fontweight','bold');
    ylabel('Depth (m)','fontsize',24,'fontweight','bold');
    title(['Age'],'fontsize',24,'fontweight','bold');
    set(gcf,'color','white');figsize(11.69,8.26);
    % suptitle(['Simulated vs TTD-based for: ',ocnname{n-1}],'fontsize',20,'fontweight','bold');
    %suptitle([ocnname1{n-1}],'fontsize',24,'fontweight','bold');
    
    Path=['/Users/heyc/Micom/output/carbon/NOIIANOC_mr03_hist2/profile/diff/'];
%     figname=['NOIIANOC_mr03_profile_',ocnname2{n-1},'_',TRC,'_r10_s100-999_T2002_r0812mage'];
    figname=['NOIIANOC_mr03_profile_',ocnname2{n-1},'_CFC12-SF6_r10_s999vs86_T2002_r0812mage'];
    save2pdf([Path,figname]);

    mfile=[mfilename('fullpath'),'.m'];
    fighist(mfile,Path,figname)

end
