%   Powered by He Yanchun, PKU,NZC
%   $Revision: 1.0  $  $Date: 2010/04/25 19:13:41 $
%   $Revision: 2.0  $  $Date: 2011/12/04 19:13:41 $

clear all;
close all;

%voltr(time, type, layer, section) ;
% section_name =
	%"Drake passage                 ",
	%"South Africa - Antarctica     ",
	%"Australia - Antarctica        ",
	%"Atlantic Ocean 24N            ",
	%"Atlantic Ocean 48N            " ;
%layer =
%	1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
%	19, 20, 21, 22,23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35 ;

%clear all;
%close all;
%% Define Constants
ExpID='H04S5';
%YEAR='2005.01';
%define the isopycnal layers
theta=[29.514,30.119,30.724,31.329,31.934,...
       32.539,33.144,33.749,34.354,34.947,...
       35.507,35.954,36.285,36.530,36.712,...
       36.846,36.946,37.020,37.074,37.115,...
       37.145,37.175,37.205,37.235,37.265,...
       37.295,37.325,37.355,37.385,37.415,...
       37.440,37.460,37.480,37.580,37.800];
% keep two significance digits,and transfer it into string
% theta=round(theta*100)/100;
%sigma=zeros(size(theta,2),5);
%convert theta from num to str for using in 'xticklabel'
kdm=length(theta);
for i=1:kdm
	sigma(i,:)=num2str(theta(i),'%3.3f');
end

% label={'(a)','(b)','(c)','(d)','(e)','(f)','(g)','(h)'};
label={'a)','b)','c)','d)','e)','f)','g)','h)'};
%% Define Vars
% unix(['ncra -v voltr ',ExpID,'_voltr_ave.*.nc ',ExpID,'_voltr_ave_1948-2007.nc']);
%% Load Data
unix(['ls ~/Micom/output/',ExpID,'/voltr_24.5N/',ExpID,'_voltr_ave_1948-2007.nc>filelist']);
files=importdata('filelist'); %including H04S5_voltr_ave_1948-2007.nc
unix('rm -rf filelist');
lenfile=length(files);
sec_name=nc_varget([ExpID,'_voltr.2000.05.nc'],'section_name');
for j=1:lenfile
voltr=nc_varget(files{j},'voltr');
voltr=squeeze(voltr(1,:,:));         %only keep the net transport
sec=size(sec_name,1);                %get the number of sections
fig=1;                                          %id of subplot

secid=[20];
% for secid=[2 3] %secid: secid=2,24N; secid=3, 48N
for kk=[20] %secid: section secidd : secid=2,24N; i=3, 48N

	maxim=fix(max(max(abs(voltr(2:36,kk)))));  %2:36, layer 2-36, note the first is the sum
    if kk==2||kk==3
        maxim=fix(max(max(abs(voltr(2:36,2:3)))));
    end
%     ylim=[-20 20];
    if maxim>10
        ytick=(-fix(maxim/5):fix(maxim/5))*5;
        ylim=[-maxim-2 maxim+2];
%         textypos=-maxim; %postion of the 'text string'
        textypos=0.9*ylim(1);
    elseif maxim<=10 && maxim>5
        ytick=(-fix(maxim/2):fix(maxim/2))*2;
        ylim=[-maxim-1.2 maxim+1.2];
%         textypos=-maxim-1.0;
        textypos=0.9*ylim(1);
    elseif maxim<=5 && maxim>2
        ytick=(-maxim-1:maxim+1);
        ylim=[-maxim-1.1 maxim+1.1];
        textypos=0.9*ylim(1);
    else
        maxim=max(max(abs(voltr(2:36,kk))));
        ytick=(-ceil(maxim):0.5:ceil(maxim));
        ylim=[-ceil(maxim) ceil(maxim)];
        textypos=0.9*ylim(1);
    end
%     subplot(3,2,[fig,fig+2]); 
%     hf1=subplot(3,1,[1,2]); 
    hf1=subplot(4,1,[1,2]); 

    %set the subplot figure's size
    pos=get(gca,'position');
    pos(2)=pos(2)-0.02;
    if fig==1
        pos(1)=pos(1)+0.04;
    else
        pos(1)=pos(1)-0.04;
    end
%     set(gca,'position',pos);
    

    %H_ax1: the layer axis
    H_ax1=gca;
    kdm=22;
	set(H_ax1,'xlim',[1 kdm],'ylim',ylim,'xaxislocation','bottom',...
		'box','on','ygrid','on','fontsize',12,'xtick',[1:kdm],...
		'xticklabel',{'ML',2:kdm},'ytick',ytick,'yticklabel',{ytick},...
        'view',[90 90],'nextplot','add','fontweight','bold');
%     set(H_ax1,'ytick',[]); %delete the bottom axis' tick
    if fig==2
%         set(H_ax1,'xtick',[]);
    end
    if fig==1
        xlabel('Model layer','fontsize',12);
    end

	H_ax2=axes('position',get(H_ax1,'position'));
	set(H_ax2,'xlim',[1 kdm],'ylim',ylim,'xaxislocation','top',...
		'box','on','xgrid','on','ygrid','on','fontsize',12,'xtick',[1:kdm],...
		'xticklabel',{['<',num2str(0.5*(theta(1)+theta(2)),'%2.2f')],...
        sigma(2:kdm-1,:),['>',num2str(0.5*(theta(kdm-1)+theta(kdm)),...
        '%2.2f')]},'ytick',ytick,'view',[90 90],'nextplot','add',...
        'fontweight','bold');
%     set(H_ax2,'xtick',[]);%delete the right tick
     ylabel(H_ax2,'Volume transport (Sv)','fontsize',12,'fontweight','bold');
     
	H_ax3=axes('position',get(H_ax2,'position'));
	set(H_ax3,'xlim',[1 kdm],'ylim',ylim,'xaxislocation','top',...
		'box','on','xgrid','on','ygrid','on','fontsize',12,'xtick',[1:kdm],...
		'xticklabel',{['<',num2str(0.5*(theta(1)+theta(2)),'%2.2f')],...
        sigma(2:kdm-1,:),['>',num2str(0.5*(theta(kdm-1)+theta(kdm)),...
        '%2.2f')]},'ytick',ytick,'view',[90 90],'nextplot','add','fontweight','bold');
     set(H_ax3,'yticklabel',[],'tickdir','in','linewidth',1,'gridlinestyle',':'); %delete the botom tick
    if fig==1
        set(H_ax3,'xticklabel',[]);
    end
     %draw a tick
     for tt=ytick
         plot([34,35],[tt tt],'k-');
     end
     %draw a rectangular to enclose the deep flow
%      if i==2
%       plot(H_ax3,[16,21],[1,1],'-k','linewidth',3);
%       plot(H_ax3,[16,16],[0.3,1],'-k','linewidth',3');
%       plot(H_ax3,[21,21],[0.3,1],'-k','linewidth',3);
%       
%       plot(H_ax3,[8,10],[1,1],'-k','linewidth',3);
%       plot(H_ax3,[8,8],[0.3,1],'-k','linewidth',3');
%       plot(H_ax3,[10,10],[0.3,1],'-k','linewidth',3);
%      else
%       plot(H_ax3,[16,21],[1,1],'-k','linewidth',1);
%       plot(H_ax3,[16,16],[0.3,1],'-k','linewidth',1');
%       plot(H_ax3,[21,21],[0.3,1],'-k','linewidth',1);
%      end
% 	xlabel(H_ax2,'layer density (\sigma_2)','fontsize',18,'fontweight','bold');
% 	ylabel(H_ax2,'Transport (Sv, positive eastward/northward)',...
%         'fontsize',18,'fontweight','bold');
% 	H_title=title(H_ax1,strcat('Density-Partitioned Net Transport across: ',...
%         sec_name(i,:)));
% 	H_title=title(H_ax1,sec_name(i,:));
% 	set(H_title,'fontsize',18,'fontweight','bold');
	H_bar=bar(H_ax3,[1:kdm],voltr(2:kdm+1,kk));
% 	set(H_bar,'Barwidth',1,'FaceColor',[0.1 0.1 0.1]);
%     set(H_bar,'Barwidth',0.8,'FaceColor',[80 182 200]/255,...
%         'edgecolor',[0.2 0.2 0.2]);
    set(H_bar,'Barwidth',0.8,'FaceColor',[0.1 0.1 0.1],...
        'edgecolor',[0.5 0.5 0.5]);
    
%accumlated vertical profile
 %     for kk=1:35;
%         acc(kk)=sum(voltr(2:kk+1,i),1);
%     end
%     plot(acc);
      text(23,7.5,'$$\sigma_2$$','fontsize',15,'fontweight','bold',...
          'interpreter','latex');  

%     text(34,textypos,['total=',num2str(voltr(1,i),'%2.3f'),'Sv'],'fontsize',18,...
%         'fontweight','bold')
    text(20.5,-7,['net = ',num2str(voltr(1,kk),'%2.2f'),' Sv'],'fontsize',11,'fontweight','bold')
%     text(2.5,textypos,label{fig},'fontsize',14,'fontweight','bold')
    text(0.4,-7.2,label{fig},'fontsize',14,'fontweight','bold')
	set(gcf,'Position',[0 400 1024 768],'PaperPositionMode','auto','Renderer','painters');
   
    if fig==2
        xlabel('\sigma_2 Potential density','fontsize',12);
    end
%	print -painter -dpng -r300 strcat(EXPID,'_voltr_ave.','YEAR','png');
%	saveas(gcf, strcat(EXPID,'_voltr_ave_sec',int2str(i),'.',int2str(YEAR)),'png');
% 	saveas(gcf, int2str(i),'png');
% 	oldname=strcat(int2str(i),'.png');
% 	newname=strcat(EXPID,'_voltr_sec',int2str(i),'.',YEAR,'nomxl','.png');
% 	movefile(oldname,newname);
     grid off %close it when print in eps
      set(H_ax3,'xgrid','on');
% 	close;
    fig=fig+1;
end
%     filename=strcat(EXPID,'_voltr.',files{j}(17:20));
% %    print('-dpng','-r300',[filename,'.png']);
% %     print('-deps','-r300',[filename,'.eps']);
% % %     close;
% % print -painters-dpng -r300 filename.png
% % print -painters -deps -r300 filename.eps
end
%%

% clear all;
% close all;
%load rapidmoc obs. data
load ~/output/tracer/H04S5/amoc/moc_transports.mat

% %load simuated amoc between 20N-50N
% load ~/output/H04S5/amoc/H04S5_amoc.mat

% % load transport across 24N ?total southward transport whole column,layer
% layer 8-10 and 16-21
% load ~/output/H04S5/voltr/voltr_southward_times_layer8-10_16-21.mat

% % load transport across 24N ?total southward transport whole column,layer
% layer 8-10 and 16-21
% still keep 24.5N
load ~/output/tracer/H04S5/voltr_24.5N/voltr_southward_times_all_southward_layers.mat

amoc=-vol_24N;

% 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
hf2=subplot(4,1,3);
% [xhi,amoc2]=myfilter(amoc,12,3);
amoc2=amoc;
% h1=plot(amoc2,'color',[.5 .5 .5]);
h1=plot(amoc2,'color','k','linewidth',0.5);
hold on;
% [xhi,moc2]=myfilter(moc,12,3);
%plot Rapid MOC obs.
% 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([min(moc),min(amoc2(mon+1:mon+length(moc)))]),...
%     1.1*max([max(moc),max(amoc2(mon+1:mon+length(moc)))])],...
%     'xminortick','on');
set(gca,'xlim',[0,684+1],'xtick',[1:120:671,672],...
    'xticklabel',[1948:10:2003,2004],...
    'ylim',[0.9*min([min(moc),min(amoc2(mon+1:mon+length(moc)))]),...
    1.1*max([max(moc),max(amoc2(mon+1:mon+length(moc)))])],...
    '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,'k','linewidth',1.5); %set june 1948, as the postion of annual mean.
plot(6:12:12*yr,amoc3,'k','linewidth',2.5); 

%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','k');
text(0.1,31,'b)','fontsize',14,'fontweight','bold');

disp('min voltr(1:672), 1948-2004');
min(amoc2(1:672))
disp('max voltr(1:672), 1948-2004');
max(amoc2(1:672))
%% use 26.5N
load ~/output/tracer/H04S5/voltr_26.5N/voltr_southward_times_all_southward_layers.mat
amoc=-vol_24N;
%zoomed in short period, 3 moth filtered
hf3=subplot(4,1,4);
% [xhi,amoc2]=myfilter(amoc,3,3);
% amoc2=(amoc-mean(amoc))/std(amoc);
amoc2=amoc;
plot(amoc2(mon+1:mon+length(moc)),'k','linewidth',1.5);
% plot(amoc2(mon+1:mon+length(moc)),'k');

hold on;
% [xhi,moc2]=myfilter(moc,3,3);
% moc2=(moc-mean(moc))/std(moc);
moc2=moc;
%plot RAPID MOC obs.
% plot(moc2,'-','color','r','linewidth',1.5);
plot(moc2,'k:','linewidth',2);

% 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
set(gca,'xlim',[0,length(moc2)+1],'xtick',1:12:length(moc2),...
    'xticklabel',2004:1:2004+fix(length(moc2)/12),...
    'ylim',[0.9*min([min(moc2),min(amoc2(mon+1:mon+length(moc2)))]),...
    1.1*max([max(moc2),max(amoc2(mon+1:mon+length(moc2)))])],...
    'xminortick','on');
set(gca,'fontsize',12,'fontweight','bold','linewidth',1.0);
xlabel('Year','fontsize',12,'fontweight','bold');
ylabel('Sv','fontsize',12,'fontweight','bold');
% text(0.1,27,'(c)','fontsize',14,'fontweight','bold');
text(0.1,31,'c)','fontsize',14,'fontweight','bold');

% r=corr(moc2',amoc2(mon+1:mon+length(moc))')
[r,P,RLO,RUP]=corrcoef(reshape(moc2,[length(moc2),1]),...
    reshape(amoc2(mon+1:mon+length(moc)),[length(mon+1:mon+length(moc)),1]));

text(55,27,['r=',num2str(r(2), '%3.2f')],'fontsize',14,'fontweight','bold');

disp('mean voltr, 2004-2009');
mean(amoc2(mon+1:mon+length(moc)))

disp('min voltr, 2004-2009');
min(amoc2(mon+1:mon+length(moc)))
disp('max voltr, 2004-2009');
max(amoc2(mon+1:mon+length(moc)))


disp('min obs, 2004-2009');
min(moc2)
disp('max obs, 2004-2009');
max(moc2)

%%
south=sum(squeeze(voltr([9:11,17:22],2)));
percent=voltr(18,2)/south*100;
% percent=41.0947
disp('mean simulated mean southward trans=');
mean(amoc2(mon+1:mon+length(moc)))
disp('correlation=');
r(2)
disp('significance=');
1-P(2)

%%
set(gcf,'color','w');
set(gcf,'Position',[0 400 1024 768],'PaperPositionMode','auto','Renderer','painters');
figsize(8.26,11.69)

% print -painters -dpng -r300 micom_transports_atl24N-50N_rapidmoc.png
% print -painters -depsc -r300 micom_transports_atl24N-50N_rapidmoc.eps
% saveas(gcf,'micom_transports_atl24N_rapidmoc_v3.fig'); 
% 
% export_fig (gcf, 'micom_transports_atl265N_rapidmoc_v3',...
%     '-pdf', '-png','-m2','-r300', '-a2', '-q101', '-painters', '-rgb');