clc;clear all;close all
addpath('/data3/yiwen/matlab_toolbox/tight_subplot')
addpath('/data3/yiwen/matlab_toolbox/colorbrewer')

%========== Model output ===========%
filename="/data3/yiwen/cesm2/JGR_ocean/Model_output/G_JRA55_test1.pop.h.0031-04.nc"
depth_moc=ncread(filename,'moc_z');
lat_depth=ncread(filename,'lat_aux_grid');

pathi1='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/';
pathi2='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/'
pathi3='/data3/yiwen/cesm2/JGR_ocean/Model_output/';

%------ TSD variables ------%
depth=ncread(filename,'moc_z');depth_km=depth./1e5;

%------(2) from NCL mask-------%
load([pathi3,'TOPO_mask'])
load([pathi1,'MOC_spatial']);for ii=1:6,MOC_test1{ii}=MOC_every10{ii}.*mask;end
load([pathi2,'MOC_spatial']);for ii=1:6,MOC_vdc{ii}=MOC_every10{ii}.*mask;end

for n=1:6 ; 
diff{n} = MOC_vdc{n}-MOC_test1{n};
end

MOC=diff;
set(0,'defaultfigurecolor','w')
[x2,y2]=meshgrid(lat_depth,depth_km);

time_labels = {'Difference [1958-1967]',...
               'Difference [1968-1977]',...
               'Difference [1978-1987]',...
               'Difference [1988-1997]',...
               'Difference [1998-2007]',...
               'Difference [2009-2018]'};
[ha, hfigure]= tight_subplot(2, 3, [.1 .05], [.25 .08], [.08 .05]);

for ii=1:6
  axes(ha(ii));
s1{ii}=MOC{ii};s2{ii}=MOC{ii};
s1{ii}(s1{ii}<0)=0;s2{ii}(s2{ii}>0)=0;

  pcolor(x2,y2,MOC{ii}');shading interp;hold on
  [C,H]=contourf(x2,y2,MOC{ii}',[-5:.2:5],'Showtext','on','fill','off','Linestyle','-','LabelSpacing',500);hold on;
  clabel(C,H,[-5:.4:5],'Fontsize',14,'fontweight', 'bold')

  % caxis([-5 5]);
  caxis([-3 3]);
  c=colorbar('location','southoutside','Position', [0.15 .075 0.75 .02],...
           'colormap', colormap(flipud(brewermap(30, 'RdBu'))))
  c.Title.String = "SV";

if ii==1,
    [C,H]=contourf(x2,y2,MOC{ii}',[-5:.05:5],'Showtext','on','fill','off','Linestyle','-','LabelSpacing',500);hold on;
    clabel(C,H,'Fontsize',14,'fontweight', 'bold')
end

if ii==1|| ii==4,ylabel({'depth [km]'});end
if ii==4|| ii==5|| ii==6,xlabel('lat');end

set(gca,'YDir','reverse')
hold on;xlim([-30 90]);
set(gca,'Fontsize',16,'Linewidth',1.5,'XMinorTick','on',...
    'YMinorTick','on','TickDir','out','fontweight', 'bold');
title({time_labels{ii}},'Interpreter','none');
colormap(flipud(brewermap(30, 'RdBu')));

end
