% scale of Ubar, Ustar and u
file0='All_WA1Jan35_EQ1M40_0.nc';
file='All_WA1Jan35_EQ1M40_5.nc';
lat=double(ncread(file,'lat'));
lev=double(ncread(file,'lev'));
ilev=double(ncread(file,'ilev'));
lon=double(ncread(file,'lon'));
potentialT=permute(repmat((1000./lev).^(2./7),[1,length(lat),length(lon)]),[3,2,1]);
U=double(squeeze(ncread(file,'U')));
UU=double(squeeze(ncread(file,'UU')));
V=double(squeeze(ncread(file,'V')));
VV=double(squeeze(ncread(file,'VV')));
T=double(squeeze(ncread(file,'T'))).*potentialT;
Z=double(squeeze(ncread(file,'Z3')));
OMEGA=double(squeeze(ncread(file,'OMEGA')));
VT=double(squeeze(ncread(file,'VT'))).*potentialT;
UV=double(squeeze(ncread(file,'VU')));

U0=double(squeeze(ncread(file0,'U')));
UU0=double(squeeze(ncread(file0,'UU')));
V0=double(squeeze(ncread(file0,'V')));
VV0=double(squeeze(ncread(file0,'VV')));
T0=double(squeeze(ncread(file0,'T'))).*potentialT;
Z0=double(squeeze(ncread(file0,'Z3')));
OMEGA0=double(squeeze(ncread(file0,'OMEGA')));
VT0=double(squeeze(ncread(file0,'VT'))).*potentialT;
UV0=double(squeeze(ncread(file0,'VU')));

%% [U,V]{bar,star}, [u,v]
Ubar=squeeze(mean(U));
Ubar0=squeeze(mean(U0));
Ustar=sqrt(squeeze(squeeze(mean(U.^2))-Ubar.^2));
Ustar0=sqrt(squeeze(squeeze(mean(U0.^2))-Ubar0.^2));
u=sqrt(squeeze(mean(UU-U.^2)));

Vbar=squeeze(mean(V));
Vbar0=squeeze(mean(V0));
Vstar=sqrt(squeeze(squeeze(mean(V.^2))-Vbar.^2));
Vstar0=sqrt(squeeze(squeeze(mean(V0.^2))-Vbar0.^2));
v=sqrt(squeeze(mean(VV-V.^2)));

Tbar=squeeze(mean(T));
Tbar0=squeeze(mean(T0));
Tstar=sqrt(squeeze(squeeze(mean(T.^2))-Tbar.^2));
Tstar0=sqrt(squeeze(squeeze(mean(T0.^2))-Tbar0.^2));

dUbar=Ubar-Ubar0;
dUstar=sqrt(squeeze(mean((U-U0).^2))-dUbar.^2);
dVbar=Vbar-Vbar0;
dVstar=sqrt(squeeze(mean((V-V0).^2))-dVbar.^2);
dTbar=Tbar-Tbar0;
dTstar=sqrt(squeeze(mean((T-T0).^2))-dTbar.^2);

U_V=struct('Ubar0',Ubar0,'Ustar0',Ustar0,'u',u,...
    'Vbar0',Vbar0,'Vstar0',Vstar0,'v',v,...
    'dUbar',dUbar,'dUstar',dUstar,'dVbar',dVbar,'dVstar',dVstar);
varnames=fieldnames(U_V);
for iplt=1:length(varnames)
figure(20+iplt)
clf
colormap(jet(20))
contourf(lat,lev,(getfield(U_V,char(varnames(iplt)))'),[0:1:20])
caxis([0,20])
title(varnames(iplt))
set(gca,'yscale','log')
set(gca,'ydir','reverse')
xlabel('lat')
ylabel('lev')
xlim([-15,45])
ylim([100,1000]);
set(gca,'FontSize',18)
colorbar
end

%% Terms_U0
Re=6371000.;
f0=2*2*pi/86400*sin(10/180*pi);
beta=2*2*pi/86400/Re;
g=9.8;
tau=40.*86400;
y=lat*Re/180*pi;
pbot=900.;
ptop=120.;
[~,ipbot]=min(abs(lev-pbot));
[~,iptop]=min(abs(lev-ptop));
dp=repmat(ilev(2:end)-ilev(1:end-1),[1,length(lat)])';
dp=dp(:,iptop:ipbot);
pltcontour=false;

TU0_fdV0=repmat(beta*y,[1,length(lev)]).*dVbar;
TU0_fdV0_int=sum(TU0_fdV0(:,iptop:ipbot).*dp,2);
[~,TU0_uv]=gradient(-squeeze(mean((UV-U.*V)-(UV0-U0.*V0))),lev,y);
TU0_uv_int=sum(TU0_uv(:,iptop:ipbot).*dp,2);
TU0_dU0_tau=-dUbar./tau;
TU0_dU0_tau_int=sum(TU0_dU0_tau(:,iptop:ipbot).*dp,2);
[~,Uy_temp]=gradient(Ubar0,lev,y);
TU0_dV0_U0=-dVbar.*Uy_temp;
TU0_dV0_U0_int=sum(TU0_dV0_U0(:,iptop:ipbot).*dp,2);
[~,Uy_temp]=gradient(dUbar,lev,y);
TU0_V0_dU0=-Vbar0.*Uy_temp;
TU0_V0_dU0_int=sum(TU0_V0_dU0(:,iptop:ipbot).*dp,2);

Terms_U0=struct('fdV0',TU0_fdV0,'uv',TU0_uv,'dU0_tau',TU0_dU0_tau,...
    'dV0_U0',TU0_dV0_U0,'V0_dU0',TU0_V0_dU0);
TermsInt_U0=struct('fdV0',TU0_fdV0_int,'uv',TU0_uv_int,'dU0_tau',TU0_dU0_tau_int,...
    'dV0_U0',TU0_dV0_U0_int,'V0_dU0',TU0_V0_dU0_int);
varnames=fieldnames(Terms_U0);


if pltcontour
for iplt=1:length(varnames)
figure(100+iplt)
clf
colormap(jet(20))
contourf(lat,lev,(getfield(Terms_U0,char(varnames(iplt)))'.*86400),[-10:0.2:10])
caxis([-10,10])
title(varnames(iplt))
set(gca,'yscale','log')
set(gca,'ydir','reverse')
xlabel('lat')
ylabel('lev')
xlim([-15,45])
ylim([100,1000]);
set(gca,'FontSize',18)
colorbar
end
else
figure(100)
clf
hold on
for iplt=1:length(varnames)
plot(lat,(getfield(TermsInt_U0,char(varnames(iplt))).*86400),'LineWidth',2)
end
legend(varnames)
xlabel('lat')
ylabel('response of U tendency m/s/day')
set(gca,'FontSize',18)
end

%% Terms_V0

TV0_fdU0=-repmat(beta*y,[1,length(lev)]).*dUbar;
TV0_fdU0_int=sum(TV0_fdU0(:,iptop:ipbot).*dp,2);
[~,TV0_vv]=gradient(-squeeze(mean((VV-V.*V)-(VV0-V0.*V0))),lev,y);
TV0_vv_int=sum(TV0_vv(:,iptop:ipbot).*dp,2);
TV0_dV0_tau=-dVbar./tau;
TV0_dV0_tau_int=sum(TV0_dV0_tau(:,iptop:ipbot).*dp,2);
dZbar=squeeze(mean(Z-Z0));
[~,TV0_dpgfbar]=gradient(dZbar,lev,y);
TV0_dpgfbar=-TV0_dpgfbar.*g;
TV0_dpgfbar_int=sum(TV0_dpgfbar(:,iptop:ipbot).*dp,2);
[~,Vy_temp]=gradient(Vbar0,lev,y);
TV0_dV0_V0=-dVbar.*Vy_temp;
TV0_dV0_V0_int=sum(TV0_dV0_V0(:,iptop:ipbot).*dp,2);
[~,Vy_temp]=gradient(dVbar,lev,y);
TV0_V0_dV0=-Vbar0.*Vy_temp;
TV0_V0_dV0_int=sum(TV0_V0_dV0(:,iptop:ipbot).*dp,2);

Terms_V0=struct('fdU0',TV0_fdU0,'vv',TV0_vv,'dV0_tau',TV0_dV0_tau,...
    'PGF',TV0_dpgfbar,'dV0_V0',TV0_dV0_V0,'V0_dV0',TV0_V0_dV0);
TermsInt_V0=struct('fdU0',TV0_fdU0_int,'vv',TV0_vv_int,'dV0_tau',TV0_dV0_tau_int,...
    'PGF',TV0_dpgfbar_int','dV0_V0',TV0_dV0_V0_int,'V0_dV0',TV0_V0_dV0_int);
varnames=fieldnames(Terms_V0);


if pltcontour
for iplt=1:length(varnames)
figure(200+iplt)
clf
colormap(jet(20))
contourf(lat,lev,(getfield(Terms_V0,char(varnames(iplt)))'.*86400),[-10:0.2:10])
caxis([-10,10])
title(varnames(iplt))
set(gca,'yscale','log')
set(gca,'ydir','reverse')
xlabel('lat')
ylabel('lev')
xlim([-15,45])
ylim([100,1000]);
set(gca,'FontSize',18)
colorbar
end
else
figure(200)
clf
hold on
for iplt=1:length(varnames)
plot(lat,(getfield(TermsInt_V0,char(varnames(iplt))).*86400),'LineWidth',2)
end
legend(varnames)
xlabel('lat')
ylabel('response of V tendency m/s/day')
set(gca,'FontSize',18)
end

%% Terms_H0 (Temperature)

[~,TH0_vt]=gradient(-squeeze(mean((VT-V.*T)-(VT0-T0.*V0))),lev,y);
TH0_vt_int=sum(TH0_vt(:,iptop:ipbot).*dp,2);
TH0_dT0_tau=-dTbar./tau;
TH0_dT0_tau_int=sum(TH0_dT0_tau(:,iptop:ipbot).*dp,2);
[~,Ty_temp]=gradient(Tbar0,lev,y);
TH0_dV0_T0=-dVbar.*Ty_temp;
TH0_dV0_T0_int=sum(TH0_dV0_T0(:,iptop:ipbot).*dp,2);
[~,Ty_temp]=gradient(dTbar,lev,y);
TH0_V0_dT0=-Vbar0.*Ty_temp;
TH0_V0_dT0_int=sum(TH0_V0_dT0(:,iptop:ipbot).*dp,2);
[Tp_temp,~]=gradient(Tbar0,lev.*100,y);
TH0_dOmega0_Tp0=-squeeze(mean(OMEGA-OMEGA0)).*Tp_temp;
TH0_dOmega0_Tp0_int=sum(TH0_dOmega0_Tp0(:,iptop:ipbot).*dp,2);

Terms_H0=struct('vt',TH0_vt,'dT0_tau',TH0_dT0_tau,...
    'dV0_T0',TH0_dV0_T0,'V0_dT0',TH0_V0_dT0,'dOmega0_Tp0',TH0_dOmega0_Tp0);
TermsInt_H0=struct('vt',TH0_vt_int,'dT0_tau',TH0_dT0_tau_int,...
    'dV0_T0',TH0_dV0_T0_int,'V0_dT0',TH0_V0_dT0_int,'dOmega0_Tp0',TH0_dOmega0_Tp0_int);
varnames=fieldnames(Terms_H0);


if pltcontour
for iplt=1:length(varnames)
figure(300+iplt)
clf
colormap(jet(20))
contourf(lat,lev,(getfield(Terms_H0,char(varnames(iplt)))'.*86400),[-10:0.2:10])
caxis([-10,10])
title(varnames(iplt))
set(gca,'yscale','log')
set(gca,'ydir','reverse')
xlabel('lat')
ylabel('lev')
xlim([-15,45])
ylim([100,1000]);
set(gca,'FontSize',18)
colorbar
end
else
figure(300)
clf
hold on
for iplt=1:length(varnames)
plot(lat,(getfield(TermsInt_H0,char(varnames(iplt))).*86400),'LineWidth',2)
end
legend(varnames)
xlabel('lat')
ylabel('response of T tendency K/day')
set(gca,'FontSize',18)
end


%% Terms_U1
fft_multi=(2/length(lon)).*repmat(exp(-1i*pi/180.*lon),[1,length(lat),length(lev)]);
dV1=squeeze(sum((V-V0).*fft_multi,1));
TU1_fdV1=repmat(beta*y,[1,length(lev)]).*dV1;
TU1_fdV1_int=sum(TU1_fdV1(:,iptop:ipbot).*dp,2);

duv1=squeeze(sum(((UV-U.*V)-(UV0-U0.*V0)).*fft_multi,1));
[~,TU1_uv]=gradient(-duv1,lev,y);
TU1_uv_int=sum(TU1_uv(:,iptop:ipbot).*dp,2);

duu1=squeeze(sum(((UU-U.*U)-(UU0-U0.*U0)).*fft_multi,1));
TU1_uu=-duu1.*1i*(1/Re);
TU1_uu_int=sum(TU1_uu(:,iptop:ipbot).*dp,2);

dU1=squeeze(sum((U-U0).*fft_multi,1));
TU1_dU1_tau=-dU1./tau;
TU1_dU1_tau_int=sum(TU1_dU1_tau(:,iptop:ipbot).*dp,2);

dZ1=squeeze(sum((Z-Z0).*fft_multi,1));
TU1_dpgf1=-1i.*dZ1.*g./Re;
TU1_dpgf1_int=sum(TU1_dpgf1(:,iptop:ipbot).*dp,2);

[~,Uy_temp]=gradient(squeeze(sum((U0).*fft_multi,1)),lev,y);
TU1_dV0_U1=-dVbar.*Uy_temp;
TU1_dV0_U1_int=sum(TU1_dV0_U1(:,iptop:ipbot).*dp,2);

[~,Uy_temp]=gradient(Ubar0,lev,y);
TU1_dV1_U0=-squeeze(sum((V-V0).*fft_multi,1)).*Uy_temp;
TU1_dV1_U0_int=sum(TU1_dV1_U0(:,iptop:ipbot).*dp,2);

[~,Uy_temp]=gradient(dU1,lev,y);
TU1_V0_dU1=-Vbar0.*Uy_temp;
TU1_V0_dU1_int=sum(TU1_V0_dU1(:,iptop:ipbot).*dp,2);

[~,Uy_temp]=gradient(dUbar,lev,y);
TU1_V1_dU0=-squeeze(sum(V0.*fft_multi,1)).*Uy_temp;
TU1_V1_dU0_int=sum(TU1_V1_dU0(:,iptop:ipbot).*dp,2);



Terms_U1=struct('fdV1',TU1_fdV1,'uv',TU1_uv,'dU1_tau',TU1_dU1_tau,'PGF',TU1_dpgf1,...
    'dV1_U0',TU1_dV1_U0,'V1_dU0',TU1_V1_dU0,'dV0_U1',TU1_dV0_U1,'V0_dU1',TU1_V0_dU1);
TermsInt_U1=struct('fdV1',TU1_fdV1_int,'uv',TU1_uv_int,'dU1_tau',TU1_dU1_tau_int,'PGF',TU1_dpgf1_int,...
    'dV1_U0',TU1_dV1_U0_int,'V1_dU0',TU1_V1_dU0_int,'dV0_U1',TU1_dV0_U1_int,'V0_dU1',TU1_V0_dU1_int);

varnames=fieldnames(Terms_U1);


if pltcontour
for iplt=1:length(varnames)
figure(400+iplt)
clf
colormap(jet(20))
contourf(lat,lev,real(getfield(Terms_U1,char(varnames(iplt)))'.*86400),[-10:0.2:10])
caxis([-10,10])
title(varnames(iplt))
set(gca,'yscale','log')
set(gca,'ydir','reverse')
xlabel('lat')
ylabel('lev')
xlim([-15,45])
ylim([100,1000]);
set(gca,'FontSize',18)
colorbar
end
else
figure(400)
clf
hold on
for iplt=1:length(varnames)
plot(lat,real(getfield(TermsInt_U1,char(varnames(iplt))).*86400),'LineWidth',2)
end
legend(varnames)
xlabel('lat')
ylabel('response of real(U1) tendency m/s/day')
set(gca,'FontSize',18)
figure(499)
clf
hold on
for iplt=1:length(varnames)
plot(lat,imag(getfield(TermsInt_U1,char(varnames(iplt))).*86400),'LineWidth',2)
end
legend(varnames)
xlabel('lat')
ylabel('response of imag(U1) tendency m/s/day')
set(gca,'FontSize',18)
end

%% Terms_V1
dU1=squeeze(sum((U-U0).*fft_multi,1));
TV1_fdU1=-repmat(beta*y,[1,length(lev)]).*dU1;
TV1_fdU1_int=sum(TV1_fdU1(:,iptop:ipbot).*dp,2);

dvv1=squeeze(sum(((VV-V.*V)-(VV0-V0.*V0)).*fft_multi,1));
[~,TV1_vv]=gradient(-dvv1,lev,y);
TV1_vv_int=sum(TV1_vv(:,iptop:ipbot).*dp,2);

dvu1=squeeze(sum(((UV-U.*V)-(UV0-U0.*V0)).*fft_multi,1));
TV1_vu=-dvu1.*1i*(1/Re);
TV1_vu_int=sum(TV1_vu(:,iptop:ipbot).*dp,2);

dV1=squeeze(sum((U-U0).*fft_multi,1));
TV1_dV1_tau=-dV1./tau;
TV1_dV1_tau_int=sum(TV1_dV1_tau(:,iptop:ipbot).*dp,2);

dZ1=squeeze(sum((Z-Z0).*fft_multi,1));
[~,TV1_dpgf1]=gradient(-dZ1*g,lev,y);
TV1_dpgf1_int=sum(TV1_dpgf1(:,iptop:ipbot).*dp,2);

[~,Vy_temp]=gradient(squeeze(sum((V0).*fft_multi,1)),lev,y);
TV1_dV0_V1=-dVbar.*Vy_temp;
TV1_dV0_V1_int=sum(TV1_dV0_V1(:,iptop:ipbot).*dp,2);

[~,Vy_temp]=gradient(Vbar0,lev,y);
TV1_dV1_V0=-squeeze(sum((V-V0).*fft_multi,1)).*Vy_temp;
TV1_dV1_V0_int=sum(TV1_dV1_V0(:,iptop:ipbot).*dp,2);

[~,Vy_temp]=gradient(dV1,lev,y);
TV1_V0_dV1=-Vbar0.*Vy_temp;
TV1_V0_dV1_int=sum(TV1_V0_dV1(:,iptop:ipbot).*dp,2);

[~,Vy_temp]=gradient(dVbar,lev,y);
TV1_V1_dV0=-squeeze(sum(V0.*fft_multi,1)).*Vy_temp;
TV1_V1_dV0_int=sum(TV1_V1_dV0(:,iptop:ipbot).*dp,2);



Terms_V1=struct('fdU1',TV1_fdU1,'vv',TV1_vv,'dV1_tau',TV1_dV1_tau,'PGF',TV1_dpgf1,...
    'dV1_V0',TV1_dV1_V0,'V1_dV0',TV1_V1_dV0,'dV0_V1',TV1_dV0_V1,'V0_dV1',TV1_V0_dV1);
TermsInt_V1=struct('fdU1',TV1_fdU1_int,'vv',TV1_vv_int,'dV1_tau',TV1_dV1_tau_int,'PGF',TV1_dpgf1_int,...
    'dV1_V0',TV1_dV1_V0_int,'V1_dV0',TV1_V1_dV0_int,'dV0_V1',TV1_dV0_V1_int,'V0_dV1',TV1_V0_dV1_int);

varnames=fieldnames(Terms_V1);


if pltcontour
for iplt=1:length(varnames)
figure(400+iplt)
clf
colormap(jet(20))
contourf(lat,lev,real(getfield(Terms_V1,char(varnames(iplt)))'.*86400),[-10:0.2:10])
caxis([-10,10])
title(varnames(iplt))
set(gca,'yscale','log')
set(gca,'ydir','reverse')
xlabel('lat')
ylabel('lev')
xlim([-15,45])
ylim([100,1000]);
set(gca,'FontSize',18)
colorbar
end
else
figure(500)
clf
hold on
for iplt=1:length(varnames)
plot(lat,real(getfield(TermsInt_V1,char(varnames(iplt))).*86400),'LineWidth',2)
end
legend(varnames)
xlabel('lat')
ylabel('response of real(V1) tendency m/s/day')
set(gca,'FontSize',18)
figure(599)
clf
hold on
for iplt=1:length(varnames)
plot(lat,imag(getfield(TermsInt_V1,char(varnames(iplt))).*86400),'LineWidth',2)
end
legend(varnames)
xlabel('lat')
ylabel('response of imag(V1) tendency m/s/day')
set(gca,'FontSize',18)
end

