clear all;
close all;
%load rapidmoc obs. data
load moc_transports.mat
load H04S5_amoc.mat

% rapid moc data is twice daily data starting on April 2nd 2004 
% 60=data for one month
mon=fix(length(moc_mar_hc10)/60);
moc=nan(1,mon);
for i=1:mon
moc(i)=mean(moc_mar_hc10((i-1)*60+1:i*60));
end

mon=56*12+3; %year*mon/per+april
%%
% raw simu. and obs. amoc 
figure(1)
plot(amoc(1:end))
hold on;
plot(mon+1:mon+length(moc),moc,'r');


% zoomed in short period
figure(2)
plot(amoc(mon+1:mon+length(moc)),'b');
% plot(amoc(mon+1:end),'b');
hold on
plot(moc,'r');

%% using twice daily data
%zoomed in short period, 12 moth filtered, full length
figure(3);
[xhi,amoc2]=myfilter(amoc,12,3);
h1=plot(amoc2,'b');
hold on;
[xhi,moc2]=myfilter(moc_mar_hc10,720,3);
plot(mon+1:1/60:mon+1+(length(moc2)-1)/60,moc2,'r');

%zoomed in short period, 12 moth filtered
figure(4);
[xhi,amoc2]=myfilter(amoc,3,3);
plot(amoc2(mon+1:mon+length(moc)),'g');
hold on;
[xhi,moc2]=myfilter(moc_mar_hc10,180,3);
plot(0:1/60:(length(moc_mar_hc10)-1)/60,moc2,'r');

%% using monthly rapidmoc data
%zoomed in short period, 12 moth filtered, full length
figure(5);
[xhi,amoc2]=myfilter(amoc,12,3);
h1=plot(amoc2,'b');
hold on;
[xhi,moc2]=myfilter(moc,12,3);
plot(mon+1:mon+length(moc),moc,'r');

%zoomed in short period, 12 moth filtered
figure(6);
[xhi,amoc2]=myfilter(amoc,12,3);
plot(amoc2(mon+1:mon+length(moc)),'g');
hold on;
[xhi,moc2]=myfilter(moc,12,3);
% plot(0:1/60:(length(moc_mar_hc10)-1)/60,moc2,'r');
plot(moc2,'r');
%% using monthly rapidmoc data
%zoomed in short period, 12 moth filtered, full length
hf7=figure(7);
[xhi,amoc2]=myfilter(amoc,12,3);
h1=plot(amoc2,'b');
hold on;
[xhi,moc2]=myfilter(moc,12,3);
plot(mon+1:mon+length(moc),moc,'r');

%amoc: ref. cummingham, grl, 2007,bryden,science,2005
% 1957.10, 22.9 Sv, corresponds data point: 12*(1957-1948)+10=118
% 1981.08, 18.7 Sv, 12*(1981-1948)+8=404
% 1992.07, 19.4 Sv, 12*(1992-1948)+7=535
% 1998.02, 16.1 Sv, 12*(1998-1948)+2=602
% 2004.04, 14.8 Sv, 12*(2004-1948)+4=676
plot([118,404,535,602,676],[22.9,18.7,19.4,16.1,14.8],'.',...
    'markersize',15,'color','g');

%zoomed in short period, 3 moth filtered
hf8=figure(8);
[xhi,amoc2]=myfilter(amoc,3,3);
plot(amoc2(mon+1:mon+length(moc)),'g');
hold on;
[xhi,moc2]=myfilter(moc,3,3);
plot(moc2,'r');

corr(moc2',amoc2(mon+1:mon+length(moc))')

[h_m h_i]=inset(hf7,hf8);

line([600,100],[28,27],'LineStyle','--')
set(h_m,'ylim',[8 35]);
set(h_i,'xlim',[1,length(moc2)]);


%% using monthly rapidmoc data

clear all;
close all;
%load rapidmoc obs. data
load ~/output/H04S5/amoc/moc_transports.mat
% load ~/output/H04S5/amoc/H04S5_amoc.mat
load ~/output/NOIAF_SPINUP2/amoc/NOIAF_SPINUP2_amoc.mat
amoc=amoc(5040:5040+719)+amoc(5760:5760+719)+amoc(6480:6480+719)+amoc(7200:7200+719);
amoc=amoc/4;

amoc=reshape(amoc,1,length(amoc));
% rapid moc data is twice daily data starting on April 2nd 2004 
% 60=data for one month
mon=fix(length(moc_mar_hc10)/60);
moc=nan(1,mon);
for i=1:mon
moc(i)=mean(moc_mar_hc10((i-1)*60+1:i*60));
end

mon=56*12+3; %year*mon/per+april

% using monthly rapidmoc data
%zoomed in short period, 12 moth filtered, full length
figure(9);
subplot(2,1,1);
[xhi,amoc2]=myfilter(amoc,12,3);
amoc2=amoc;
h1=plot(amoc2,'b');
hold on;
[xhi,moc2]=myfilter(moc,12,3);
plot(mon+1:mon+length(moc),moc,'r-');
% plot(mon+1:mon+length(moc),moc,'-','color',[.5 .5 .5]);
% set(gca,'xlim',[0,length(amoc)+1],'xtick',1:120:length(amoc),...
%     'xticklabel',1948:10:1948+fix(length(amoc)/12),...
%     'ylim',[0.9*min([moc,amoc2(mon+1:mon+length(moc))]),...
%     1.1*max([moc,amoc2(mon+1:mon+length(moc))])],...
%     'xminortick','on');
axis tight;
set(gca,'xtick',1:120:length(amoc),...
    'xticklabel',1948:10:1948+fix(length(amoc)/12),...
    'ylim',[0.9*min([moc,amoc2(mon+1:end)]),...
    1.1*max([moc,amoc2(mon+1:end)])],...
    'xminortick','on');
set(gca,'fontsize',12,'fontweight','bold','linewidth',1.0);
ylabel('Sv','fontsize',12,'fontweight','bold');

yr=fix(length(amoc)/12);
amoc3=nan(1,yr);
for i=1:yr
amoc3(i)=mean(amoc((i-1)*12+1:i*12));
end
plot(6:12:12*yr,amoc3,'b','linewidth',1.5); %set june 1948, as the postion of annual mean.

%amoc: ref. cummingham, grl, 2007,bryden,science,2005
% 1957.10, 22.9 Sv, corresponds data point: 12*(1957-1948)+10=118
% 1981.08, 18.7 Sv, 12*(1981-1948)+8=404
% 1992.07, 19.4 Sv, 12*(1992-1948)+7=535
% 1998.02, 16.1 Sv, 12*(1998-1948)+2=602
% 2004.04, 14.8 Sv, 12*(2004-1948)+4=676
plot([118,404,535,602,676],[22.9,18.7,19.4,16.1,14.8],'.',...
    'markersize',15,'color','g');

%zoomed in short period, 3 moth filtered
subplot(2,1,2);
[xhi,amoc2]=myfilter(amoc,3,3);
amoc2=(amoc-mean(amoc))/std(amoc);
% amoc2=amoc;
% plot(amoc2(mon+1:mon+length(moc)),'b','linewidth',1.5);
plot(amoc2(mon+1:end),'b','linewidth',1.5);

% plot(amoc2(mon+1:mon+length(moc)),'k');

hold on;
[xhi,moc2]=myfilter(moc,3,3);

len=length(moc);
if len>length(amoc2(mon+1:end))
    len=length(amoc2(mon+1:end));
end
% moc2=moc(1:len);
moc2=moc(1:end);
moc2=(moc2-mean(moc2))/std(moc2);
% moc2=moc;

plot(moc2,'-','color','r','linewidth',1.5);
% plot(moc2,'-','color',[0.5 0.5 0.5]);
% set(gca,'xlim',[0,length(moc2)+1],'xtick',10:12:length(moc2),...
%     'xticklabel',2005:1:2004+fix(length(moc2)/12),...
%     'ylim',[0.9*min([moc2,amoc2(mon+1:mon+length(moc2))]),...
%     1.1*max([moc2,amoc2(mon+1:mon+length(moc2))])],...
%     'xminortick','on');
% the obs. data start on April, 2004, however, here take first data as Jan. 2004
% only for a nice axis
axis tight;
set(gca,'xtick',1:12:length(moc2),...
    'xticklabel',2004:1:2004+fix(length(moc2)/12),...
    'ylim',[0.9*min([moc2,amoc2(mon+1:mon+len)]),...
    1.1*max([moc2,amoc2(mon+1:mon+len)])],...
    'xminortick','on');
set(gca,'fontsize',12,'fontweight','bold','linewidth',1.0);
xlabel('Year','fontsize',12,'fontweight','bold');
ylabel('Sv','fontsize',12,'fontweight','bold');
text(2.5,59,'(b)','fontsize',14,'fontweight','bold');
text(2.5,27,'(c)','fontsize',14,'fontweight','bold');


r=corr(moc2(1:len)',amoc2(mon+1:mon+len)')

mean(amoc2(mon+1:mon+len))
text(50,1.8,['r=',num2str(r, '%3.2f')],'fontsize',14,'fontweight','bold');

savefig('NOIAF_SPINUP2_moc_rapidmoc', 'png', 'eps', '-rgb', '-lossless', '-r200')