function [ha] = Plot_Acceleration_LocalRates(LALT,TACCBud,M,t,sF,HRrel,SDSL,BarySL,IBAint,GIA_Field_KS,SE,SE_SDSL,SE_Bary,SE_GIA)

% Load GIA Fields

load('C:\KS Trash\1880-2021 EOF Based\GIA_Fields.mat')
load('final_estimates.mat', 'p')
pGIA = sum(reshape(p(end,:),40,100),2);

lon = [-179.5:179.5]';
lat = [-89.5:89.5]';
[lat,lon] = meshgrid(lat,lon);
L1(:,1) = reshape(lon,size(lon,1)*size(lon,2),1);
L1(:,2) = reshape(lat,size(lon,1)*size(lon,2),1);
idx = knnsearch(LALT,L1);

T1 = NaN(size(L1,1),1);
T1 = TACCBud(idx,:)*1000;

T1 = NaN(size(L1,1),1);
T1 = TACCBud(idx,1)*1000;
for i = 1:numel(idx);
    [d1km(i,1) d2km(i,1)]=lldistkm(L1(i,:),LALT(idx(i),:));
end
sdis = find(d1km>100);
T1(sdis,:) = NaN;

ha = tight_subplot(3,2,[.12 .08],[.1 .1],[.1 .1]);

axes(ha(1));
T2 = reshape(T1(:,1),size(lon));
m_proj('robinson','lon',[-180 180]);
m_pcolor(lon,lat,T2);shading flat
m_coast('patch',[1 1 1],'edgecolor','k');
m_grid('XTick',[],'YTick',[]);
caxis([-.25 .25])
load('roma.mat')
colormap(roma);
title('Acceleration since 1970','Fontsize',10)
cb = colorbar('location','Southoutside','YTick',[-.2:.1:.2])
hold on
m_plot([-65 -50 -25 -25 -65 -65],[70 58 58 80 80 70],'Color',[0.49,0.18,0.56],'Linewidth',1.5)
m_plot([-100 -45 -45 -100 -100],[17 17 44 44 17],'Color',[0.49,0.18,0.56],'Linewidth',1.5)
m_plot([-165 -132 -132 -135 -155.8 -165 -165],[-59 -59 -56 -53 -53 -57 -59],'Color',[0.49,0.18,0.56],'Linewidth',1.5)
m_plot([130 170 170 130 130],[28 28 44 44 28],'Color',[0.49,0.18,0.56],'Linewidth',1.5)
m_plot([145 180 180 145 145],[-10 -10 10 10 -10],'Color',[0.49,0.18,0.56],'Linewidth',1.5)

stg1 = inpolygon(LALT(:,1),LALT(:,2),[-65 -50 -25 -25 -65 -65],[70 58 58 80 80 70]);
stg2 = inpolygon(LALT(:,1),LALT(:,2),[-100 -45 -45 -100 -100],[17 17 44 44 17]);
stg3 = inpolygon(LALT(:,1),LALT(:,2),[-165 -132 -132 -135 -155.8 -165 -165],[-59 -59 -56 -53 -53 -57 -59]);
stg4 = inpolygon(LALT(:,1),LALT(:,2),[130 170 170 130 130],[28 28 44 44 28]);
stg5 = inpolygon(LALT(:,1),LALT(:,2),[145 180 180 145 145],[-10 -10 10 10 -10]);


load('final_estimates.mat', 'xs_multi')
xs_multi(:,1:20) = [];
xs_multi(sF,:) = [];
M1 = M;
isn = find(isnan(M1));
Mxs = xs_multi(1:end-1,:)';
M1(isn) = Mxs(isn);

GIA = GIA_Field_KS;

[MHRrelsm(:,1),~,~,~,~,MHRrelsmEns]=ssatrend(mean(HRrel(stg1,:))'*1000,15,mean(SE(stg1,:))'*1000);
[MSDSLsm(:,1),~,~,~,~,MSDSLsmEns]=ssatrend(mean(SDSL(stg1,:))'*1000,15,mean(SE_SDSL(stg1,:))'*1000);
[MBarysm(:,1),~,~,~,~,MBarysmEns]=ssatrend(mean(BarySL(stg1,:))'*1000,15,mean(SE_Bary(stg1,:))'*1000);
yfit = 1*t;yfit = yfit-yfit(end);
[MIBEsm(:,1),~,~,~,~,MIBEsmEns]=ssatrend(mean(IBAint(stg1,:))'*1000+yfit,15,std(mean(IBAint(stg1,:))')*1000);
MIBEsm = MIBEsm-yfit; MIBEsmEns = MIBEsmEns-yfit;

[MHRrelsm(:,2),~,~,~,~,MHRrelsmEns(:,:,2)]=ssatrend(mean(HRrel(stg2,:))'*1000,15,mean(SE(stg2,:))'*1000);
[MSDSLsm(:,2),~,~,~,~,MSDSLsmEns(:,:,2)]=ssatrend(mean(SDSL(stg2,:))'*1000,15,mean(SE_SDSL(stg2,:))'*1000);
[MBarysm(:,2),~,~,~,~,MBarysmEns(:,:,2)]=ssatrend(mean(BarySL(stg2,:))'*1000,15,mean(SE_Bary(stg2,:))'*1000);
yfit = 1*t;yfit = yfit-yfit(end);
[MIBEsm(:,2),~,~,~,~,MIBEsmEns(:,:,2)]=ssatrend(mean(IBAint(stg2,:))'*1000+yfit,15,std(mean(IBAint(stg2,:))')*1000);
MIBEsm(:,2) = MIBEsm(:,2)-yfit; MIBEsmEns(:,:,2) = MIBEsmEns(:,:,2)-yfit;

[MHRrelsm(:,3),~,~,~,~,MHRrelsmEns(:,:,3)]=ssatrend(mean(HRrel(stg3,:))'*1000,15,mean(SE(stg3,:))'*1000);
[MSDSLsm(:,3),~,~,~,~,MSDSLsmEns(:,:,3)]=ssatrend(mean(SDSL(stg3,:))'*1000,15,mean(SE_SDSL(stg3,:))'*1000);
[MBarysm(:,3),~,~,~,~,MBarysmEns(:,:,3)]=ssatrend(mean(BarySL(stg3,:))'*1000,15,mean(SE_Bary(stg3,:))'*1000);
yfit = 1*t;yfit = yfit-yfit(end);
[MIBEsm(:,3),~,~,~,~,MIBEsmEns(:,:,3)]=ssatrend(mean(IBAint(stg3,:))'*1000+yfit,15,std(mean(IBAint(stg3,:))')*1000);
MIBEsm(:,3) = MIBEsm(:,3)-yfit; MIBEsmEns(:,:,3) = MIBEsmEns(:,:,3)-yfit;

[MHRrelsm(:,4),~,~,~,~,MHRrelsmEns(:,:,4)]=ssatrend(mean(HRrel(stg4,:))'*1000,15,mean(SE(stg4,:))'*1000);
[MSDSLsm(:,4),~,~,~,~,MSDSLsmEns(:,:,4)]=ssatrend(mean(SDSL(stg4,:))'*1000,15,mean(SE_SDSL(stg4,:))'*1000);
[MBarysm(:,4),~,~,~,~,MBarysmEns(:,:,4)]=ssatrend(mean(BarySL(stg4,:))'*1000,15,mean(SE_Bary(stg4,:))'*1000);
yfit = 1*t;yfit = yfit-yfit(end);
[MIBEsm(:,4),~,~,~,~,MIBEsmEns(:,:,4)]=ssatrend(mean(IBAint(stg4,:))'*1000+yfit,15,std(mean(IBAint(stg4,:))')*1000);
MIBEsm(:,4) = MIBEsm(:,4)-yfit; MIBEsmEns(:,:,4) = MIBEsmEns(:,:,4)-yfit;

[MHRrelsm(:,5),~,~,~,~,MHRrelsmEns(:,:,5)]=ssatrend(mean(HRrel(stg5,:))'*1000,15,mean(SE(stg5,:))'*1000);
[MSDSLsm(:,5),~,~,~,~,MSDSLsmEns(:,:,5)]=ssatrend(mean(SDSL(stg5,:))'*1000,15,mean(SE_SDSL(stg5,:))'*1000);
[MBarysm(:,5),~,~,~,~,MBarysmEns(:,:,5)]=ssatrend(mean(BarySL(stg5,:))'*1000,15,mean(SE_Bary(stg5,:))'*1000);
yfit = 1*t;yfit = yfit-yfit(end);
[MIBEsm(:,5),~,~,~,~,MIBEsmEns(:,:,5)]=ssatrend(mean(IBAint(stg5,:))'*1000+yfit,15,std(mean(IBAint(stg5,:))')*1000);
MIBEsm(:,5) = MIBEsm(:,5)-yfit; MIBEsmEns(:,:,5) = MIBEsmEns(:,:,5)-yfit;

axes(ha(2))
i = 1;
ciplot(diff(MHRrelsm(:,i))+std(diff(MHRrelsmEns(:,:,i)),0,2)*2,diff(MHRrelsm(:,i))-std(diff(MHRrelsmEns(:,:,i)),0,2)*2,t(2:end),[0 .45 .74])
hold on
ciplot(diff(MSDSLsm(:,i))+std(diff(MSDSLsmEns(:,:,i)),0,2)*2,diff(MSDSLsm(:,i))-std(diff(MSDSLsmEns(:,:,i)),0,2)*2,t(2:end),[.93 .69 .13])
ciplot(diff(MBarysm(:,i))+std(diff(MBarysmEns(:,:,i)),0,2)*2,diff(MBarysm(:,i))-std(diff(MBarysmEns(:,:,i)),0,2)*2,t(2:end),[.85 .33 .1])
RMGIAwm = sum(mean(GIA_Fields_all(stg1,:)).*pGIA',2);
RMGIAwmSE = sqrt(sum((mean(GIA_Fields_all(stg1,:))-RMGIAwm).^2.*pGIA',2));
ciplot((diff(mean(GIA(stg1,:))')+RMGIAwmSE*2)*1000,(diff(mean(GIA(stg1,:))')-+RMGIAwmSE*2)*1000,t(2:end),[0 .5 0])
ciplot(diff(MIBEsm(:,i))+std(diff(MIBEsmEns(:,:,i)),0,2)*2,diff(MIBEsm(:,i))-std(diff(MIBEsmEns(:,:,i)),0,2)*2,t(2:end),[0.70,0.38,0.57])
plot(t(2:end),diff(MHRrelsm(:,i)),'Linewidth',1.5,'Color',[0 .45 .74])
plot(t(2:end),diff(MSDSLsm(:,i)),'Linewidth',1.5,'Color',[.93 .69 .13])
plot(t(2:end),diff(MBarysm(:,i)),'Linewidth',1.5,'Color',[.85 .33 .1])
plot(t(2:end),diff(mean(GIA(stg1,:))')*1000,'Linewidth',1.5,'Color',[0 .5 0])
plot(t(2:end),diff(MIBEsm(:,i)),'Linewidth',1.5,'Color',[0.70,0.38,0.57])

ciplot([11 11],[9 9],[1901 1908],[0 .45 .74])
plot([1901 1908],[10 10],'Color',[0 .45 .74],'Linewidth',1.5)
text(1909,10,'Reconstruction','Fontsize',8)
ciplot([8 8],[6 6],[1901 1908],[.85 .33 .1])
plot([1901 1908],[7 7],'Color',[.85 .33 .1],'Linewidth',1.5)
text(1909,7,'Barystatic GRD','Fontsize',8)
ciplot([11 11],[9 9],[1951 1958],[.93 .69 .13])
plot([1951 1958],[10 10],'Color',[.93 .69 .13],'Linewidth',1.5)
text(1959,10,'Sterodynamic','Fontsize',8)
ciplot([8 8],[6 6],[1951 1958],[0.70,0.38,0.57])
plot([1951 1958],[7 7],'Color',[0.70,0.38,0.57],'Linewidth',1.5)
text(1959,7,'IBE','Fontsize',8)
ciplot([8 8],[6 6],[2001 2008],[0 .5 0])
plot([2001 2008],[7 7],'Color',[0 .5 0],'Linewidth',1.5)
text(2009,7,'GIA','Fontsize',8)

set(gca,'XTick',[1920:20:2000],'XTicklabel',[1920:20:2000],'Xlim',[1899 2023],...
    'YTick',[0:4:8],'YTicklabel',[0:4:8],'Ylim',[-4 12],...
    'Fontsize',8,'Layer','top');
ylabel('Rates [mmyr^{-1}]','FontSize',10)
xlabel('Year','FontSize',10)
title('Greenland','FontSize',10)


axes(ha(3))
i = 3;
ciplot(diff(MHRrelsm(:,i))+std(diff(MHRrelsmEns(:,:,i)),0,2)*2,diff(MHRrelsm(:,i))-std(diff(MHRrelsmEns(:,:,i)),0,2)*2,t(2:end),[0 .45 .74])
hold on
ciplot(diff(MSDSLsm(:,i))+std(diff(MSDSLsmEns(:,:,i)),0,2)*2,diff(MSDSLsm(:,i))-std(diff(MSDSLsmEns(:,:,i)),0,2)*2,t(2:end),[.93 .69 .13])
ciplot(diff(MBarysm(:,i))+std(diff(MBarysmEns(:,:,i)),0,2)*2,diff(MBarysm(:,i))-std(diff(MBarysmEns(:,:,i)),0,2)*2,t(2:end),[.85 .33 .1])
RMGIAwm = sum(mean(GIA_Fields_all(stg3,:)).*pGIA',2);
RMGIAwmSE = sqrt(sum((mean(GIA_Fields_all(stg3,:))-RMGIAwm).^2.*pGIA',2));
ciplot((diff(mean(GIA(stg3,:))')+RMGIAwmSE*2)*1000,(diff(mean(GIA(stg3,:))')-+RMGIAwmSE*2)*1000,t(2:end),[0 .5 0])
ciplot(diff(MIBEsm(:,i))+std(diff(MIBEsmEns(:,:,i)),0,2)*2,diff(MIBEsm(:,i))-std(diff(MIBEsmEns(:,:,i)),0,2)*2,t(2:end),[0.70,0.38,0.57])
plot(t(2:end),diff(MHRrelsm(:,i)),'Linewidth',1.5,'Color',[0 .45 .74])
plot(t(2:end),diff(MSDSLsm(:,i)),'Linewidth',1.5,'Color',[.93 .69 .13])
plot(t(2:end),diff(MBarysm(:,i)),'Linewidth',1.5,'Color',[.85 .33 .1])
plot(t(2:end),diff(mean(GIA(stg3,:))')*1000,'Linewidth',1.5,'Color',[0 .5 0])
plot(t(2:end),diff(MIBEsm(:,i)),'Linewidth',1.5,'Color',[0.70,0.38,0.57])
set(gca,'XTick',[1920:20:2000],'XTicklabel',[1920:20:2000],'Xlim',[1899 2023],...
    'YTick',[-3:3:3],'YTicklabel',[-3:3:3],'Ylim',[-6 6],...
    'Fontsize',8,'Layer','top');
ylabel('Rates [mmyr^{-1}]','FontSize',10)
xlabel('Year','FontSize',10)
title('South Pacific','FontSize',10)

axes(ha(4))
i = 2;
ciplot(diff(MHRrelsm(:,i))+std(diff(MHRrelsmEns(:,:,i)),0,2)*2,diff(MHRrelsm(:,i))-std(diff(MHRrelsmEns(:,:,i)),0,2)*2,t(2:end),[0 .45 .74])
hold on
ciplot(diff(MSDSLsm(:,i))+std(diff(MSDSLsmEns(:,:,i)),0,2)*2,diff(MSDSLsm(:,i))-std(diff(MSDSLsmEns(:,:,i)),0,2)*2,t(2:end),[.93 .69 .13])
ciplot(diff(MBarysm(:,i))+std(diff(MBarysmEns(:,:,i)),0,2)*2,diff(MBarysm(:,i))-std(diff(MBarysmEns(:,:,i)),0,2)*2,t(2:end),[.85 .33 .1])
RMGIAwm = sum(mean(GIA_Fields_all(stg2,:)).*pGIA',2);
RMGIAwmSE = sqrt(sum((mean(GIA_Fields_all(stg2,:))-RMGIAwm).^2.*pGIA',2));
ciplot((diff(mean(GIA(stg2,:))')+RMGIAwmSE*2)*1000,(diff(mean(GIA(stg2,:))')-+RMGIAwmSE*2)*1000,t(2:end),[0 .5 0])
ciplot(diff(MIBEsm(:,i))+std(diff(MIBEsmEns(:,:,i)),0,2)*2,diff(MIBEsm(:,i))-std(diff(MIBEsmEns(:,:,i)),0,2)*2,t(2:end),[0.70,0.38,0.57])
plot(t(2:end),diff(MHRrelsm(:,i)),'Linewidth',1.5,'Color',[0 .45 .74])
plot(t(2:end),diff(MSDSLsm(:,i)),'Linewidth',1.5,'Color',[.93 .69 .13])
plot(t(2:end),diff(MBarysm(:,i)),'Linewidth',1.5,'Color',[.85 .33 .1])
plot(t(2:end),diff(mean(GIA(stg2,:))')*1000,'Linewidth',1.5,'Color',[0 .5 0])
plot(t(2:end),diff(MIBEsm(:,i)),'Linewidth',1.5,'Color',[0.70,0.38,0.57])
set(gca,'XTick',[1920:20:2000],'XTicklabel',[1920:20:2000],'Xlim',[1899 2023],...
    'YTick',[-2:2:6],'YTicklabel',[-2:2:6],'Ylim',[-4 8],...
    'Fontsize',8,'Layer','top');
ylabel('Rates [mmyr^{-1}]','FontSize',10)
xlabel('Year','FontSize',10)
title('Gulf Stream','FontSize',10)

axes(ha(5))
i = 4;
ciplot(diff(MHRrelsm(:,i))+std(diff(MHRrelsmEns(:,:,i)),0,2)*2,diff(MHRrelsm(:,i))-std(diff(MHRrelsmEns(:,:,i)),0,2)*2,t(2:end),[0 .45 .74])
hold on
ciplot(diff(MSDSLsm(:,i))+std(diff(MSDSLsmEns(:,:,i)),0,2)*2,diff(MSDSLsm(:,i))-std(diff(MSDSLsmEns(:,:,i)),0,2)*2,t(2:end),[.93 .69 .13])
ciplot(diff(MBarysm(:,i))+std(diff(MBarysmEns(:,:,i)),0,2)*2,diff(MBarysm(:,i))-std(diff(MBarysmEns(:,:,i)),0,2)*2,t(2:end),[.85 .33 .1])
RMGIAwm = sum(mean(GIA_Fields_all(stg4,:)).*pGIA',2);
RMGIAwmSE = sqrt(sum((mean(GIA_Fields_all(stg4,:))-RMGIAwm).^2.*pGIA',2));
ciplot((diff(mean(GIA(stg4,:))')+RMGIAwmSE*2)*1000,(diff(mean(GIA(stg4,:))')-+RMGIAwmSE*2)*1000,t(2:end),[0 .5 0])
ciplot(diff(MIBEsm(:,i))+std(diff(MIBEsmEns(:,:,i)),0,2)*2,diff(MIBEsm(:,i))-std(diff(MIBEsmEns(:,:,i)),0,2)*2,t(2:end),[0.70,0.38,0.57])
plot(t(2:end),diff(MHRrelsm(:,i)),'Linewidth',1.5,'Color',[0 .45 .74])
plot(t(2:end),diff(MSDSLsm(:,i)),'Linewidth',1.5,'Color',[.93 .69 .13])
plot(t(2:end),diff(MBarysm(:,i)),'Linewidth',1.5,'Color',[.85 .33 .1])
plot(t(2:end),diff(mean(GIA(stg4,:))')*1000,'Linewidth',1.5,'Color',[0 .5 0])
plot(t(2:end),diff(MIBEsm(:,i)),'Linewidth',1.5,'Color',[0.70,0.38,0.57])
set(gca,'XTick',[1920:20:2000],'XTicklabel',[1920:20:2000],'Xlim',[1899 2023],...
    'YTick',[-2:2:6],'YTicklabel',[-2:2:6],'Ylim',[-4 8],...
    'Fontsize',8,'Layer','top');
ylabel('Rates [mmyr^{-1}]','FontSize',10)
xlabel('Year','FontSize',10)
title('Kuroshio Region','FontSize',10)

axes(ha(6))
i = 5;
ciplot(diff(MHRrelsm(:,i))+std(diff(MHRrelsmEns(:,:,i)),0,2)*2,diff(MHRrelsm(:,i))-std(diff(MHRrelsmEns(:,:,i)),0,2)*2,t(2:end),[0 .45 .74])
hold on
ciplot(diff(MSDSLsm(:,i))+std(diff(MSDSLsmEns(:,:,i)),0,2)*2,diff(MSDSLsm(:,i))-std(diff(MSDSLsmEns(:,:,i)),0,2)*2,t(2:end),[.93 .69 .13])
ciplot(diff(MBarysm(:,i))+std(diff(MBarysmEns(:,:,i)),0,2)*2,diff(MBarysm(:,i))-std(diff(MBarysmEns(:,:,i)),0,2)*2,t(2:end),[.85 .33 .1])
RMGIAwm = sum(mean(GIA_Fields_all(stg5,:)).*pGIA',2);
RMGIAwmSE = sqrt(sum((mean(GIA_Fields_all(stg5,:))-RMGIAwm).^2.*pGIA',2));
ciplot((diff(mean(GIA(stg5,:))')+RMGIAwmSE*2)*1000,(diff(mean(GIA(stg5,:))')-+RMGIAwmSE*2)*1000,t(2:end),[0 .5 0])
ciplot(diff(MIBEsm(:,i))+std(diff(MIBEsmEns(:,:,i)),0,2)*2,diff(MIBEsm(:,i))-std(diff(MIBEsmEns(:,:,i)),0,2)*2,t(2:end),[0.70,0.38,0.57])
plot(t(2:end),diff(MHRrelsm(:,i)),'Linewidth',1.5,'Color',[0 .45 .74])
plot(t(2:end),diff(MSDSLsm(:,i)),'Linewidth',1.5,'Color',[.93 .69 .13])
plot(t(2:end),diff(MBarysm(:,i)),'Linewidth',1.5,'Color',[.85 .33 .1])
plot(t(2:end),diff(mean(GIA(stg5,:))')*1000,'Linewidth',1.5,'Color',[0 .5 0])
plot(t(2:end),diff(MIBEsm(:,i)),'Linewidth',1.5,'Color',[0.70,0.38,0.57])
set(gca,'XTick',[1920:20:2000],'XTicklabel',[1920:20:2000],'Xlim',[1899 2023],...
    'YTick',[0:3:6],'YTicklabel',[0:3:6],'Ylim',[-3 9],...
    'Fontsize',8,'Layer','top');
ylabel('Rates [mmyr^{-1}]','FontSize',10)
xlabel('Year','FontSize',10)
title('Tropical Western Pacific','FontSize',10)

end