clear
close all

addpath /data/mgoes/matlab/Seawater/                    

%Calculated in save_coast.jnl
ncload('/data/mgoes/SACZ/coast_32_15S.nc')




[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));
ATU = angle(dx-dy*sqrt(-1));%-180;
% 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,:,:)); 
   F = scatteredInterpolant(X(:),Y(:),taux1(:));
   tauxii = F(x,y);
   taux(ii,:,:) = tauxii;%.*u;
   F = scatteredInterpolant(X(:),Y(:),tauy1(:));
   tauyii = F(x,y);
   tauy(ii,:,:) = tauyii;%.*v;   
  %    cost = 
  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

ii=10;
figure,
quiver(X,Y,squeeze(TAUX(ii,:,:)),squeeze(TAUY(ii,:,:)),4)
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])

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



