clear
%close all

addpath /data/mgoes/matlab/Seawater/                    
addpath /data/mgoes/matlab/timeplt/                    
addpath /data/mgoes/AX8/Time_transform/                    
addpath /data/mgoes/matlab/                    

%Calculated in save_coast.jnl
ncload('/data/mgoes/SACZ/coast_32_15S.nc')
oi = find(ETOPO05_Y697_901>=-30);
ETOPO05_Y697_901 = ETOPO05_Y697_901(oi);
ZERO=ZERO(oi);



[dist,phaseangle] = sw_dist(ETOPO05_Y697_901,ZERO,'km');
figure,plot((ZERO(2:end)+ZERO(1:end-1))/2,phaseangle)


x = ZERO;                
y = ETOPO05_Y697_901;    
u = gradient(ZERO);                
v = gradient(ETOPO05_Y697_901);    
mag = sqrt(u.^2+v.^2);
u=u./mag;v=v./mag;

% [~,dx,dy] = sw_area(x,y);
% dx = diag(dx);dy = diag(dy);dy = dy(:);
% dRT = sqrt(dx.^2+dy.^2);
% dy_deg = gradient(y);
%    dy = 4.e7*dy_deg/360.;
% ATU = angle(dx+dy*sqrt(-1));
% Ul = U_VEL*0;
%  for i=1:nx
%   Ul(:,:,i) = -1*sin(angleU(:,i*ones(1,nz))).*squeeze(U_VEL(:,:,i))...
%    + cos(angleU(:,i*ones(1,nz))).*squeeze(V_VEL(:,:,i));
%  end
% end  %0


ncload('TAUX.nc')
ncload('TAUY.nc')
TAUX = taux;
TAUY = tauy;
clear taux tauy
miss = max(TAUX(:));
TAUX(TAUX==miss)=nan;
TAUY(TAUY==miss)=nan;
time = time + julian(1950,1,1);
month = julian2month(time);
yfrac = month2yfraction(month);
longitude(longitude>180)=longitude(longitude>180)-360;

%REDUCE AREA
loi = longitude>-60&longitude<20;
lai = latitude>-40&latitude<0;
TAUX = TAUX(:,lai,loi);
TAUY =TAUY(:,lai,loi);
longitude=longitude(loi);
latitude=latitude(lai);


nt = length(time);
nx = length(u);
[X,Y] = meshgrid(longitude,latitude);
taux = nan*ones(nt,nx);
tauy = nan*ones(nt,nx);
tauxu = nan*ones(nt,nx);
tauyv = nan*ones(nt,nx);
tau   = nan*ones(nt,nx);
for ii = 1:nt
   taux1 = squeeze(TAUX(ii,:,:)); 
   tauy1 = squeeze(TAUY(ii,:,:)); 
   oi = ~isnan(taux1);
   F = scatteredInterpolant(X(oi),Y(oi),taux1(oi));
   tauxii = F(x,y);
   taux(ii,:,:) = tauxii;%.*u;
   F = scatteredInterpolant(X(oi),Y(oi),tauy1(oi));
   tauyii = F(x,y);
   tauy(ii,:,:) = tauyii;%.*v;   

  tauxu(ii,:,:) = tauxii.*u;
  tauyv(ii,:,:) = tauyii.*v;   
  %tauxu(ii,:,:) = cos(ATU).*tauxii - sin(ATU).*tauyii;
  %tauyv(ii,:,:) = sin(ATU).*tauxii + cos(ATU).*tauyii;

  
%   tauxu(ii,:,:) = tauxii.*u-tauyii.*v;
%   tauyv(ii,:,:) = tauyii.*v+tauxii.*u;   
  tau(ii,:) = tauxii.*u + tauyii.*v;
end


f = sw_f(y)';f = f(ones(1,nt),:);
rho = 1025;
upwelling = tau./(rho*f);


upwelling = detrend(upwelling);

clim = make_clim(upwelling,month,'clim');
clim = squeeze(clim);
anom = upwelling - make_clim(upwelling,month);

lai = y>-23&y<-18;
CB  = nanmean(upwelling(:,lai),2);
CBa = nanmean(anom(:,lai),2);

if 0
ii=10;
figure,
quiver(X,Y,squeeze(TAUX(ii,:,:)),squeeze(TAUY(ii,:,:)),2)
hold on
quiver(x,y,squeeze(tauxu(ii,:))',squeeze(tauyv(ii,:))',2)
%quiver(x,y,squeeze(taux(ii,:))',squeeze(tauy(ii,:))',2)
quiver(x,y,u,v,2,'k')
axis([-55 -35 -34 -14])


figure,plot(yfrac,CB)

figure,pcolor(yfrac,y,anom')
shading interp
caxis([-.8 .8])
xlim([1993 2011])
colormap(mapa)

figure,plot(yfrac,weim(3,'rect',CBa))


load /data/mgoes/AVISO/SA_Atl_Sep8_2015/EOF/mapa.mat
figure,contourf(1:12,y,clim',-1:.1:1);
caxis([-.8 .8])
colormap(mapa)
end
stop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load(['/data/mgoes/SA_gyre/AX97/Data/transport_from_altimetry_nomess.mat']) 
 bc    = [tsp_bc2.transp_400_minus]/1e6;
[bcmo,~,ymo] = make_monthly_mean(bc,year_mon); 
bcmoc        = make_clim(bcmo',ymo);
bcmoa        = bcmo - bcmoc';
 timey = month2yfraction(ymo);%month2yfraction( julian2month(time+14)); %center in 15 days instead of 1
 save('upwelling_BC_oaflux.mat','bcmoa','CBa','timey','ymo','yfrac','month')

%mywidefigure(4)
%plot(timey,weim(3,'rect',bcmoa),'k')
%hold on
%plot(yfrac,-weim(5,'rect',CBa)*10,'r')
%xlim([1993 2015])
%ylim([-3 3])


mywidefigure(4)
%plot(timey,weim(3,'rect',bcmoa),'k')
plot(timey,bcmoa,'k')
hold on
plot(yfrac,-weim(3,'rect',CBa)*10,'r')
xlim([1993 2015])
ylim([-4 4])
legend('BC transport','Upwelling')
ylabel('Anomalies')


load /data/mgoes/SA_gyre/Wind/Residual_22S_5mo.mat

mywidefigure(5)
plot(timeyR,residual,'k')
hold on
plot(yfrac,-weim(3,'rect',CBa)*12,'b')
xlim([1993 2015])
ylim([-5 5])
legend('BC transport','Upwelling')
ylabel('Anomalies')

clear
load upwelling_BC_oaflux.mat
CBo = CBa;
yfraco=yfrac;
load upwelling_BC_era.mat
io = yfraco>2012;
ie = yfrac>2012;
yfrac = [yfraco(~io) yfrac(ie)];
CB    = [CBo(~io); CBa(ie)];
mywidefigure(4)
plot(timey,bcmoa,'k')
hold on
plot(yfrac,-weim(3,'rect',CB)*10,'r')
xlim([1993 2015])
ylim([-4 4])
legend('BC transport','Upwelling')
ylabel('Anomalies')
grid on
set(gca,'fontsize',16)
print -djpeg -r300 upwelling_era_tropflux.jpg
