%%%

close all
clear all

dataPath='UVic_dic_c14';

load(fullfile(dataPath,'data_for_plots'))

semilogy(sphist.nfevalhist,cellfun(@(c)abs(c.Fdicglob),sphist.convdatahist),'o-')
hold on
semilogy(aahist.nfevalhist,cellfun(@(c)abs(c.Fdicglob),aahist.convdatahist),'rx-')
semilogy([0 7500]',[.01 .01]','k-','linewidth',2)
xlim([0 7500])
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('|Net air-sea CO_2 flux| [PgC/y]')  
legend('Direct Integration','Anderson Acceleration','OMIP criterion')
title('Abiotic carbon cycle model - UVic online')
print -dpng uvic_online_dic_omip.png

figure
semilogy(sphist.nfevalhist,cellfun(@(c)abs(c.c14driftfrac),sphist.convdatahist),'o-')
hold on
semilogy(aahist.nfevalhist,cellfun(@(c)abs(c.c14driftfrac),aahist.convdatahist),'rx-')
semilogy([0 7500]',[.98 .98]','k-','linewidth',2)                                                   
xlim([0 7500])
ylim([0 1.5])
grid                
ax=gca;             
ax.YMinorGrid='off';
xlabel('Years')
ylabel(['Fraction of ocean volume with \Delta^{14}C drift<0.001 ' char(8240) ' per year'],'fontsize',8) 
legend('Direct Integration','Anderson Acceleration','OMIP criterion','location','southeast')
title('Abiotic carbon cycle model - UVic online')
print -dpng uvic_online_c14_omip.png

xs=sphist.xhist{end};
xs=reshape(xs,[nb numTracers]);
vs=xs*XTovScaleFac;
xa=aahist.xhist{end};
xa=reshape(xa,[nb numTracers]);
va=xa*XTovScaleFac;

Rocns=vs(:,3:4)./vs(:,1:2);
Rocna=va(:,3:4)./va(:,1:2);
DC14s=(Rocns/rc14std-1)*1e3;
DC14a=(Rocna/rc14std-1)*1e3;

figure
tmps=vs(:,1:2);
tmpa=va(:,1:2);
plot(tmps(:),tmpa(:),'b.')
xlim([1.8 2.2]),ylim([1.8 2.2])
hold on
plot(xlim,xlim,'k-')
grid
xlabel('Direct Integration solution [mol/m^3]')
ylabel('Anderson Acceleration solution [mol/m^3]')
title('DIC - UVic online')
print -dpng uvic_online_dic_comparison.png

figure
tmps=DC14s(:,1:2);
tmpa=DC14a(:,1:2);
plot(tmps(:),tmpa(:),'b.')
xlim([-250 0]),ylim([-250 0])
hold on
plot(xlim,xlim,'k-')
grid
xlabel(['Direct Integration solution [' char(8240) ']'])
ylabel(['Anderson Acceleration solution [' char(8240) ']'])
title(['\Delta^{14}C - UVic online'])
print -dpng uvic_online_c14_comparison.png

%%%

close all
clear all

dataPath='tmm_ocmipabioticcarbon_ecco';

load(fullfile(dataPath,'data_for_plots'))

figure
semilogy(sphist.nfevalhist,cellfun(@(c)abs(c.Fdicglob),sphist.convdatahist),'o-')
hold on
semilogy(aahist.nfevalhist,cellfun(@(c)abs(c.Fdicglob),aahist.convdatahist),'rx-')
semilogy([0 8000]',[.01 .01]','k-','linewidth',2)
xlim([0 8000])
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('|Net air-sea CO_2 flux| [PgC/y]')  
legend('Direct Integration','Anderson Acceleration','OMIP criterion')
title('Abiotic carbon cycle model - ECCO TMM')
print -dpng tmm_ocmipabioticcarbon_ecco_dic_omip.png

figure
semilogy(sphist.nfevalhist,cellfun(@(c)abs(c.c14driftfrac),sphist.convdatahist),'o-')
hold on
semilogy(aahist.nfevalhist,cellfun(@(c)abs(c.c14driftfrac),aahist.convdatahist),'rx-')
semilogy([0 8000]',[.98 .98]','k-','linewidth',2)                                                   
xlim([0 8000])
ylim([0 1.5])
grid                
ax=gca;             
ax.YMinorGrid='off';
xlabel('Years')
ylabel(['Fraction of ocean volume with \Delta^{14}C drift<0.001 ' char(8240) ' per year'],'fontsize',8) 
legend('Direct Integration','Anderson Acceleration','OMIP criterion','location','southeast')
title('Abiotic carbon cycle model - ECCO TMM')
print -dpng tmm_ocmipabioticcarbon_ecco_c14_omip.png

%%%

close all
clear all

dataPath='tmm_ocmipabioticcarbon_eccov4';

load(fullfile(dataPath,'data_for_plots'))

figure
semilogy(sphist.nfevalhist,cellfun(@(c)abs(c.Fdicglob),sphist.convdatahist),'o-')
hold on
semilogy(aahist.nfevalhist,cellfun(@(c)abs(c.Fdicglob),aahist.convdatahist),'rx-')
semilogy([0 8000]',[.01 .01]','k-','linewidth',2)
xlim([0 8000])
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('|Net air-sea CO_2 flux| [PgC/y]')  
legend('Direct Integration','Anderson Acceleration','OMIP criterion')
title('Abiotic carbon cycle model - ECCOv4 TMM')
print -dpng tmm_ocmipabioticcarbon_eccov4_dic_omip.png

figure
semilogy(sphist.nfevalhist,cellfun(@(c)abs(c.c14driftfrac),sphist.convdatahist),'o-')
hold on
semilogy(aahist.nfevalhist,cellfun(@(c)abs(c.c14driftfrac),aahist.convdatahist),'rx-')
semilogy([0 8000]',[.98 .98]','k-','linewidth',2)                                                   
xlim([0 8000])
ylim([0 1.5])
grid                
ax=gca;             
ax.YMinorGrid='off';
xlabel('Years')
ylabel(['Fraction of ocean volume with \Delta^{14}C drift<0.001 ' char(8240) ' per year'],'fontsize',8) 
legend('Direct Integration','Anderson Acceleration','OMIP criterion','location','southeast')
title('Abiotic carbon cycle model - ECCOv4 TMM')
print -dpng tmm_ocmipabioticcarbon_eccov4_c14_omip.png

xs=sphist.xhist{end};
xs=reshape(xs,[nb numTracers]);
vs=xs*XTovScaleFac;
xa=aahist.xhist{end};
xa=reshape(xa,[nb numTracers]);
va=xa*XTovScaleFac;

Rocns=vs(:,2)./vs(:,1);
Rocna=va(:,2)./va(:,1);
DC14s=(Rocns-1)*1e3;
DC14a=(Rocna-1)*1e3;

figure
plot(vs(:,1),va(:,1),'b.')
xlim([1.4 2.3]),ylim([1.4 2.3])
hold on
plot(xlim,xlim,'k-')
grid
xlabel('Direct Integration solution [mol/m^3]')
ylabel('Anderson Acceleration solution [mol/m^3]')
title('DIC - ECCOv4 TMM')
print -dpng tmm_ocmipabioticcarbon_eccov4_dic_comparison.png

figure
plot(DC14s,DC14a,'b.')
xlim([-450 0]),ylim([-450 0])
hold on
plot(xlim,xlim,'k-')
grid
xlabel(['Direct Integration solution [' char(8240) ']'])
ylabel(['Anderson Acceleration solution [' char(8240) ']'])
title(['\Delta^{14}C - ECCOv4 TMM'])
print -dpng tmm_ocmipabioticcarbon_eccov4_c14_comparison.png

%%%

close all
clear all

dataPath='tmm_idealage_ecco';

load(fullfile(dataPath,'data_for_plots'))

figure
semilogy(sphist.nfevalhist,cellfun(@(c)c(1),sphist.vnormshist),'o-')                   
hold on
semilogy(aahist.nfevalhist,cellfun(@(c)c(1),aahist.vnormshist),'rx-')
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('||residual||')  
legend('Direct Integration','Anderson Acceleration')
title('Ideal age - ECCO TMM')
print -dpng tmm_idealage_ecco.png

figure
semilogy(sphist.nfevalhist,cellfun(@(c)c.agedriftfrac,sphist.convdatahist),'o-')                   
hold on
semilogy(aahist.nfevalhist,cellfun(@(c)c.agedriftfrac,aahist.convdatahist),'rx-')
semilogy([0 8000]',[.98 .98]','k-','linewidth',2)                                                   
xlim([0 8000])
ylim([0 1.5])
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('Fraction of ocean volume with ideal age drift<0.01 year per year','fontsize',8) 
legend('Direct Integration','Anderson Acceleration','OMIP criterion','location','southeast')
title('Ideal age - ECCO TMM')
print -dpng tmm_idealage_drift_ecco.png

xs=sphist.xhist{end};
xs=reshape(xs,[nb numTracers]);
vs=xs*XTovScaleFac;
xa=aahist.xhist{end};
xa=reshape(xa,[nb numTracers]);
va=xa*XTovScaleFac;

figure
plot(vs,va,'.')
xlim([0 4000]),ylim([0 4000])
hold on
plot(xlim,xlim,'k-')
grid
xlabel('Direct Integration solution [y]')
ylabel('Anderson Acceleration solution [y]')
% gtext(['RMS difference: ' sprintf('%4.3f y',sqrt(mean((va-vs).^2)))],'fontsize',12,'fontweight','bold')
title('Ideal age - ECCO TMM')
print -dpng tmm_idealage_ecco_comparison.png

%%%

close all
clear all

dataPath='tmm_path';

load(fullfile(dataPath,'data_for_plots'))

xs=sphist.xhist{end};
xs=reshape(xs,[nb 2]);

for iJob=1:maxNumTracers
  numTracers=numTracersJobs(iJob);
  vscaleFac=vscaleJobs{iJob};
  XTovScaleFac=diag(vscaleFac);
  vToXScaleFac=diag(1./vscaleFac);
  vs{iJob}=xs(:,iJob)*XTovScaleFac;
  xa=aa{iJob}.aahist.xhist{end};
  xa=reshape(xa,[nb numTracers]);
  va{iJob}=xa*XTovScaleFac;
  if iJob==1
	vs{iJob}=vs{iJob}*fmoltodpmPa; % dpm/m^3
	va{iJob}=va{iJob}*fmoltodpmPa; % dpm/m^3
  elseif iJob==2	
	vs{iJob}=vs{iJob}*fmoltodpmTh; % dpm/m^3
	va{iJob}=va{iJob}*fmoltodpmTh; % dpm/m^3
  end           
end

iJob=1;
figure
semilogy(sphist.nfevalhist,cellfun(@(c)c(iJob),sphist.vnormshist),'o-') 
hold on
semilogy(aa{iJob}.aahist.nfevalhist,cellfun(@(c)c(1),aa{iJob}.aahist.vnormshist),'rx-')
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('||residual||')  
legend('Direct Integration','Anderson Acceleration')
title('^{231}Pa - UVic TMM')
print -dpng tmm_path_uvic_pa231_separate.png

iJob=2;
figure
semilogy(sphist.nfevalhist,cellfun(@(c)c(iJob),sphist.vnormshist),'o-') 
hold on
semilogy(aa{iJob}.aahist.nfevalhist,cellfun(@(c)c(1),aa{iJob}.aahist.vnormshist),'rx-')
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('||residual||')  
legend('Direct Integration','Anderson Acceleration')
title('^{230}Th - UVic TMM')
print -dpng tmm_path_uvic_th230_separate.png

iJob=1;
figure
plot(vs{iJob},va{iJob},'.')
xlim([0 1]),ylim(xlim)
hold on
plot(xlim,xlim,'k-')
grid
xlabel('Direct Integration solution [dpm/m^3]')
ylabel('Anderson Acceleration solution [dpm/m^3]')
title('^{231}Pa - UVic TMM')
print -dpng tmm_path_uvic_pa231_comparison_separate.png

iJob=2;
figure
plot(vs{iJob},va{iJob},'.')
xlim([0 3]),ylim(xlim)
hold on
plot(xlim,xlim,'k-')
grid
xlabel('Direct Integration solution [dpm/m^3]')
ylabel('Anderson Acceleration solution [dpm/m^3]')
title('^{230}Th - UVic TMM')
print -dpng tmm_path_uvic_th230_comparison_separate.png

%%%

close all
clear all

dataPath='tmm_znpo4dop';

load(fullfile(dataPath,'data_for_plots'))

figure
semilogy(sphist.nfevalhist,cellfun(@(c)c(1),sphist.vnormshist),'o-')                   
hold on
semilogy(aa{1}.aahist.nfevalhist,cellfun(@(c)c(1),aa{1}.aahist.vnormshist),'rx-')
semilogy(aa{2}.aahist.nfevalhist,cellfun(@(c)c(1),aa{2}.aahist.vnormshist),'gx-')
semilogy(aa{3}.aahist.nfevalhist,cellfun(@(c)c(1),aa{3}.aahist.vnormshist),'cx-')
semilogy(aa{4}.aahist.nfevalhist,cellfun(@(c)c(1),aa{4}.aahist.vnormshist),'mx-')
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('||residual||')  
legend('Direct Integration','Anderson Acceleration (m=50)','Anderson Acceleration (m=40, restart)','Anderson Acceleration (m=30, restart)','Anderson Acceleration (m=20, restart)')
title('PO_4 - MITgcm2.8deg TMM')
print -dpng tmm_znpo4dop_mitgcm2p8_po4.png

figure
semilogy(sphist.nfevalhist,cellfun(@(c)c(2),sphist.vnormshist),'o-')                   
hold on
semilogy(aa{1}.aahist.nfevalhist,cellfun(@(c)c(2),aa{1}.aahist.vnormshist),'rx-')
semilogy(aa{2}.aahist.nfevalhist,cellfun(@(c)c(2),aa{2}.aahist.vnormshist),'gx-')
semilogy(aa{3}.aahist.nfevalhist,cellfun(@(c)c(2),aa{3}.aahist.vnormshist),'cx-')
semilogy(aa{4}.aahist.nfevalhist,cellfun(@(c)c(2),aa{4}.aahist.vnormshist),'mx-')
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('||residual||')  
legend('Direct Integration','Anderson Acceleration (m=50)','Anderson Acceleration (m=40, restart)','Anderson Acceleration (m=30, restart)','Anderson Acceleration (m=20, restart)')
title('DOP - MITgcm2.8deg TMM')
print -dpng tmm_znpo4dop_mitgcm2p8_dop.png

figure
semilogy(sphist.nfevalhist,cellfun(@(c)c(3),sphist.vnormshist),'o-')                   
hold on
semilogy(aa{1}.aahist.nfevalhist,cellfun(@(c)c(3),aa{1}.aahist.vnormshist),'rx-')
semilogy(aa{2}.aahist.nfevalhist,cellfun(@(c)c(3),aa{2}.aahist.vnormshist),'gx-')
semilogy(aa{3}.aahist.nfevalhist,cellfun(@(c)c(3),aa{3}.aahist.vnormshist),'cx-')
semilogy(aa{4}.aahist.nfevalhist,cellfun(@(c)c(3),aa{4}.aahist.vnormshist),'mx-')
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('||residual||')  
legend('Direct Integration','Anderson Acceleration (m=50)','Anderson Acceleration (m=40, restart)','Anderson Acceleration (m=30, restart)','Anderson Acceleration (m=20, restart)')
title('Zn - MITgcm2.8deg TMM')
print -dpng tmm_znpo4dop_mitgcm2p8_zn.png

vs=cellfun(@(xs)reshape(xs,[nb numTracers])*XTovScaleFac,sphist.xhist,'UniformOutput',0);
va=cellfun(@(xa)reshape(xa,[nb numTracers])*XTovScaleFac,aa{1}.aahist.xhist,'UniformOutput',0);

figure
plot(vs{end}(:,1),va{end}(:,1),'.')
xlim([0 3]),ylim(xlim)
hold on
plot(xlim,xlim,'k-')
grid
xlabel('Direct Integration solution [{\mu}mol/l]')
ylabel('Anderson Acceleration solution [{\mu}mol/l]')
title('PO_4 - MITgcm2.8deg TMM')
print -dpng tmm_znpo4dop_mitgcm2p8_po4_comparison.png

figure
plot(vs{end}(:,2),va{end}(:,2),'.')
xlim([0 .1]),ylim(xlim)
hold on
plot(xlim,xlim,'k-')
grid
xlabel('Direct Integration solution [{\mu}mol/l]')
ylabel('Anderson Acceleration solution [{\mu}mol/l]')
title('DOP - MITgcm2.8deg TMM')
print -dpng tmm_znpo4dop_mitgcm2p8_dop_comparison.png

figure
plot(vs{end}(:,3),va{end}(:,3),'.')
xlim([0 .035]),ylim(xlim)
hold on
plot(xlim,xlim,'k-')
grid
xlabel('Direct Integration solution [{\mu}mol/l]')
ylabel('Anderson Acceleration solution [{\mu}mol/l]')
title('Zn - MITgcm2.8deg TMM')
print -dpng tmm_znpo4dop_mitgcm2p8_zn_comparison.png

%%%

close all
clear all

dataPath='tmm_preformedtracer_uvic';

load(fullfile(dataPath,'data_for_plots'))

for iJob=1:maxNumTracers
  numTracers=numTracersJobs(iJob);
  vscaleFac=vscaleJobs{iJob};
  XTovScaleFac=diag(vscaleFac);
  vToXScaleFac=diag(1./vscaleFac);
  xs=sp{iJob}.sphist.xhist{end};
  xs=reshape(xs,[nb numTracers]);
  vs{iJob}=xs*XTovScaleFac;
  xa=aa{iJob}.aahist.xhist{end};
  xa=reshape(xa,[nb numTracers]);
  va{iJob}=xa*XTovScaleFac;
end

iJob=1;
figure
semilogy(sp{iJob}.sphist.nfevalhist,cellfun(@(c)c(1),sp{iJob}.sphist.vnormshist),'o-') 
hold on
semilogy(aa{iJob}.aahist.nfevalhist,cellfun(@(c)c(1),aa{iJob}.aahist.vnormshist),'rx-')
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('||residual||')  
legend('Direct Integration','Anderson Acceleration')
title('Preformed O_2 - UVic TMM')
print -dpng tmm_preformedtracer_uvic_oxy.png

iJob=2;
figure
semilogy(sp{iJob}.sphist.nfevalhist,cellfun(@(c)c(1),sp{iJob}.sphist.vnormshist),'o-') 
hold on
semilogy(aa{iJob}.aahist.nfevalhist,cellfun(@(c)c(1),aa{iJob}.aahist.vnormshist),'rx-')
grid
ax=gca;
ax.YMinorGrid='off';
xlabel('Years')
ylabel('||residual||')  
legend('Direct Integration','Anderson Acceleration')
title('Preformed PO_4 - UVic TMM')
print -dpng tmm_preformedtracer_uvic_po4.png

iJob=1;
figure
plot(vs{iJob},va{iJob},'.')
xlim([150 450]),ylim(xlim)
hold on
plot(xlim,xlim,'k-')
grid
xlabel('Direct Integration solution [mmol/m^3]')
ylabel('Anderson Acceleration solution [mmol/m^3]')
title('Preformed O_2 - UVic TMM')
print -dpng tmm_preformedtracer_uvic_oxy_comparison.png

iJob=2;
figure
plot(vs{iJob},va{iJob},'.')
xlim([0 3]),ylim(xlim)
hold on
plot(xlim,xlim,'k-')
grid
xlabel('Direct Integration solution [mmol/m^3]')
ylabel('Anderson Acceleration solution [mmol/m^3]')
title('Preformed PO_4 - UVic TMM')
print -dpng tmm_preformedtracer_uvic_po4_comparison.png
