% estimate sss gradients using Dong (2006) methodology
%
% Comente o programa

  
load sssCAP_2014.mat
sss=[sss(:,181:360) sss(:,1:180)];  
lon=lon-180;

gap=[.08,.08];
marg_h=[.08,.06];
marg_v=[.06,.04];
% Interpolation

ilat=[lat(1):2.5:lat(end)];
ilon=[lon(1):2.5:lon(end)]';
sssi=griddata(lon,lat,sss,ilon,ilat);

% -----------------------------------------------------------
% Distância 
% Em que unidades estão as distâncias?

dx=111*ones(size(ilon));
dy=111*cosd(ilat);
[dy,dx]=meshgrid(dx,dy);

% -----------------------------------------------------------
% Estimate the absolute gradient

aux=sssi;
aux1=smoo(aux,nan,5,1);
aux=a2b(aux1,-9999,0);
[dsalx,dsaly]=gradient(aux);
% get rid off high values near the coastlines
L=abs(dsalx)>1;
dsalx(L)=nan;
L=abs(dsaly)>1;
dsaly(L)=nan;
dsalx=a2b(dsalx,0,nan);
dsaly=a2b(dsaly,0,nan);
dSx=dsalx./dx;
dSy=dsaly./dy;
dsal=sqrt(dSx.*dSx+dSy.*dSy);
clear L aux* i

% -----------------------------------------------------------
% Make plot

subtightplot(1,1,1,gap,marg_h,marg_v); 
imagesc(ilon,ilat,dsal,[0 .0025]);axis('xy');c=colorbar;
cc=get(c,'title');
set(cc,'string','dS/km')
t=title(['Gradiente de SSS']);
set(t,'fontsize',10);
xlabel('longitude')
ylabel('latitude')
set(gca,'fontsize',8)

print -dpng dsss