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

%========== Model output ===========%
filename="/data3/yiwen/cesm2/JGR_ocean/Model_output/G_JRA55_test1.pop.h.0031-04.nc"
pathi='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/'

file=[pathi,'TSRyr_000101-006112.nc']
LAT=ncread(file,'TLAT');lat=LAT(1,:);        
LON=ncread(file,'TLONG');lon=LON(:,1);       
depth=ncread(file,'z_t');depth_m=depth./100; %unit=cm


salt=ncread(file,'SALT');temp=ncread(file,'TEMP');
mask=temp;mask(~isnan(mask))=1; %1=ocean & nan=land
%====[our study] black box ====%
loc1=find(LON<=360-37.5 & LON>=360-55 & LAT<=57.5 & LAT>=47.5);
depthIN=find(depth_m<=300);

%-------- volume mean ----------%
area=ncread(filename,'TAREA'); 
for k=1:60
areamask{k}=area.*mask(:,:,k).*1e-4;    %2D
areamask_IN{k}=areamask{k}(loc1); %1D
sumarea(k)=squeeze(nansum(areamask_IN{k}));
%--------------%
 for count=1:61
 temp_k{k}{count}=temp(:,:,k,count);
 salt_k{k}{count}=salt(:,:,k,count);
 temp_IN{k}{count}=temp_k{k}{count}(loc1);
 salt_IN{k}{count}=salt_k{k}{count}(loc1);

 temp_ls(k,count)=nansum(nansum(temp_IN{k}{count}.*areamask_IN{k}))/sumarea(k);
 salt_ls(k,count)=nansum(nansum(salt_IN{k}{count}.*areamask_IN{k}))/sumarea(k);


 pd_ls(k,count) = sw_dens0(salt_ls(k,count),temp_ls(k,count));  %kg/m^3
 end
end

temp_ls=temp_ls(depthIN,:);
salt_ls=salt_ls(depthIN,:);
pd_ls=pd_ls(depthIN,:);


m_map=colormap(flipud(jet(size(temp_ls,1))));
Msize=length(sumarea);
%------ TS diagram ------%
set(0,'defaultfigurecolor','w')

theta_sdiag(temp_ls,salt_ls);hold on;
for k=1:size(salt_ls,1)
scatter(salt_ls(k,:),temp_ls(k,:),Msize,pd_ls(k,:),'MarkerFaceColor',m_map(k,:),'MarkerEdgeColor','none')
end
colormap(m_map);colorbar
axis square;ytickformat('%.2f')
set(gca,'Fontsize',20,'Linewidth',3,'XMinorTick','on','YMinorTick','on','TickDir','out','fontweight', 'bold');
cbh = colorbar('XTickLabel',{'0','100','200','300'}, ...
               'XTick', 0:20:length(sumarea),'location','eastoutside')
cbh.LineWidth = 2;cbh.FontSize = 18;  
cbh.Title.String='depth(m)'
caxis([0 length(sumarea)]);
text(34.7,6.2,{'Salinity'; 'dominant'},'Fontsize',20,'fontweight', 'bold','Color',[0.4660 0.6740 0.1880])
text(34.9,3,{'Temperature'; 'dominant'},'Fontsize',20,'fontweight', 'bold','Color','b') 
title('KPP_{ICEL} T-S diagram','Interpreter','tex')


