clear

%grid
lats=-89:89; nz=90;
fn='/nobackup/hzhang1/pub/llc2160/grid/RF.data';
RF=readbin(fn,nz+1);
fn='/nobackup/hzhang1/pub/llc2160/grid/RC.data';
RC=readbin(fn,nz);

load MOC_year.mat
amoc_yr=amoc; gmoc_yr=gmoc; pmoc_yr=pmoc;
amoc_yr(amoc_yr==0)=nan;
gmoc_yr(gmoc_yr==0)=nan;
pmoc_yr(pmoc_yr==0)=nan;

kk=find(lats<-35|lats>70); amoc_yr(kk,:)=nan;
kk=find(lats<-35|lats>65); pmoc_yr(kk,:)=nan;


figure(1)
fig_ps=[752 56 1537 1038];
set(gcf,'Position',fig_ps)

subplot(311)
pcolor(lats,RF,gmoc_yr')
caxis([-1 1]*50)
shading flat, colorbar
ylim([-6000 0])
title('GMOC')

subplot(312)
pcolor(lats,RF,amoc_yr')
caxis([-1 1]*50)
shading flat, colorbar
ylim([-6000 0])
title('AMOC')

subplot(313)
pcolor(lats,RF,pmoc_yr')
caxis([-1 1]*50)
shading flat, colorbar
ylim([-6000 0])
title('PMOC')

colormap(jet)


%monthly
load MOC_month.mat
amoc_yr2=mean(amoc,3); gmoc_yr2=mean(gmoc,3); pmoc_yr2=mean(pmoc,3);
amoc_yr2(amoc_yr2==0)=nan;
gmoc_yr2(gmoc_yr2==0)=nan;
pmoc_yr2(pmoc_yr2==0)=nan;
kk=find(lats<-35|lats>70); amoc_yr2(kk,:)=nan;
kk=find(lats<-35|lats>65); pmoc_yr2(kk,:)=nan;

amoc_yr3=std(amoc,0,3); gmoc_yr3=std(gmoc,0,3); pmoc_yr3=std(pmoc,0,3);
amoc_yr3(amoc_yr3==0)=nan;
gmoc_yr3(gmoc_yr3==0)=nan;
pmoc_yr3(pmoc_yr3==0)=nan;
kk=find(lats<-35|lats>70); amoc_yr3(kk,:)=nan;
kk=find(lats<-35|lats>65); pmoc_yr3(kk,:)=nan;
%compare w/ annual mean
figure(2)
subplot(311)
pcolor(lats,RF,gmoc_yr'-gmoc_yr2')
caxis([-1 1])
shading flat, colorbar
ylim([-6000 0])
title('GMOC')

subplot(312)
pcolor(lats,RF,amoc_yr'-amoc_yr2')
caxis([-1 1])
shading flat, colorbar
ylim([-6000 0])
title('AMOC')

subplot(313)
pcolor(lats,RF,pmoc_yr'-pmoc_yr2')
caxis([-1 1])
shading flat, colorbar
ylim([-6000 0])
title('PMOC')

colormap(bluewhitered)


figure(3)
set(gcf,'Position',fig_ps)

subplot(311)
pcolor(lats,RF,gmoc_yr3')
caxis([0 50])
shading flat, colorbar
ylim([-6000 0])
title('Std GMOC')

subplot(312)
pcolor(lats,RF,amoc_yr3')
caxis([0 10])
shading flat, colorbar
ylim([-6000 0])
title('Std AMOC')

subplot(313)
pcolor(lats,RF,pmoc_yr3')
caxis([0 50])
shading flat, colorbar
ylim([-6000 0])
title('Std PMOC')

colormap(jet)

%1000m time series
figure(4)
set(gcf,'Position',fig_ps)
subplot(211)
kk=find(lats<-35|lats>70); amoc(kk,:)=nan;
kk=find(lats<-35|lats>65); pmoc(kk,:)=nan;
tmp1=abs(-RC-1000); kk=find(tmp1==min(tmp1)); kk=kk(1);
TT=3:14;
gmoc1k=squeeze(gmoc(:,kk,:)); 
plot(TT,[gmoc1k(90+25,:);gmoc1k(90+35,:);gmoc1k(90+45,:);gmoc1k(90+55,:)],'linew',2)
grid
xlim([2 15])
xlabel('month from 2020/01')
title('annual global overturning at \approx 1000m depth (Sv)')
legend('25N','35N','45N','55N')

subplot(212)
amoc1k=squeeze(amoc(:,kk,:)); 
plot(TT,[amoc1k(90+25,:);amoc1k(90+35,:);amoc1k(90+45,:);amoc1k(90+55,:)],'linew',2)
grid
xlim([2 15])
xlabel('month from 2020/01')
title('annual Atlantic overturning at \approx 1000m depth (Sv)')
legend('25N','35N','45N','55N')



figure(1)
print -dpng lk_trsp_moc_mean
figure(3)
print -dpng lk_trsp_moc_std
figure(4)
print -dpng lk_trsp_moc_time