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 ------%
lat2=-80:1:90; %1X1
depth=ncread(filename,'z_t');depth_km=depth./1e5; 
%--------(1) MOC data ---------%
load([pathi1,'MOC_spatial']);MOC_test1=MOC_every10;  
load([pathi2,'MOC_spatial']);MOC_vdc=MOC_every10;
for ii=1:6,MOC_diff{ii}=MOC_vdc{ii}-MOC_test1{ii};end
%------(2)-MOC mask ------%
load([pathi3,'TSD_lat_mask_1X1']);

depth2=linspace(depth_moc(1),depth_moc(end),60);
[x1,k1]=meshgrid(lat_depth,depth_moc);x1=double(x1);k1=double(k1);
[xq,kq]=meshgrid(lat2,depth);xq=double(xq);kq=double(kq);
for ii=1:6
MOC_N{ii}=griddata(x1,k1,MOC_diff{ii}',xq,kq);
MOC{ii}=MOC_N{ii}'.*mask;
end
%--------(2) SALT & PD ---------%
load([pathi1,'SALT_WE.mat']);      
SALT_test1=SALT_M;PD_test1=PD_M;
load([pathi2,'SALT_WE.mat']);
SALT_vdc=SALT_M;PD_vdc=PD_M;
for ii=1:6
SALT_diff{ii}=(SALT_vdc{ii}-SALT_test1{ii});
PD_diff{ii}=(PD_vdc{ii}-PD_test1{ii});
end
set(0,'defaultfigurecolor','w')

depth_m=depth./100;
[xx,yy]=meshgrid(lat2,depth_m);

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


for ii=1:6
  axes(ha(ii)); 
[C,h]=contourf(xx,yy,SALT_diff{ii}',[-.1:.01:.1],'ShowText','off','Fill','on');

set(h,'LineColor','none')
cRange=caxis;drawnow;hold on
set(gca,'Fontsize',16,'Linewidth',3,'XMinorTick','on',...
    'YMinorTick','on','TickDir','out','fontweight', 'bold');
set(gca,'YDir','reverse')
xlim([30 89])
%------------------------------------------------%
step=.2; 
[C1,H1]=contourf(xx,yy,MOC{ii}',[-5:step:5],'ShowText','on','Fill','off');
clabel(C1,H1,[0:.4:4],'Fontsize',10,'fontweight', 'bold')
%----- 27.8 /27.95 -----%    
[C2,H2]=contourf(xx,yy,PD_test1{ii}'-1000,[27.8 27.95],...
            'LineColor','b','LineWidth',1.5,'Fill','off','Showtext','off','LabelSpacing',400);
[C3,H3]=contourf(xx,yy,PD_vdc{ii}'-1000,[27.8 27.95],...
            'LineColor','g','LineWidth',2.0,'Fill','off','Showtext','on','LabelSpacing',400); 
clabel(C3,H3,'Fontsize',14,'fontweight', 'bold','Color',[0.4660 0.6740 0.1880])

caxis(cRange);hold on
set(gca,'Fontsize',16,'Linewidth',3,'XMinorTick','on',...
    'YMinorTick','on','TickDir','out','fontweight', 'bold');
set(gca,'YDir','reverse')
xlim([30 89])
text(59,3400,'GISR','Fontsize',12,'fontweight', 'bold')
text(68,3800,'KPP_{ICEL} - CTRL','Fontsize',10,'fontweight', 'bold','Interpreter','tex')
%-------------------------------------------%
if ii==1|| ii==4,ylabel({'depth [km]'});end
if ii==4|| ii==5|| ii==6,xlabel('lat');end

title({time_labels{ii}},'Interpreter','none');  


caxis([-.08 .08]); 
c=colorbar('location','southoutside','Position', [0.15 .075 0.75 .02],...
           'colormap', colormap(flipud(brewermap(30, 'RdBu'))))
c.Title.String = "Psu";
ylim([0 4000])
end