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

%========== Model output ===========%
filename="/data3/yiwen/cesm2/JGR_ocean/Model_output/G_JRA55_test1.pop.h.0031-04.nc"
LAT=ncread(filename,'TLAT');LON=ncread(filename,'TLONG');    
depth=ncread(filename,'z_t');                 %unit=cm

%========= EXP ==========%
for Model_load=1
pathi1='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/';      %CTRL
pathi2='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/' %KPP_ICEL
pathi3='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_R/'      %KPP_R  

file1=[pathi1,'TSRyr_000101-006112.nc']
file2=[pathi2,'TSRyr_000101-006112.nc']
file3=[pathi3,'TSRyr_000101-006112.nc']

% choose diff depth 
layer = [8 24]; %75m/ 250m
n=48:57;        %2005-2014[WOA time]

%-------------[time avg]--------------------%
salt_test1={};temp_test1={};
salt_ice={};temp_ice={};
salt_ice01={};temp_ice01={};
for ii=1:2
S1=ncread(file1,'SALT');salt_test1{ii}=squeeze(nanmean(S1(:,:,layer(ii),n),4));
T1=ncread(file1,'TEMP');temp_test1{ii}=squeeze(nanmean(T1(:,:,layer(ii),n),4));
clear S1 T1
rho_test1{ii}=sw_dens(salt_test1{ii},temp_test1{ii},0); %kg/m^3
%--- linear 
S2=ncread(file2,'SALT');salt_ice{ii}=squeeze(nanmean(S2(:,:,layer(ii),n),4));
T2=ncread(file2,'TEMP');temp_ice{ii}=squeeze(nanmean(T2(:,:,layer(ii),n),4));
clear S2 T2
rho_ice{ii}=sw_dens(salt_ice{ii},temp_ice{ii},0); %kg/m^3
%--- exp 
S3=ncread(file3,'SALT');salt_ice01{ii}=squeeze(nanmean(S3(:,:,layer(ii),n),4));
T3=ncread(file3,'TEMP');temp_ice01{ii}=squeeze(nanmean(T3(:,:,layer(ii),n),4));
clear S3 T3
rho_ice01{ii}=sw_dens(salt_ice01{ii},temp_ice01{ii},0); %kg/m^3
end
end

%===== {layer}*{EXP} =====%
for ii=1:length(layer)
Var{ii}={salt_ice{ii} salt_ice01{ii} salt_test1{ii}...
         temp_ice{ii} temp_ice01{ii} temp_test1{ii}...
          rho_ice{ii}  rho_ice01{ii}  rho_test1{ii}};

%===== difference[salt(1-3)/temp(4-6)] =====%
Var_diff{ii}={Var{ii}{1}-Var{ii}{3} Var{ii}{2}-Var{ii}{3}...
              Var{ii}{1}-Var{ii}{2}...
              Var{ii}{4}-Var{ii}{6} Var{ii}{5}-Var{ii}{6}...
              Var{ii}{4}-Var{ii}{5}...
              Var{ii}{7}-Var{ii}{9} Var{ii}{8}-Var{ii}{9}...
              Var{ii}{7}-Var{ii}{8}};
end


%------- interplation -------%
lat_real=-90:1:90;
lon_real=0:1:360;
[xq,yq] = meshgrid(lon_real,lat_real);

for N=1:size(Var_diff{1},2)  
  for ii=1:size(Var_diff,2)  
Var_N{ii}{N}=griddata(LON,LAT,Var_diff{ii}{N},xq,yq); %change

Var_N2{ii}{N}=circshift(Var_N{ii}{N},-180,2);
%-------1X1-------%
Var_N2{ii}{N}(:,181:183)=[];
 end
end
%====== study region ======%
Xw=[-55 -55 -37.5 -37.5 -55]
Yw=[47.5 57.5 57.5 47.5 47.5]

%--range_total--%
lon2=linspace(-180,180,size(Var_N2{1}{1},2));
lat2=linspace(-90,90,size(Var_N2{1}{1},1));
lon_loc=lon2>=-180 & lon2<=180;
lat_loc=lat2>=-90 & lat2<=90;
[xx,yy]=meshgrid(lon2(lon_loc),lat2(lat_loc));

n=2;  %layer[75/250 m]
%---- contour line ----%
N=Var_N2{n};
P=Var_N2{n};
C_range=[-.5:.05:.5];


for ii=1:9 
N{ii}(N{ii}>0)=nan;
P{ii}(P{ii}<0)=nan;
end

set(0,'defaultfigurecolor','w')
[ha, hfigure]= tight_subplot(2, 3, [.07 .01], [.06 .06], [.08 .1]);

for ii=1:6
  axes(ha(ii));  
m_proj('azimuthal equal-area','latitude',80,'radius',80,'rectbox','on');%45/60
m_pcolor(xx,yy,Var_N2{n}{ii}(lat_loc,lon_loc));hold on  
m_coast('patch',[.7 .7 .7],'edgecolor','k');
m_grid('xtick',9,'tickdir','out','ytick',[20 45 65 75 85],'Fontsize',10);
m_line(Xw,Yw,'linewi',2,'color','b');hold on
if ii==1 | ii==4
    indx=7;
    [cs,h]=m_contour(xx,yy,P{indx}(lat_loc,lon_loc),C_range,'color','k',...
           'LineWidth',1.5,'LineStyle','-','ShowText','off');hold on 
    [cs2,h2]=m_contour(xx,yy,N{indx}(lat_loc,lon_loc),C_range,'color','k',...
           'LineStyle','--','ShowText','off');hold on
    text(-.2,-.83,{'KPP_{ICEL} minus CTRL'},'Fontsize',12,...
         'fontweight', 'bold','Interpreter','tex')
elseif ii==2 | ii==5
    indx=8;
    [cs,h]=m_contour(xx,yy,P{indx}(lat_loc,lon_loc),C_range,'color','k',...
           'LineWidth',1.5,'LineStyle','-','ShowText','off');hold on 
    [cs2,h2]=m_contour(xx,yy,N{indx}(lat_loc,lon_loc),C_range,'color','k',...
           'LineStyle','--','ShowText','off');hold on
    text(-.2,-.83,{'KPP_{R} minus CTRL'},'Fontsize',12,...
    'fontweight', 'bold','Interpreter','tex')
else
    indx=9;
    [cs,h]=m_contour(xx,yy,P{indx}(lat_loc,lon_loc),C_range,'color','k',...
           'LineWidth',1.5,'LineStyle','-','ShowText','off');hold on 
    [cs2,h2]=m_contour(xx,yy,N{indx}(lat_loc,lon_loc),C_range,'color','k',...
           'LineStyle','--','ShowText','off');hold on
    text(-.2,-.83,{'KPP_{ICEL} minus KPP_{R}'},'Fontsize',12,...
     'fontweight', 'bold','Interpreter','tex')
end
if ii==1,y=text(-1.2,-.3,'SALT diff','Fontsize',20,'fontweight', 'bold')
set(y,'Rotation',90);end
if ii==4,y=text(-1.2,-.3,'TEMP diff','Fontsize',20,'fontweight', 'bold')
set(y,'Rotation',90);end

if ii <= 3
    caxis([-.45 .45]); %250m
else
    caxis([-.6 .6]);
end


if ii==3
c=colorbar('location','eastoutside','Position', [.9 .55 .015 .35],...
         'colormap', colormap(flipud(brewermap(50, 'RdBu'))))
c.Title.String = {"[250m diff.]";"psu"};
c.FontSize = 14;
caxis([-.45 .45])
end

if ii==6
c=colorbar('location','eastoutside','Position', [.9 .08 .015 .35],...
         'colormap', colormap(flipud(brewermap(50, 'RdBu'))))
c.Title.String = {"degree"};
c.FontSize = 14;
caxis([-.6 .6])
end
end





