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

%========== 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

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

%=================================%
file1=[pathi1,'TSRyr_000101-006112.nc']
file2=[pathi2,'TSRyr_000101-006112.nc']
salt_test1=ncread(file1,'SALT');temp_test1=ncread(file1,'TEMP');
salt_vdc=ncread(file2,'SALT');temp_vdc=ncread(file2,'TEMP');

Var={salt_test1 temp_test1 salt_vdc temp_vdc};

for ii=1:4
Var_1st{ii}=squeeze(nanmean(Var{ii}(:,:,:,1:10),4));  %1958~1967
Var_2nd{ii}=squeeze(nanmean(Var{ii}(:,:,:,11:20),4)); %1968~1977
Var_3rd{ii}=squeeze(nanmean(Var{ii}(:,:,:,21:30),4)); %1978~1987
Var_4th{ii}=squeeze(nanmean(Var{ii}(:,:,:,31:40),4)); %1988~1997
Var_5th{ii}=squeeze(nanmean(Var{ii}(:,:,:,41:50),4)); %1998~2007
Var_10{ii} =squeeze(nanmean(Var{ii}(:,:,:,52:61),4)); %last-10yrs
end
salt_CTRL={Var_1st{1} Var_2nd{1} Var_3rd{1} Var_4th{1} Var_5th{1} Var_10{1}}
temp_CTRL={Var_1st{2} Var_2nd{2} Var_3rd{2} Var_4th{2} Var_5th{2} Var_10{2}}

salt_VDC={Var_1st{3} Var_2nd{3} Var_3rd{3} Var_4th{3} Var_5th{3} Var_10{3}}
temp_VDC={Var_1st{4} Var_2nd{4} Var_3rd{4} Var_4th{4} Var_5th{4} Var_10{4}}

%==============================%
%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_CTRL,2)
temp_k1{ii}{k}=temp_CTRL{ii}(:,:,k);
temp_k2{ii}{k}=temp_VDC{ii}(:,:,k);
temp_IN1{ii}{k}=temp_k1{ii}{k}(loc1);
temp_IN2{ii}{k}=temp_k2{ii}{k}(loc1);

salt_k1{ii}{k}=salt_CTRL{ii}(:,:,k);
salt_k2{ii}{k}=salt_VDC{ii}(:,:,k);
salt_IN1{ii}{k}=salt_k1{ii}{k}(loc1);
salt_IN2{ii}{k}=salt_k2{ii}{k}(loc1);

TEMP_N1{ii}(k)=nansum(nansum(temp_IN1{ii}{k}.*areamask_IN{k}))/sumarea(k);
TEMP_N2{ii}(k)=nansum(nansum(temp_IN2{ii}{k}.*areamask_IN{k}))/sumarea(k);
SALT_N1{ii}(k)=nansum(nansum(salt_IN1{ii}{k}.*areamask_IN{k}))/sumarea(k);
SALT_N2{ii}(k)=nansum(nansum(salt_IN2{ii}{k}.*areamask_IN{k}))/sumarea(k);

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

end
end

C = {[0 0 1], [1 0 0], [0 0.4470 0.7410], [0.8500 0.3250 0.0980]};
set(0, 'defaultfigurecolor', 'w');

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


for ii=1:6
axes(ha(ii));
ax1 = gca;
line(TEMP_N1{ii}(depth_loc),depth_m(depth_loc),'Color',C{1},...
    'Linewidth',4,'Parent',ax1,'Linewidth',2.5);hold on
line(TEMP_N2{ii}(depth_loc),depth_m(depth_loc),'Color',C{2},...
    'Linewidth',4,'Parent',ax1,'Linewidth',2.5);hold on

set(ax1,'XColor','k','YColor','k')
if ii==4|| ii==5 || ii==6,xlabel('deg C','fontsize',16);end
if ii==1|| ii==4,ylabel('depth(m)','fontsize',16);end
xticklabel_fs = get(ax1,'XTickLabel');
xlim(ax1,[3 7]);ylim(ax1,[0 400])
set(ax1, 'XTick', 3:2:7);
set(ax1, 'XTickLabel', arrayfun(@(x) sprintf('%d', x), 3:2:7, 'UniformOutput', false));
set(gca, 'YDir', 'reverse', 'FontSize', 16, 'LineWidth', 2, 'FontWeight', 'bold', ...
    'XMinorTick', 'on', 'YMinorTick', 'on', 'TickDir', 'out');
ytickformat('%.2f');
grid on;axis square;hold on

%=============================================%
ax2 = axes('Position', get(ax1, 'Position'), 'XAxisLocation', 'top', ...
               'YAxisLocation', 'right', 'Color', 'none');
hold(ax2, 'on');

line(SALT_N1{ii}(depth_loc),depth_m(depth_loc),'LineStyle','--',...
     'Color',C{1},'Parent',ax2,'Linewidth',2.5);hold on
line(SALT_N2{ii}(depth_loc),depth_m(depth_loc),'LineStyle','--',...
     'Color',C{2},'Parent',ax2,'Linewidth',2.5);hold on

if ii==1|| ii==2 || ii==3,xlabel('psu','fontsize',16);end
xticklabel_fs2 = get(ax2,'XTickLabel');
xlim([34.1 35.2]);ylim(ax2,[0 400])
set(gca, 'YDir', 'reverse', 'FontSize', 16, 'LineWidth', 2, 'FontWeight', 'bold', ...
    'XMinorTick', 'on', 'YMinorTick', 'on', 'TickDir', 'out');
ytickformat('%.2f');
grid on;axis square;hold on
set(ax2,'YTick',[])
h=legend({'CTRL','KPP_{ICEL}'},'Location','southwest','Interpreter','tex');
legend('boxoff');set(h,'Fontsize',16);

title({time_labels{ii}});
end


