%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           MICOM DIAGNOSTICS (Matlab package)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Contributed by Mats Bentsen 
% Contact: Detelina Ivanova, detelina.ivanova@nersc.no
% 24/03/2015
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plots and saves plot of the Meridional Overturning Circulation (MOC)
% for choice of region: 1 - Atlantic; 2 - Indo-Pacific; 3- Global
% for choice of vertical coordinates: 0 -isopycnal layers; 1 - depth levels
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; close all;

% Grunch & Hexagon (Change if local set-up)
addpath /work-common/shared/bjerknes/diagnostics/Packages/MICOM_DIAG/matlab_tools/mats
addpath /work-common/shared/bjerknes/diagnostics/Packages/MICOM_DIAG/matlab_tools/m_map1.4
addpath /work-common/shared/bjerknes/diagnostics/Packages/MICOM_DIAG/matlab_tools/m_map1.4f
addpath /work-common/shared/bjerknes/diagnostics/Packages/MICOM_DIAG/matlab_tools/matlab_netcdf_5_0
addpath /work-common/shared/bjerknes/diagnostics/Packages/MICOM_DIAG/matlab_tools/eosben07
addpath /work-common/shared/bjerknes/diagnostics/Packages/MICOM_DIAG/matlab_tools/eoslib05
addpath /work-common/shared/bjerknes/diagnostics/Packages/MICOM_DIAG/matlab_tools/evenmat


% Set a structure with experiments' information:
expr(4)=struct('name','NBF1850OC_f19_tn11_01_130ka', ...
               'name_disp','130ka', ...
               'first_year',1900, ...
               'last_year',1999);

expr(3)=struct('name','NBF1850OC_f19_tn11_01_125ka', ...
               'name_disp','125ka', ...
               'first_year',1900, ...
               'last_year',1999);

expr(2)=struct('name','NBF1850OC_f19_tn11_01_120ka', ...
               'name_disp','120ka', ...
               'first_year',1900, ...
               'last_year',1999);

expr(1)=struct('name','NBF1850OC_f19_tn11_01_115ka', ...
               'name_disp','115ka', ...
               'first_year',1900, ...
               'last_year',1999);

expr(5)=struct('name','NBF1850OC_f19_tn11_02', ...
               'name_disp','PI', ...
               'first_year',1400, ...
               'last_year',1499);



% Set the path to the processed results from this script. You need to create the directory in advance
workpath=['/fimm/work/cgu025/noresm/micom_diag/Eemian/'];

% Set the path to the plots directory. You need to create the directory in advance
picpath=['/fimm/work/cgu025/noresm/micom_diag/Eemian/'];

% Set the depth coodinates: 0 - isopycnal layers; 1 - depth levels
depth_coord = 1;

% Set the region: 1 - Atlantic; 2- Indo-Pacific; 3 - Global Ocean;
region = 1;


mmm=zeros(166,70,3);

% Reading the time averaged model output
for nexp=1:5
  fyear=expr(nexp).first_year;
  lyear=expr(nexp).last_year;
  exp=expr(nexp).name;

  fname = [workpath '/moc/plot_moc_' exp '_hy.'  num2str(fyear) '-' num2str(lyear) 'y.nc' ];

       ncid=netcdf.open(fname,'NC_NOWRITE');

       varid=netcdf.inqVarID(ncid,'lat');
       lat=netcdf.getVar(ncid,varid);
       varid=netcdf.inqVarID(ncid,'depth');
       depth=netcdf.getVar(ncid,varid);
       varid=netcdf.inqVarID(ncid,'region');
       region_name=netcdf.getVar(ncid,varid);

       varid=netcdf.inqVarID(ncid,'mmflxd');
       tmp=netcdf.getVar(ncid,varid);
       fill_value=netcdf.getAtt(ncid,varid,'_FillValue');
       mmflxd=zeros(size(tmp));
       mmflxd_num=zeros(size(tmp));
       ind=find(tmp~=fill_value);
       mmflxd(ind)=tmp(ind);

       varid=netcdf.inqVarID(ncid,'mmflxl');
       tmp=netcdf.getVar(ncid,varid);
       fill_value=netcdf.getAtt(ncid,varid,'_FillValue');
       mmflxl=zeros(size(tmp));
       mmflxl_num=zeros(size(tmp));
       ind=find(tmp~=fill_value);
       mmflxl(ind)=tmp(ind);

       netcdf.close(ncid)



if     region==1
%  vDr = [-4:2:34];
  vDr = [-6:2:28];
  vDr1 = [-3:0.5:3];
  title_name='Atlantic';
  fname_prefix='atlantic';
elseif region==2
  vDr = [-36:4:36];
  title_name='Pacific-Indian';
  fname_prefix='pacific_indian';
elseif region==3
  vDr = [-36:4:36];
  title_name='Global';
  fname_prefix='global';
end
%colormap('jet')
cmap1=cbrewer('div','RdYlBu',17);
colormap(cmap1)



if depth_coord
  f=squeeze(mmflxd(:,:,region))*1e-9;
  [nx ny]=size(f);
  y=-depth;
  ylabel_text='Depth (m)';
  vctag='d';
else
  f=squeeze(mmflxl(:,:,region))*1e-9;
  [nx ny]=size(f);
  y=-(1:ny)';
  ylabel_text='Layer index';
  vctag='l';
end

mask=zeros(size(f));
f(f==0)=NaN;
ind=find(~isnan(f));
mask(ind)=1;
mask1=sum(mask,2);
i1=find(mask1>0,1,'first');
i2=find(mask1>0,1,'last');

mask(:,end)=1;
mask(end,:)=1;
ind=find(mask==0);
lat2=lat*ones(1,ny);
y2=ones(nx,1)*y';

maskx(1,:)=lat2(ind);masky(1,:)=y2(ind);
maskx(2,:)=lat2(ind+1);masky(2,:)=y2(ind+1);
maskx(3,:)=lat2(ind+nx+1);masky(3,:)=y2(ind+nx+1);
maskx(4,:)=lat2(ind+nx);masky(4,:)=y2(ind+nx);

%maskx=[];masky=[];  % modified by Chuncheng

f=extrap(f);
f=max(min(vDr),f);

mmm(:,:,nexp)=f;

end


% -------------------------------------------------------------
tt={'(a)','(b)','(c)','(d)','(e)','(f)'};
%
ff=figure;
set(ff,'paperunits','inches','paperposition',[1,1,7,5],'position',[100 100 1200 800]);
set(0, 'DefaultAxesFontName', 'AvantGarde');hold on;

fon=4;
% -- (a) --
subplot(2,2,1)

[c h]=contourf(lat,y,squeeze(mmm(:,:,5))',vDr,'linewidth',0.3);
h=patch(maskx,masky,'k');
set(h,'edgecolor','k')
set(gca,'xlim',[lat(i1) lat(i2)],'fontsize',2)

hold on;
[c h]=contour(lat,y,squeeze(mmm(:,:,5))',[0 0],'linewidth',1.0,'color','k');hold on;
hold on;pp=plot(gca,[40 40],[-6750 0],'k','linewidth',1.3);
pp=plot(57.44,-3442,'ro','markerfacecolor',[1 0 1]);text(57.44-9,-3000,'MD03-2664','color','w','fontsize',fon);set(pp,'color',[1 0 1]);
pp=plot(33.69,-4584,'ro','markerfacecolor','r');text(33.69-12,-4000,'ODP Site 1063','color','w','fontsize',fon);
pp=plot(30.5,-3480,'rd','markerfacecolor','r');text(30.5-12,-2900,'ODP Site 1060','color','w','fontsize',fon);
pp=plot(5,-4056,'gd','markerfacecolor','g');text(5-7,-3600,'EW9209','color','w','fontsize',fon);
pp=plot(-25.5,-4194,'m^','markerfacecolor',[127 0 255]/255,'markeredgecolor',[127 0 255]/255);text(-25.5-7,-3700,'RC13-229','color','w','fontsize',fon);

ylabel(ylabel_text,'FontSize',6);
xlabel('Latitude','FontSize',6);

cc=ecolorbar(vDr,'r','',1/20);set(cc,'position',[0.46 0.6 0.01 0.3],'fontsize',4);
title(cc,'(Sv)','fontsize',8)

text(55,-6000,'(a)','fontsize',8,'color','w');

set(gca,'LineWidth',1,'FontSize',6) 
title('AMOC (PI)', 'FontSize',6,'FontWeight','Bold') 

% -- (b) --
subplot(2,2,2)
mmmv=squeeze(mmm(121,:,[1 3 5]));
ppp=plot(mmmv,y,'linewidth',1.3);%axis ij;
ylim([-6500 0]);

set(ppp(1),'color','c');
set(ppp(2),'color',[255 102 102]/255);
set(ppp(3),'color','k');

legend(ppp,{'115 ka' '125 ka' 'PI'},'fontsize',6,'location','south');
text(22,-6000,'(b)','fontsize',8,'color','k');
xlabel('Sv','FontSize',6);
title('AMOC (40^{\circ} N)','FontSize',6);
ylabel('Depth (m)','FontSize',6);
set(gca,'LineWidth',1,'FontSize',6)


% -- (c) --
subplot(2,2,3)
dd=squeeze(mmm(:,:,1))'-squeeze(mmm(:,:,5))';
[c h]=contourf(lat,y,dd,vDr1,'linewidth',0.3);
h=patch(maskx,masky,'k');
set(h,'edgecolor','k')
set(gca,'xlim',[lat(i1) lat(i2)],'fontsize',2)

hold on;gg=0.5;
[c h]=contour(lat,y,dd,[0 0],'linewidth',1.0,'color','k');hold on;
hold on;pp=plot(gca,[40 40],[-6750 0],'c','linewidth',1.3);
pp=plot(57.44,-3442,'ro','markerfacecolor',[1 0 1]);text(57.44-9,-3000,'MD03-2664','color',gg*[1 1 1],'fontsize',fon);set(pp,'color',[1 0 1]);
pp=plot(33.69,-4584,'ro','markerfacecolor','r');text(33.69-12,-4000,'ODP Site 1063','color',gg*[1 1 1],'fontsize',fon);
pp=plot(30.5,-3480,'rd','markerfacecolor','r');text(30.5-12,-2900,'ODP Site 1060','color',gg*[1 1 1],'fontsize',fon);
pp=plot(5,-4056,'gd','markerfacecolor','g');text(5-7,-3600,'EW9209','color',gg*[1 1 1],'fontsize',fon);
pp=plot(-25.5,-4194,'m^','markerfacecolor',[127 0 255]/255,'markeredgecolor',[127 0 255]/255);text(-25.5-7,-3700,'RC13-229','color',gg*[1 1 1],'fontsize',fon);
ylabel(ylabel_text,'FontSize',6);
xlabel('Latitude','FontSize',6);

cc=ecolorbar(vDr1,'r','',1/20);set(cc,'position',[0.9 0.1 0.01 0.4],'fontsize',4);
title(cc,'(Sv)','fontsize',8)

text(55,-6000,'(c)','fontsize',8,'color','w');

set(gca,'LineWidth',1,'FontSize',6) 
title('AMOC (115 ka - PI)', 'FontSize',6,'FontWeight','Bold') 

% -- (d) --
subplot(2,2,4)
dd=squeeze(mmm(:,:,3))'-squeeze(mmm(:,:,5))';
[c h]=contourf(lat,y,dd,vDr1,'linewidth',0.3);
h=patch(maskx,masky,'k');
set(h,'edgecolor','k')
set(gca,'xlim',[lat(i1) lat(i2)],'fontsize',2)

hold on;
[c h]=contour(lat,y,dd,[0 0],'linewidth',1.0,'color','k');hold on;
hold on;pp=plot(gca,[40 40],[-6750 0],'color',[255 102 102]/255,'linewidth',1.3);
pp=plot(57.44,-3442,'ro','markerfacecolor',[1 0 1]);text(57.44-9,-3000,'MD03-2664','color',gg*[1 1 1],'fontsize',fon);set(pp,'color',[1 0 1]);
pp=plot(33.69,-4584,'ro','markerfacecolor','r');text(33.69-12,-4000,'ODP Site 1063','color',gg*[1 1 1],'fontsize',fon);
pp=plot(30.5,-3480,'rd','markerfacecolor','r');text(30.5-12,-2900,'ODP Site 1060','color',gg*[1 1 1],'fontsize',fon);
pp=plot(5,-4056,'gd','markerfacecolor','g');text(5-7,-3600,'EW9209','color',gg*[1 1 1],'fontsize',fon);
pp=plot(-25.5,-4194,'m^','markerfacecolor',[127 0 255]/255,'markeredgecolor',[127 0 255]/255);text(-25.5-7,-3700,'RC13-229','color',gg*[1 1 1],'fontsize',fon);
ylabel(ylabel_text,'FontSize',6);
xlabel('Latitude','FontSize',6);

cc=ecolorbar(vDr1,'r','',1/20);set(cc,'position',[0.46 0.1 0.01 0.3],'fontsize',4);
title(cc,'(Sv)','fontsize',8)
colormap(flipud(cmap1));

text(55,-6000,'(d)','fontsize',8,'color','w');

set(gca,'LineWidth',1,'FontSize',6) 
title('AMOC (125ka - PI)', 'FontSize',6,'FontWeight','Bold') 


  

% Save plot
print (gcf,'-dtiff', '-r400', [picpath '/moc/zz_plot_moc_paper_yiming/fig2_version2.tiff'])
