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"

pathi='/data3/yiwen/cesm2/JGR_ocean/Model_output/CTRL/'
% pathi='/data3/yiwen/cesm2/JGR_ocean/Model_output/KPP_linear/'

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

%------ every 10-yrs -------% 
n={1:10 11:20 21:30 31:40 41:50 52:61}  

for ii=1:size(n,2)
salt_M{ii}=squeeze(nanmean(salt(:,:,:,n{ii}),4));
temp_M{ii}=squeeze(nanmean(temp(:,:,:,n{ii}),4));
pd_M{ii}=sw_dens0(salt_M{ii},temp_M{ii});
end

%------ Interplation ------%
depth=ncread(filename,'z_t');depth_km=depth./1e5; 
lat_real=-80:1:90;
lon_real=0:1:360;
[xq,yq] = meshgrid(lon_real,lat_real);

V={salt_M temp_M}
for ii=1:size(V,2)
 for n=1:6
  for k=1:60
V_N{ii}{n}(:,:,k)=griddata(LON,LAT,V{ii}{n}(:,:,k),xq,yq);
  end
V_N2{ii}{n}=circshift(V_N{ii}{n},-180,2);
V_N2{ii}{n}(:,181:183,:)=[];
 end
end

lon2=linspace(-180,180,size(V_N2{1}{1},2));
%====(1) W to E ====%
lonIN=find(lon2>=-60 & lon2<=20); %LS+ overflow

for n=1:6
SALT_M{n}=squeeze(nanmean(V_N2{1}{n}(:,lonIN,:),2));
TEMP_M{n}=squeeze(nanmean(V_N2{2}{n}(:,lonIN,:),2));

PD_M{n}=sw_dens0(SALT_M{n},TEMP_M{n});
end

%----- save to mat -------%
patho=pathi;

fn1 = [patho 'SALT_WE.mat'];
save(fn1,'SALT_M','PD_M','-v7.3')


