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

%========== 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');               

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']

salt_test1=ncread(file1,'SALT');temp_test1=ncread(file1,'TEMP');
salt_kpp=ncread(file2,'SALT');temp_kpp=ncread(file2,'TEMP');
salt_kpp_ts=ncread(file3,'SALT');temp_kpp_ts=ncread(file3,'TEMP');

Var={salt_test1 temp_test1...
     salt_kpp temp_kpp...
     salt_kpp_ts temp_kpp_ts};

n=52:61; %last-10yrs
for ii=1:size(Var,2)
Var_10{ii} =squeeze(nanmean(Var{ii}(:,:,:,n),4)); 
end

SALT_V={Var_10{1} Var_10{3} Var_10{5}}
TEMP_V={Var_10{2} Var_10{4} Var_10{6}}

%==============================%
%range of depth(cm->m)
%==============================%
depth_m=depth./100;
depth_loc=depth_m<=410;       %choose depth reange
loc1=find(LON<=360-37.5 & LON>=360-55 & LAT<=57.5 & LAT>=47.5);

area=ncread(filename,'TAREA');%every layer has diff area
temp1=ncread(filename,'TEMP');mask=temp1;mask(~isnan(mask))=1;

for k=1:60
areamask{k}=area.*mask(:,:,k);    %2D
areamask_IN{k}=areamask{k}(loc1); %1D
sumarea(k)=squeeze(nansum(areamask_IN{k}));

for ii=1:size(SALT_V,2)
temp_k{ii}{k}=TEMP_V{ii}(:,:,k);
temp_IN{ii}{k}=temp_k{ii}{k}(loc1);
salt_k{ii}{k}=SALT_V{ii}(:,:,k);
salt_IN{ii}{k}=salt_k{ii}{k}(loc1);

TEMP_N1{ii}(k)=nansum(nansum(temp_IN{ii}{k}.*areamask_IN{k}))/sumarea(k);
SALT_N1{ii}(k)=nansum(nansum(salt_IN{ii}{k}.*areamask_IN{k}))/sumarea(k);

RHO_N1{ii}=sw_dens(SALT_N1{ii},TEMP_N1{ii},0); %kg/m^3
end
end

C = {[0 0 1], [1 0 0], [0.9290 0.6940 0.1250], [0.8500 0.3250 0.0980]};
set(0, 'defaultfigurecolor', 'w');


subplot(1,2,1) 
ax1 = gca;
for ii=1:size(SALT_V,2) 
line(TEMP_N1{ii}(depth_loc),depth_m(depth_loc),'Color',C{ii},...
    'Linewidth',4,'Parent',ax1,'Linewidth',2.5);hold on
end
set(ax1,'XColor','k','YColor','k')
xlabel('deg C','fontsize',16);

xticklabel_fs = get(ax1,'XTickLabel');
set(gca, 'YDir', 'reverse', 'FontSize', 16, 'LineWidth', 2, 'FontWeight', 'bold', ...
    'XMinorTick', 'on', 'YMinorTick', 'on', 'TickDir', 'out');ytickformat('%.2f');
grid on;axis square;hold on;box on;
ylim([0 400])
set(ax1,'xticklabels')
ylabel('Depth(m)')
%=============================================%
ax2 = axes('Position', get(ax1, 'Position'), 'XAxisLocation', 'top', ...
               'YAxisLocation', 'right', 'Color', 'none');
hold(ax2, 'on');
for ii=1:size(SALT_V,2)  
line(SALT_N1{ii}(depth_loc),depth_m(depth_loc),'LineStyle','--',...
     'Color',C{ii},'Parent',ax2,'Linewidth',2.5);hold on
end
xlabel('psu','fontsize',16);
xticklabel_fs2 = get(ax2,'XTickLabel');
xlim([34.3 35.2])
set(gca, 'YDir', 'reverse', 'FontSize', 16, 'LineWidth', 2, 'FontWeight', 'bold', ...
    'XMinorTick', 'on', 'YMinorTick', 'on', 'TickDir', 'out');ytickformat('%.2f');
grid on;axis square;hold on;box on;
ylim([0 400])

set(ax2,'YTick',[])
h=legend({'CTRL','KPP_{ICEL}','KPP_{R}'},'Location','southwest','Interpreter','tex');
legend('boxoff');set(h,'Fontsize',16);


subplot(1,2,2)
for ii=1:size(SALT_V,2)  
line(RHO_N1{ii}(depth_loc),depth_m(depth_loc),...
     'Color',C{ii},'Linewidth',2.5);hold on
end
xlabel('kg/m^{3}','fontsize',16,'Interpreter','tex');
set(gca, 'YDir', 'reverse', 'FontSize', 16, 'LineWidth', 2, 'FontWeight', 'bold', ...
    'XMinorTick', 'on', 'YMinorTick', 'on', 'TickDir', 'out');ytickformat('%.2f');
grid on;axis square;hold on;box on;
ylim([0 400])

title('analyzed region [2009-2018]');
xlim([1026.95 1027.85])

