%% E vector DIV ~ dU
% need the output from drift_Ubudget.m, 'ENS_Ubudget.mat', which is 
% calculated on cluster and downloaded. 
levshowtop=100.;
levshowbot=600.;
scale=4;

MJONAME='MJOC40';
FRCNAME='EQ1M40lon180300';
offsetplt=15;
Day0=1;
Day1=20;
vinthname='vinth';
leftrightplt=2;


%  MJONAME='MJOB';
%  FRCNAME='EQ1M40lon60180';
%  offsetplt=0;
%  Day0=1;
%  Day1=20;
%  vinthname='vinth3';
% leftrightplt=1;
%vinth3 in CTRL, MJOB, MJOC are logp interpolated, so is vinth in MJOC40.

file0='vinth_All_SP1DJF35_FREE.nc';
file1=['vinth_All_SP1DJF35_' FRCNAME '_10.nc'];
file_trans0=file0;
file_trans1=file1;


% Day0=18;
% Day1=20;
% file0='vinth_timeavg_CTRL_ENSavg.nc';
% file1='vinth_lasttwomean_MJOC_reynold_term_ENSavg.nc';
% file_trans0='vinth_CTRL_H2_ENS800.nc';
% file_trans1='vinth_MJOC_H2_ENS800.nc';

%% calculate E vector
L=Day1-Day0+1;
L0=20;
T1day=86400.;
d2r=pi/180.;
a=6371000.;
f0=4*pi/86400;
g=9.8;
tau=40*86400;
kappa=2./7.;

% prepare coordinate
lat=squeeze(double(ncread(file0,'lat')));
lon=squeeze(double(ncread(file0,'lon')));
lev=squeeze(double(ncread(file0,'lev')));
ilev=squeeze(double(ncread(file0,'ilev')));
nlon=length(lon);
nlat=length(lat);
nlev=length(lev);
ip1=[2:nlev,nlev];
ip0=[1,1:nlev-1];
pres=permute(repmat(lev.*100,[1,length(lon),length(lat)]),[2,3,1]);
dp=ilev(2:end)-ilev(1:end-1);
iy1=[2:nlat,nlat];
iy0=[1,1:nlat-1];
dlat=(lat(2)-lat(1))*d2r*2;
ix1=[2:nlon,1];
ix0=[nlon,1:nlon-1]; 
dlon=(lon(2)-lon(1))*d2r*2;
coslat_=repmat(cos(lat.*d2r)',[length(lon),1,length(lev)]);
sinlat_=repmat(sin(lat.*d2r)',[length(lon),1,length(lev)]);
seclat=1./cos(lat.*d2r)';
seclat(1)=0;
seclat(end)=0;
seclat_=repmat(seclat,[nlon,1,nlev]);

%read data
uu=squeeze(double(ncread(file1,'UU')))-squeeze(double(ncread(file1,'U'))).^2 ...
    -(squeeze(double(ncread(file0,'UU')))-squeeze(double(ncread(file0,'U'))).^2);
vu=squeeze(double(ncread(file1,'VU')))-squeeze(double(ncread(file1,'V'))).*squeeze(double(ncread(file1,'U')))...
    -(squeeze(double(ncread(file0,'VU')))-squeeze(double(ncread(file0,'V'))).*squeeze(double(ncread(file0,'U'))));
vv=squeeze(double(ncread(file1,'VV')))-squeeze(double(ncread(file1,'V'))).^2 ...
    -(squeeze(double(ncread(file0,'VV')))-squeeze(double(ncread(file0,'V'))).^2);
vt=squeeze(double(ncread(file1,'VT')))-squeeze(double(ncread(file1,'T'))).*squeeze(double(ncread(file1,'V')))...
    -(squeeze(double(ncread(file0,'VT')))-squeeze(double(ncread(file0,'T'))).*squeeze(double(ncread(file0,'V'))));
T_CTRL=squeeze(double(ncread(file0,'T')));
U_CTRL=squeeze(double(ncread(file0,'U')));
V_CTRL=squeeze(double(ncread(file0,'V')));
Z_CTRL=squeeze(double(ncread(file0,'Z3')));

T_MJO=squeeze(double(ncread(file1,'T')));
U_MJO=squeeze(double(ncread(file1,'U')));
V_MJO=squeeze(double(ncread(file1,'V')));
Z_MJO=squeeze(double(ncread(file1,'Z3')));

dU=U_MJO-U_CTRL;
dV=V_MJO-V_CTRL;
dZ=smoothn(Z_MJO-Z_CTRL,'robust');

Tbar=repmat(mean(T_MJO),[nlon,1,1]);
TH_p=(Tbar(:,:,ip1)-Tbar(:,:,ip0))./(pres(:,:,ip1)-pres(:,:,ip0))+Tbar*(-kappa)./pres;

%EP flux and divergence
%EPx=(vv-uu);
EPx=-uu;
EPy=-vu;
EPz=f0.*coslat_.*vt./TH_p;

DIVEPX=1/a.*seclat_.*(EPx(ix1,:,:)-EPx(ix0,:,:))./(dlon);
DIVEPY=1/a.*seclat_.*(EPy(:,iy1,:).*coslat_(:,iy1,:)-EPy(:,iy0,:).*coslat_(:,iy0,:))./(dlat);
DIVEP=DIVEPX...
      + DIVEPY;...
      - (EPz(:,:,ip1)-EPz(:,:,ip0))./(pres(:,:,ip1)-pres(:,:,ip0));
  
%advection term of shifting terms
Adv_Ux_climadvdrift=U_CTRL.*1/a.*seclat_.*(dU(ix1,:,:)-dU(ix0,:,:))./(dlon).*T1day;
Adv_Ux_driftadvclim=dU.*1/a.*seclat_.*(U_CTRL(ix1,:,:)-U_CTRL(ix0,:,:))./(dlon).*T1day;
Adv_Ux_driftadvdrift=dU.*1/a.*seclat_.*(dU(ix1,:,:)-dU(ix0,:,:))./(dlon).*T1day;
Adv_Ux=Adv_Ux_climadvdrift+Adv_Ux_driftadvclim+Adv_Ux_driftadvdrift;

Adv_Uy_climadvdrift=V_CTRL.*1/a.*(dU(:,iy1,:)-dU(:,iy0,:))./(dlat).*T1day;
Adv_Uy_driftadvclim=dV.*1/a.*(U_CTRL(:,iy1,:)-U_CTRL(:,iy0,:))./(dlat).*T1day;
Adv_Uy_driftadvdrift=dV.*1/a.*(dU(:,iy1,:)-dU(:,iy0,:))./(dlat).*T1day;
Adv_Uy=Adv_Uy_climadvdrift+Adv_Uy_driftadvclim+Adv_Uy_driftadvdrift;

%Adv_Ux=0 ;...
%    +(1/a.*seclat_.*(dU(ix1,:,:).*U_CTRL(ix1,:,:)-dU(ix0,:,:).*U_CTRL(ix0,:,:))./(dlon)+1/a.*seclat_.*(U_CTRL(ix1,:,:).*dU(ix1,:,:)-U_CTRL(ix0,:,:).*dU(ix0,:,:))./(dlon)).*T1day; ...
%    +dU.*1/a.*seclat_.*(dU(ix1,:,:).*dU(ix1,:,:)-dU(ix0,:,:).*dU(ix0,:,:))./(dlon).*T1day;
%Adv_Uy=0 ;...
%    +(1/a.*(dU(:,iy1,:).*V_CTRL(ix1,:,:)-dU(:,iy0,:).*V_CTRL(ix0,:,:))./(dlat)+1/a.*(U_CTRL(:,iy1,:).*dV(ix1,:,:)-U_CTRL(:,iy0,:).*dV(ix0,:,:))./(dlat)).*T1day; ...
%    +1/a.*(dU(:,iy1,:).*dV(ix1,:,:)-dU(:,iy0,:).*dV(ix0,:,:))./(dlat).*T1day;

%ageostropic v
Cori_ageo=(-f0.*sinlat_.*dV+g.*1/a.*seclat_.*smoothn((dZ(ix1,:,:)-dZ(ix0,:,:)),'robust')./(dlon)).*T1day;


% \Delta U
%U1=squeeze(double(ncread(file_trans1,'U')));
%U1=squeeze(U1(:,:,:,Day1));
if Day0==1
    %U0=squeeze(double(ncread(file_trans0,'U')));
    %U0=squeeze(U0(:,:,:,L));
    dU_h2=dU;
    %dU_h2=U1-U0;
else
    %U0=squeeze(double(ncread(file_trans1,'U')));
    %U0=squeeze(U0(:,:,:,Day0));
    dU_h2=dU;
    %dU_h2=U1-U0;
end
dU_h2_=dU_h2./L;
restore_amp=lev.*0+1;
restore_amp(end-5:end)=restore_amp(end-5:end).*10;
Restore_U=dU./40.*permute(repmat(restore_amp,[1,nlon,nlat]),[2 3 1]);
dU_h2_=dU_h2_+Restore_U;

dU_h2=dU_h2_+Adv_Ux+Adv_Uy+Cori_ageo;


% eddy energy gain
Ux=1/a.*seclat_.*(U_CTRL(ix1,:,:)-U_CTRL(ix0,:,:))./(dlon);
Uy=1/a.*(U_CTRL(:,iy1,:)-U_CTRL(:,iy0,:))./(dlat);
EddyEnergyGain=EPx.*(Ux)+EPy.*(Uy);
  
%vertical integrate
struct_longterm=[];
dp_=permute(repmat(dp,[1,length(lon),length(lat)]),[2,3,1]);
indlev=find(lev<levshowbot&lev>levshowtop);
struct_longterm.du_h2=squeeze(nansum(dU_h2(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
struct_longterm.du_h2_=squeeze(nansum(dU_h2_(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
struct_longterm.divep=squeeze(nansum(DIVEP(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev)).*T1day;
struct_longterm.divepx=squeeze(nansum(DIVEPX(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev)).*T1day;
struct_longterm.divepy=squeeze(nansum(DIVEPY(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev)).*T1day;
struct_longterm.eddyenergygain=squeeze(nansum(EddyEnergyGain(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev)).*T1day./10;
struct_longterm.ex=squeeze(nansum(EPx(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev)).*T1day/1e6;
struct_longterm.ey=squeeze(nansum(EPy(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev)).*T1day/1e6;
struct_longterm.adv_ux=squeeze(nansum(Adv_Ux(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
struct_longterm.adv_uy=squeeze(nansum(Adv_Uy(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
struct_longterm.cori_ageo=squeeze(nansum(Cori_ageo(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
% fft_coria=fft(smoothn(struct_longterm.cori_ageo,'robust'),[],2);
% fft_coria(:,6:end)=0.;
% struct_longterm.cori_ageo=ifft(fft_coria,[],2,'symmetric');
struct_longterm.adv_ux_clim_drift=squeeze(nansum(Adv_Ux_climadvdrift(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
struct_longterm.adv_uy_clim_drift=squeeze(nansum(Adv_Uy_climadvdrift(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
struct_longterm.adv_ux_drift_clim=squeeze(nansum(Adv_Ux_driftadvclim(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
struct_longterm.adv_uy_drift_clim=squeeze(nansum(Adv_Uy_driftadvclim(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
struct_longterm.adv_ux_drift_drift=squeeze(nansum(Adv_Ux_driftadvdrift(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
struct_longterm.adv_uy_drift_drift=squeeze(nansum(Adv_Uy_driftadvdrift(:,:,indlev).*dp_(:,:,indlev),3))'./sum(dp(indlev));
struct_longterm.adv=struct_longterm.adv_ux+struct_longterm.adv_uy;
struct_longterm.adv_drift_drift=struct_longterm.adv_ux_drift_drift+struct_longterm.adv_uy_drift_drift;
struct_longterm.adv_drift_clim=struct_longterm.adv_ux_drift_clim+struct_longterm.adv_uy_drift_clim;
struct_longterm.adv_clim_drift=struct_longterm.adv_ux_clim_drift+struct_longterm.adv_uy_clim_drift;
struct_longterm.du_h2_other=struct_longterm.du_h2-struct_longterm.du_h2_;
struct_longterm.all_drift_tendency=struct_longterm.divep-struct_longterm.adv-struct_longterm.cori_ageo;
struct_longterm.advx_divepx=struct_longterm.adv_ux+struct_longterm.divepx;
struct_longterm.advy_divepy=struct_longterm.adv_uy+struct_longterm.divepy;


%% plot \Delta U and \Delta div(E)

figure(1+offsetplt)
clf
colormap(polarmap(21))
contourf(lon,lat,smoothn(struct_longterm.du_h2,'robust'),[-scale-10,-scale+scale/10:scale/10.:scale],'LineColor','none');
caxis([-scale,scale])

hold on
load coastlines
coastlon_=coastlon;
coastlon_(coastlon_<0)=coastlon_(coastlon_<0)+360;
plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
%colormap(jet(21))
fft_divep=fft(smoothn(struct_longterm.divep,'robust'),[],2);
fft_divep(:,end:end)=0.;
trunc_divep=ifft(fft_divep,[],2,'symmetric');
contour(lon,lat,trunc_divep,[scale/10:scale/5.:scale+10],'LineWidth',1,'Color','k','LineStyle','-');
contour(lon,lat,trunc_divep,-[scale/10.:scale/5.:scale+10],'LineWidth',1,'Color','k','LineStyle','--');
contour(lon,lat,trunc_divep,[0.],'LineWidth',1,'Color','k','LineStyle','-.');
%caxis([-scale,scale])

xlabel('lon')
ylabel('lat')
set(gca,'FontSize',18)
%title(['DIV(E) (forced-unforced)  ' MJONAME ' ' sprintf('%d',levshowbot) '-100mb'])
title(MJONAME)
set(gcf,'Position',[0 0 842 369])
yticks([-90:30:90])
xlim([0,360])
xticks([0:45:360])

colorbar


figure(6+offsetplt)
clf
colormap(polarmap(21))
contourf(lon,lat,smoothn(struct_longterm.du_h2_,'robust'),[-scale/2-10,-scale/2+scale/2/10:scale/2/10.:scale/2],'LineColor','none');
caxis([-scale/2,scale/2])

hold on
load coastlines
coastlon_=coastlon;
coastlon_(coastlon_<0)=coastlon_(coastlon_<0)+360;
plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
%colormap(jet(21))
fft_divep=fft(smoothn(struct_longterm.all_drift_tendency,'robust'),[],2);
fft_divep(:,end:end)=0.;
trunc_divep=ifft(fft_divep,[],2,'symmetric');
contour(lon,lat,trunc_divep,[scale/10:scale/2/5.:scale/2+10],'LineWidth',1,'Color','k','LineStyle','-');
contour(lon,lat,trunc_divep,-[scale/10.:scale/2/5.:scale/2+10],'LineWidth',1,'Color','k','LineStyle','--');
%contour(lon,lat,trunc_divep,[0.],'LineWidth',1,'Color','k','LineStyle','-.');
%caxis([-scale,scale])

xlabel('lon')
ylabel('lat')
set(gca,'FontSize',18)
%title(['DIV(E) (forced-unforced)  ' MJONAME ' ' sprintf('%d',levshowbot) '-100mb'])
title(MJONAME)
set(gcf,'Position',[0 0 842 369])
yticks([-90:30:90])
xlim([0,360])
xticks([0:45:360])
colorbar

%% all terms
%pltvars={'du_h2_','adv','cori_ageo','divep','du_h2_other','adv_clim_drift','adv_drift_clim','adv_drift_drift'};
%titles={'$\langle u(20)\rangle-\langle u(0)\rangle$','Advective terms','Ageo Coriolis force','Reynold Drag','Advection+ageo Coriolis','$u_0\hat{u}_x+v_0\hat{u}_y$','$\hat{u}u_{0x}+\hat{v}u_{0y}$','$\hat{u}\hat{u}_x+\hat{v}\hat{u}_y$'};

pltvars={'du_h2_','adv','cori_ageo','divep','du_h2_other','adv_clim_drift','adv_drift_clim','all_drift_tendency'};
titles={'$\langle u(20)\rangle-\langle u(0)\rangle$','Advective terms','Ageo Coriolis force','Reynold Drag','Advection+ageo Coriolis','$u_0\hat{u}_x+v_0\hat{u}_y$','$\hat{u}u_{0x}+\hat{v}u_{0y}$','all drifting terms'};


figure(2+offsetplt)
clf
for iplt=1:length(pltvars)
subplot(4,2,iplt)
colormap(polarmap(21))
contourf(lon,lat,smoothn(getfield(struct_longterm,char(pltvars(iplt))),'robust'),[-scale-10,-scale+scale/10:scale/10.:scale],'LineColor','none')
caxis([-scale,scale])
xlabel('lon')
ylabel('lat')
set(gca,'FontSize',18)
hold on
plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
title(char(titles(iplt)),'Interpreter','Latex')
%set(gcf,'Position',[0 0 842 369])
yticks([-90:30:90])
xlim([0,360])
xticks([0:45:360])
colorbar
end

%% all terms
%pltvars={'du_h2_','adv','cori_ageo','divep','du_h2_other','adv_clim_drift','adv_drift_clim','adv_drift_drift'};
%titles={'$\langle u(20)\rangle-\langle u(0)\rangle$','Advective terms','Ageo Coriolis force','Reynold Drag','Advection+ageo Coriolis','$u_0\hat{u}_x+v_0\hat{u}_y$','$\hat{u}u_{0x}+\hat{v}u_{0y}$','$\hat{u}\hat{u}_x+\hat{v}\hat{u}_y$'};

pltvars={'du_h2_','adv','cori_ageo','divep','advx_divepx','advy_divepy'};
titles={'$\langle u(20)\rangle-\langle u(0)\rangle$','Advective terms','Ageo Coriolis force','Reynold Drag','Adv_x+Renold_x','Adv_y+Renold_y'};
nvar=length(pltvars)-2;

figure(4)
if leftrightplt==1
    clf
end
for iplt=1:nvar
hold on

% if leftrightplt==2
%     hp4=get(subplot(nvar,2,(iplt-1)*2+leftrightplt),'Position');
%     colorbar('Position', [hp4(1)+hp4(3)+0.01  hp4(2)  0.01  hp4(4)])
% end

subplot(nvar,2,(iplt-1)*2+leftrightplt)
(iplt-1)*2+leftrightplt
subaxis(nvar,2,(iplt-1)*2+leftrightplt,'sh', 0.01, 'sv', 0.02, 'padding', 0.00, 'margin', 0.08)
colormap(polarmap(21))
contourf(lon,lat,smoothn(getfield(struct_longterm,char(pltvars(iplt))),'robust'),[-scale-10,-scale+scale/10:scale/10.:scale],'LineColor','none')
caxis([-scale,scale])
yc=[-90:30:90];
xc=[0:45:360];
yticks(yc)
xlim([0,360])
xticks(xc)

if iplt==nvar
    xlabel('lon')
    set(gca,'XTickLabel',string(xc));
else
    set(gca,'XTickLabel',strings(length(xc),1));
end
if leftrightplt==1
    ylabel('lat')
    set(gca,'YTickLabel',string(yc));
else
    set(gca,'YTickLabel',strings(length(xc),1));
end
set(gca,'FontSize',18)
hold on
plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
%title(char(titles(iplt)),'Interpreter','Latex')

end

%hp4=get(subplot(nvar,2,(iplt-1)*2+leftrightplt),'Position');
hp4=[0.5703    0.1100    0.3347    0.1026];
colorbar('Position', [hp4(1)+hp4(3)+0.04  hp4(2)+ hp4(4) 0.01  hp4(4)*5])

% 
% figure(3)
% colormap(polarmap(21))
% contourf(lon,lat,cori_ageo,[-scale-10,-scale+scale/10:scale/10.:scale],'LineColor','none')
% caxis([-scale,scale])
% xlabel('lon')
% ylabel('lat')
% set(gca,'FontSize',18)
% hold on
% plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
% title('Coriolis force due to ageostrophic v')
% set(gcf,'Position',[0 0 842 369])
% yticks([-90:30:90])
% xlim([0,360])
% xticks([0:45:360])
% colorbar
% 
% figure(4)
% %clf
% colormap(polarmap(21))
% contourf(lon,lat,adv_uy+adv_ux,[-scale-10,-scale+scale/10:scale/10.:scale],'LineColor','none')
% caxis([-scale,scale])
% xlabel('lon')
% ylabel('lat')
% set(gca,'FontSize',18)
% hold on
% plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
% title('advection terms')
% set(gcf,'Position',[0 0 842 369])
% yticks([-90:30:90])
% xlim([0,360])
% xticks([0:45:360])
% colorbar
% 
% figure(5)
% %clf
% colormap(polarmap(21))
% contourf(lon,lat,trunc_divep,[-scale-10,-scale+scale/10:scale/10.:scale],'LineColor','none')
% caxis([-scale,scale])
% xlabel('lon')
% ylabel('lat')
% set(gca,'FontSize',18)
% hold on
% plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
% title('DIV(E)')
% set(gcf,'Position',[0 0 842 369])
% yticks([-90:30:90])
% xlim([0,360])
% xticks([0:45:360])
% colorbar
% 
figure(7+offsetplt)
clf
colormap(polarmap(21))
contourf(lon,lat,struct_longterm.du_h2_,[-scale-10,-scale+scale/10:scale/10.:scale],'LineColor','none')
caxis([-scale,scale])
xlabel('lon')
ylabel('lat')
set(gca,'FontSize',18)
hold on
plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
title('$\frac{\partial u}{\partial t}$','Interpreter','Latex')
set(gcf,'Position',[0 0 842 369])
yticks([-90:30:90])
xlim([0,360])
xticks([0:45:360])
colorbar

%% advective and Reynold about udx
figure(8+offsetplt)
clf
colormap(polarmap(21))
contourf(lon,lat,-struct_longterm.adv_ux+struct_longterm.divepx,[-scale-10,-scale+scale/10:scale/10.:scale],'LineColor','none')
caxis([-scale,scale])
xlabel('lon')
ylabel('lat')
set(gca,'FontSize',18)
hold on
plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
title('adv_x+reynold_x')
set(gcf,'Position',[0 0 842 369])
yticks([-90:30:90])
xlim([0,360])
xticks([0:45:360])
colorbar

figure(9+offsetplt)
clf
colormap(polarmap(21))
contourf(lon,lat,-struct_longterm.adv_uy+struct_longterm.divepy,[-scale-10,-scale+scale/10:scale/10.:scale],'LineColor','none')
caxis([-scale,scale])
xlabel('lon')
ylabel('lat')
set(gca,'FontSize',18)
hold on
plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
title('adv_y+reynold_y')
set(gcf,'Position',[0 0 842 369])
yticks([-90:30:90])
xlim([0,360])
xticks([0:45:360])
colorbar


%% plot \Delta U and \Delta eddy energy gain
% figure(22+offsetplt)
% clf
% % colormap(polarmap(21))
% % contourf(lon,lat,du_h2,[-scale-10:scale/10.:scale],'LineColor','none');
% % caxis([-scale,scale])
% % 
% % hold on
% %colormap(jet(21))
% contour(lon,lat,eddyenergygain,[-scale-10,-scale+scale/10:scale/10.:scale],'LineWidth',2);
% caxis([-scale,scale])
% title('d eddy energy gain (forced-unforced, forced ref)','Interpreter','none')
% xlabel('lon')
% ylabel('lat')
% set(gca,'FontSize',18)
% 
% 
% hold on
% plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
% set(gcf,'Position',[0 0 842 369])
% colorbar

%% plot  \Delta U and Ex
% figure(33+offsetplt)
% clf
% % colormap(polarmap(21))
% % contourf(lon,lat,du_h2,[-scale-10:scale/10.:scale],'LineColor','none');
% % caxis([-scale,scale])
% % 
% % hold on
% %colormap(jet(21))
% contour(lon,lat,struct.ex,[-scale-10:scale/10.:scale],'LineWidth',2);
% caxis([-scale,scale])
% title('d Ex (forced-unforced, forced ref)','Interpreter','none')
% xlabel('lon')
% ylabel('lat')
% set(gca,'FontSize',18)
% 
% hold on
% plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
% set(gcf,'Position',[0 0 842 369])
% yticks([-90:30:90])
% xlim([0,360])
% xticks([0:45:360])
% colorbar
% 
% % plot  \Delta U and Ey
% figure(44+offsetplt)
% clf
% % colormap(polarmap(21))
% % contourf(lon,lat,du_h2,[-scale-10:scale/10.:scale],'LineColor','none');
% % caxis([-scale,scale])
% % 
% % hold on
% %colormap(jet(21))
% contour(lon,lat,struct.ey,[-scale-10:scale/10.:scale],'LineWidth',2);
% caxis([-scale,scale])
% title('d Ey (forced-unforced, forced ref)','Interpreter','none')
% xlabel('lon')
% ylabel('lat')
% set(gca,'FontSize',18)
% 
% hold on
% plot(coastlon_,coastlat,'.','MarkerSize',0.3,'Color',[0.5,0.5,0.5])
% set(gcf,'Position',[0 0 842 369])
% yticks([-90:30:90])
% xlim([0,360])
% xticks([0:45:360])
% colorbar