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{6};  
load([pathi2,'MOC_spatial']);MOC_kpp=MOC_every10{6};

MOC_diff = MOC_kpp-MOC_test1 ;
%--------(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);

MOC_N = griddata(x1,k1,MOC_diff',xq,kq);
MOC  = MOC_N'.*mask;

%--------(2) SALT & PD ---------%
load([pathi1,'SALT_WE.mat']);      
SALT_test1=SALT_M{6};PD_test1=PD_M{6};
load([pathi2,'SALT_WE.mat']);
SALT_kpp=SALT_M{6};PD_kpp=PD_M{6};

SALT_diff={SALT_kpp-SALT_test1};
PD_diff  ={PD_kpp-PD_test1 };

depth_m=depth./100;
[xx,yy]=meshgrid(lat2,depth_m);
set(0,'defaultfigurecolor','w')

[C,h]=contourf(xx,yy,SALT_diff{1}',[-3:.01:3],'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',[-5:step:5],'ShowText','on','Fill','off');
clabel(C1,H1,[0:.4:4],'Fontsize',15,'fontweight', 'bold')
%------------------------------------------------%
step=.2; 
%-----[isopycnal line] 27.8 /27.95 -----%    
[C2,H2]=contourf(xx,yy,PD_test1'-1000,[27.8 27.95],...
            'LineColor','b','LineWidth',1.5,'Fill','off','Showtext','off');
%-----------------------%

[C3,H3]=contourf(xx,yy,PD_kpp'-1000,[27.8 27.95],...
        'LineColor','g','LineWidth',2.0,'Fill','off','Showtext','on','LabelSpacing',400)
text(60,5000,'KPP_{ICEL} - CTRL','Fontsize',18,'fontweight', 'bold','Interpreter','tex')
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(60,3400,'GISR','Fontsize',14,'fontweight', 'bold')
%-------------------------------------------%
xlabel('lat');ylabel('Depth(m)')
caxis([-.1 .1]); 
c=colorbar('location','eastoutside','Position', [0.75 .12 .02 .75],...
       'colormap', colormap(flipud(brewermap(30, 'RdBu'))))
c.Title.String = "Psu"  
axis square

