%% get  MV-EOF TZ SIC WeTz
% copyright@ Fan Jia
clear
kjr=1:31;

% load 3D temperature of 31 models
load D:\southern_ocean\Tz_hr_deeper_Z_around100m_new_31.mat
Tz_h1=cat(2,Tz_h(:,22:end,:),Tz_h(:,1:20,:));
% load Ekman pumping of 31 models
load D:\southern_ocean\Ekman_pump_tau_hr_new_26.mat
load D:\southern_ocean\Ekman_pump_tau_hr_new_5f2.mat
EK_h31=cat(3,EK_h,EK_hf2);
EK_r31=cat(3,EK_r,EK_rf2);

% get heating of upwelling
dweTz=-(EK_r31-EK_h31).*permute(Tz_h1(1:180,1:359,:),[2 1 3])/60;
for i=1:31
    aa=fillnan(dweTz(:,:,i)');
    dweTz1(:,:,i)=aa';
end

% normalized by GW
load D:\southern_ocean\gw_26_new.mat
load D:\southern_ocean\gw_5f2_new.mat
gw31=cat(2,gw,gw_f2); clear gw gw_f2
dweTz2=dweTz1./permute(repmat(gw31',[1 359 180]),[2 3 1]);

load D:\southern_ocean\land_new.mat
dweTz3=dweTz2(:,:,:).*repmat(land(1:359,1:180),[1 1 31]);

% load SIC of 31 models
load D:\southern_ocean\SICchange_original_new.mat
load D:\southern_ocean\SICchange_original_new_5f2.mat
SICchange_new31=cat(3,SICchange_new,dsic); clear SICchange_new dsic
dsic=SICchange_new31./permute(repmat(gw31',[1 360 181]),[3 2 1]);
dsic1=cat(2,dsic(:,21:359,:),dsic(:,1:20,:));
for i=1:31
    aa=fillnan(dsic1(:,:,i));
    dsic2(:,:,i)=aa';
end
dsic3=dsic2(:,:,:).*repmat(land(1:359,1:181),[1 1 31]);


% load shelf water warming
load D:\southern_ocean\1116\thetao_hk.mat
load D:\southern_ocean\1116\thetao_rk.mat

dthetao0=thetao_rk1(13:end,1:151,:)-thetao_hk1(13:end,1:151,:);

LEV=[1 10:10:6000];
LEVV=LEV(1:151);
LATT=-8:1:10;

dthetao=dthetao0./permute(repmat(gw36',[1 length(LEVV) length(LATT)]),[3 2 1]);
aaa=mean(dthetao,3);
aaa(~isnan(aaa))=1;
dthetao11=dthetao.*repmat(aaa,[1 1 31]);
dthetao1=dthetao11(:,:,kjr);
dthetao01=dthetao(:,:,kjr);


lat=-90:1:89;
lon=21:1:379;
lat1=find(lat==-90);
lat2=find(lat==-50);
lon1=find(lon==21);
lon2=find(lon==379);
dweTz4=dweTz3(lon1:lon2,lat1:lat2,kjr);



lat=-90:1:89;
lon=21:1:379;
lat1=find(lat==-90);
lat2=find(lat==-50);
lon1=find(lon==21);
lon2=find(lon==379);
dsic4=dsic3(lon1:lon2,lat1:lat2,kjr);


dweTz5=reshape(dweTz4,[size(dweTz4,1)*size(dweTz4,2) size(dweTz4,3)]);
dsic5=reshape(dsic4,[size(dsic4,1)*size(dsic4,2) size(dsic4,3)]);
dsst5=reshape(dthetao1,[size(dthetao1,1)*size(dthetao1,2) size(dthetao1,3)]);


% get ocean grid only
sea=find(~isnan(dweTz5(:,1)));
land=find(isnan(dweTz5(:,1)));
dweTz6=dweTz5(sea,:);

sea1=find(~isnan(dsic5(:,1)));
land1=find(isnan(dsic5(:,1)));
dsic6=dsic5(sea1,:);


sea2=find(~isnan(dsst5(:,1)));
land2=find(isnan(dsst5(:,1)));
dsst6=dsst5(sea2,:);

% get normalized
dweTz7=nor_2(dweTz6);
dsic7=nor_2(dsic6);
dsst7=nor_2(dsst6);

dsic7(dsic7==inf)=nan;
dsic7(dsic7==-inf)=nan;
sea11=find(~isnan(dsic7(:,1)));
land11=find(isnan(dsic7(:,1)));
dsic77=dsic7(sea11,:);

sea22=find(~isnan(dweTz7(:,1)));
land22=find(isnan(dweTz7(:,1)));
dweTz77=dweTz7(sea22,:);

%% get MV-EOF
x=[dweTz77;dsic77;dsst7];

onoff=0; % indicate whether needed to perform time-space transition

[m, n]=size(x); % m: number of spacial grid points, n: number of temporal points
if m>n
    S=x'*x/m;
    onoff=1;
else
    S=x*x'/n;
end

[L,lmd]=eig(S, 'nobalance');
lmd=diag(lmd); % lmd is a column vector here.
lmd=lmd';      % lmd: a row vector here.
if onoff==1
%     L=(x*L)/((m*diag(lmd)).^0.5);   %归一化
    L=(x*L)*diag(1.0./sqrt(m*lmd));
    lmd=m/n*lmd;
end

standardized_flag=1;
M=length(kjr);
L=fliplr(L);
T=L'*x;
lmd=fliplr(lmd);
VC(1:M)=lmd(1:M)/sum(lmd);
CVC=zeros(1, M);
for i=1:M
    CVC(i)=sum(VC(1:i));
end
lmd=lmd(1:M);
L=L(:,1:M)*diag(sqrt(lmd));
if (standardized_flag==1)
   T=diag(1./sqrt(lmd))*T(1:M,:); %标准化   
elseif (standardized_flag==0)
   T=T(1:M,:); %不标准化 
end

% T is pc

%% get regression
pc1=T(1,:);
pc11=-(pc1-mean(pc1))/std(pc1);
% pc11=pc1;

lat=-90:1:89;
lon=21:1:379;
lat1=find(lat==-90);
lat2=find(lat==-50);
lon1=find(lon==21);
lon2=find(lon==379);

dweTz9=dweTz3(lon1:lon2,lat1:lat2,kjr);
dsic9=dsic3(lon1:lon2,lat1:lat2,kjr);


dweTz10=cat(1,dweTz9,dweTz9(end,:,:),dweTz9(1,:,:));
dsic10=cat(1,dsic9,dsic9(end,:,:),dsic9(1,:,:));


[f00,fa00,fsig00]=J_regress(dweTz10,pc11,0.05);
ff00=f00.*fsig00;
ff00(ff00==0)=nan;
[f11,fa11,fsig11]=J_regress(dsic10,pc11,0.05);
ff11=f11.*fsig11;
ff11(ff11==0)=nan;
[f22,fa22,fsig22,cpp]=J_regress_ice(dthetao01(:,1:end,:),pc11,0.05);
ff22=f22.*fsig22;
ff22(ff22==0)=nan;

