%% MLD temperature 
clear
clc

f1='/Volumes/Backup Plus/PAROCE/mld22.nc';
Depth = ncread(f1, 'deptht');%75 levels
Lat = ncread(f1, 'nav_lat');
Lon = ncread(f1, 'nav_lon');
index{1}=find(Lon >= -60 & Lon <= 20); %weddell
index{2}=find(Lon >= 20 & Lon <= 130); %indian ocean
index{3}=find(not(Lon >= -160 & Lon <= 130)); %wpo
index{4}=find(Lon >= -160 & Lon <= -130); %eastern ross
index{5}=find(Lon >= -130 & Lon <= -60); %abs
% index{6}=find(Lon>-180&Lon<-120&Lat>-80&Lat<-70); %polynya
index{6}=1:(1442*559); %SO

% f='/Volumes/Backup Plus/PAROCE/nemo_prod_b_atm_mod_2021_2021_1m_20221001_20221231_grid_T.nc';
f1='/Volumes/Backup Plus/PAROCE/mld22.nc';
f2='/Volumes/Backup Plus/PAROCE/mldclim.nc';
f3='/Volumes/Backup Plus/PAROCE/mld1521.nc';

% MLD22 = ncread(f1, 'mldr10_1',[1 1 1],[1442 559 12],[1 1 1]);%mld defined by sigma theta
% Lat = ncread(f1, 'nav_lat');
% Lon = ncread(f1, 'nav_lon');

PTc = ncread(f2,'thetao');%potential temperature
% MLDc = ncread(f2, 'mldr10_1');%mld defined by sigma theta

k=1;
for i=1:12:73
    
    PTa=ncread(f3,'thetao',[1 1 1 i],[1442 559 75 12],[1 1 1 1])-PTc;%potential temperature
    
    for t=1:12
        for d=1:35
            p0=PTa(:,:,d,t);
            for n=1:6
                tm(d,k,n)=nanmean(p0(cell2mat(index(n))));
            end
        end
        k=k+1;
    end
    
end

PTa=ncread(f1,'thetao',[1 1 1 1],[1442 559 75 12],[1 1 1 1])-PTc;%potential temperature
for t=1:12
    for d=1:35
        p0=PTa(:,:,d,t);
        for n=1:6
            tm(d,k,n)=nanmean(p0(cell2mat(index(n))));
        end
    end
    k=k+1;
end

PTa=ncread(f1,'thetao',[1 1 1 13],[1442 559 75 2],[1 1 1 1])-PTc(:,:,:,1:2);%potential temperature
for t=1:2
    for d=1:35
        p0=PTa(:,:,d,t);
        for n=1:6
            tm(d,k,n)=nanmean(p0(cell2mat(index(n))));
        end
    end
    k=k+1;
end


% MLD0 = ncread(f3, 'mldr10_1');%mld defined by sigma theta


% for t=1:12
% for d=1:25
%     p0=PTa(:,:,d,t);
%     tm(d,t)=mean(p0(Lon>-180&Lon<-120&Lat>-80&Lat<-70));
% end
% end
% tm0=flipud(tm);

% for t=1:12
%     m0=MLDa(:,:,t);
%     dm(t)=mean(m0(Lon>-180&Lon<-120&Lat>-80&Lat<-70));
% end

%%
clear
clc

load oceanano
f1='/Volumes/Backup Plus/PAROCE/mld22.nc';
Depth = ncread(f1, 'deptht');%75 levels
Lat = ncread(f1, 'nav_lat');
Lon = ncread(f1, 'nav_lon');
index{1}=find(Lon >= -60 & Lon <= 20); %weddell
index{2}=find(Lon >= 20 & Lon <= 130); %indian ocean
index{3}=find(not(Lon >= -160 & Lon <= 130)); %wpo
index{4}=find(Lon >= -160 & Lon <= -130); %eastern ross
index{5}=find(Lon >= -130 & Lon <= -60); %abs

% index{6}=find(Lon>-180&Lon<-120&Lat>-80&Lat<-70); %polynya
index{6}=1:(1442*559); %SO

for t=1:12
for d=1:35
    p0=PTa(:,:,d,t);
    for n=1:6
    tm(d,t,n)=nanmean(p0(cell2mat(index(n))));
%     tm(d,t)=nanmean(p0(:));
    end
end
end

% for t=1:12
%     m0=MLDa(:,:,t);
%     for n=1:6
%     dm(t,n)=nanmean(m0(cell2mat(index(n))));
% %     dm(t)=nanmean(m0(:));
%     end
% end



%%

% t=1:98;
t=1:12;
d=Depth(1:35);

ti=["WS","IO","WPO","ERS","ABS","SO"];

set(gcf,'Position',[350 200 1100 600]);
ha = tight_subplot(2,3,[0.08 0.05],[.1 .04],[.05 .03]);

for n=1:6
axes(ha(n))

[M,c]=contourf(t,d,tm(:,:,n),15);
c.LineWidth = 1;
hold on
set(gca,'YDir','reverse')
% plot(dm(:,n),'k','LineWidth',2)
ylim([0 300])
ylabel('Depth (m)','fontname','Arial','fontsize',12)
% set(gca,'xtick',1:24:98,'xticklabel',{'2015' '2017' '2019' '2021' '2023'},'fontname','Arial','fontsize',12);
set(gca,'xtick',1:3:12,'xticklabel',{'Mar' 'Jun' 'Sep' 'Dec'},'fontname','Arial','fontsize',12);

colortable = textread('/Users/jinfei/Documents/MATLAB/colorbar/ncl_default.txt');
colormap(colortable);
% caxis([-1 1]);
caxis([-0.5 0.5]);
title(ti(n),'fontsize',14)
end

c = colorbar('southoutside','FontSize',10);
set(c,'position',[0.32 0.04 0.4 0.02],'Ytick',[-0.5 0 0.5])
t1=get(c,'Yticklabel');
t1=strcat(t1,'degC');
set(c,'Yticklabel',t1);

