% extract drogued and slip-corrected undrogued drifters along 23W and 10W
% and slip-corrected Argo floats

if (0) %~exist('v')
    
    % load Argo floats
    load floats lon lat t id u v drogue Wx Wy 
    u2=u; v2=v; id2=id; t2=t; Wx2=Wx; Wy2=Wy;
    lon2=lon; lat2=lat;
    drogue2=drogue;
    
    % load drifters
    load drifters
    
    u=u_orig; v=v_orig; clear *orig;
    
    % discard drifters with NaN velocity or drogue status uncertain
    q=find(isnan(u+v) | drogue==.5); 
    truncatedrift;
    
    % combine drifter and float observations
    u=[u;u2]; v=[v;v2]; id=[id;id2]; t=[t;t2];
    drogue=[drogue;drogue2];
    lon=[lon;lon2]; lat=[lat;lat2];
    Wx=[Wx;Wx2]; Wy=[Wy;Wy2];
    clear u2 v2 id2 t2 drogue2 lon2 lat2 Wx2 Wy2
    
    
    % remove the slip. 
    %% Drogued: 7e-4 m/s per m/s of wind. (Niiler & Paduan, 1995)
    %% Undrogued: 7e-4+7.9e-3=8.6e-3 (", Pazan & Niiler 2000)  
    q=find(drogue==1);
    u(q)=u(q)-7e-4*Wx(q);
    v(q)=v(q)-7e-4*Wy(q); 
    
    q=find(drogue==0);  
    bestA=1.64e-2 + 7e-4;   % value from Perez et al. (2013)
    u(q)=u(q)-bestA*Wx(q);
    v(q)=v(q)-bestA*Wy(q);
    
    q=find(drogue==-1);  
    bestA=1.14e-2 + 7e-4;   % value from Perez et al. (2013)
    u(q)=u(q)-bestA*Wx(q);
    v(q)=v(q)-bestA*Wy(q);
    
     
    save drifters_floats_eq lon lat t id u v drogue
    
else
    load drifters_floats_eq lon lat t id u v drogue
end;

% calculate mean U and V for box with edges at 23W, 10W and 0E; 2S to 2N
Lon=[[-23-2 -23+2]; [-23 (-23-10)/2]; [(-23 -10)/2 -10]; ...
        [-10-2 -10+2];  [-10 -5]; [-5 0]; ...
        [-2 2] ];  % longitude bands for average
dLat=1;
Lat=-2:dLat/2:2; 

% store mean velocity values. 
%  U( Lat x Lon band x 3).  
%   Third index: 1=data, 2=std. dev., 3=num obs (days)
U=NaN*ones(length(Lat),size(Lon,1),3);
V=U;
FracNoDrogue=NaN*ones(length(Lat),size(Lon,1));
Ud=U; Vd=V; % mean speed of drogued drifters
Uu=U; Vu=V; % mean speed of undrogued, corrected drifters

% also save monthly values of U and V
tmon=(.5/12:1/12:11.5/12)';
Umon=NaN*ones(length(Lat),size(Lon,1),12);
Vmon=Umon;
model_mon=[1+0*tmon, sin(tmon*2*pi),cos(tmon*2*pi), ...
                sin(tmon*4*pi),cos(tmon*4*pi)];  
            
NumObs=NaN*ones(length(Lat),size(Lon,1),3);

for i=1:length(Lat);
    q=find(lat>=Lat(i)-dLat/2 & lat<=Lat(i)+dLat/2);
    for j=1:size(Lon,1);
        
        qq=find(lon(q)>=Lon(j,1) & lon(q)<=Lon(j,2));
        qqd=find(lon(q)>=Lon(j,1) & lon(q)<=Lon(j,2) & drogue(q)==1);
        qqu=find(lon(q)>=Lon(j,1) & lon(q)<=Lon(j,2) & drogue(q)==0);
        qqf=find(lon(q)>=Lon(j,1) & lon(q)<=Lon(j,2) & drogue(q)==-1);
        
        NumObs(i,j,1)=length(qqd)/4;    % drifter-days
        NumObs(i,j,2)=length(qqu)/4;
        NumObs(i,j,3)=length(qqf)/4;
        
        if (length(qq)>=12)
            
            ft=(t(q(qq))-7671)/365.24;  % day 7671: 1 Jan 2000
            fx=lon(q(qq))-mean(Lon(j,:));
            fy=lat(q(qq))-Lat(i);
            model=[1+0*ft,sin(ft*2*pi),cos(ft*2*pi) ...
                    sin(ft*4*pi) cos(ft*4*pi) fx fy]; 
                
            % regress zonal speeds in bin onto model coefficients    
            qu=pinv(model)*u(q(qq));
            U(i,j,1)=qu(1); 
            U(i,j,2)=std(u(q(qq)));
            U(i,j,3)=length(qq)/4;  % drifter-days 
             
            % save monthly values
            Umon(i,j,:)=model_mon*qu(1:5); 

            % regress meridional speeds in bin onto model coefficients    
            qu=pinv(model)*v(q(qq));
            V(i,j,1)=qu(1);
            V(i,j,2)=std(v(q(qq)));
            V(i,j,3)=length(qq)/4;
            
            % save monthly values
            Vmon(i,j,:)=model_mon*qu(1:5); 
            
            
            FracNoDrogue(i,j)=length(qqu)/length(qq);
            
            ft=t(q(qqd))/365.24;  
            fx=lon(q(qqd))-mean(Lon(j,:));
            fy=lat(q(qqd))-Lat(i);
            model=[1+0*ft,sin(ft*2*pi),cos(ft*2*pi) ...
                    sin(ft*4*pi) cos(ft*4*pi) fx fy];  
            qu=pinv(model)*u(q(qqd));
            Ud(i,j,1)=qu(1); 
            Ud(i,j,2)=std(u(q(qqd)));
            Ud(i,j,3)=length(qqd)/4;  % drifter-days  
            qu=pinv(model)*v(q(qqd));
            Vd(i,j,1)=qu(1);
            Vd(i,j,2)=std(v(q(qqd)));
            Vd(i,j,3)=length(qqd)/4;
            
            ft=t(q(qqu))/365.24;  
            fx=lon(q(qqu))-mean(Lon(j,:));
            fy=lat(q(qqu))-Lat(i);
            model=[1+0*ft,sin(ft*2*pi),cos(ft*2*pi) ...
                    sin(ft*4*pi) cos(ft*4*pi) fx fy];  
            qu=pinv(model)*u(q(qqu));
            Uu(i,j,1)=qu(1); 
            Uu(i,j,2)=std(u(q(qqu)));
            Uu(i,j,3)=length(qqu)/4;  % drifter-days  
            qu=pinv(model)*v(q(qqu));
            Vu(i,j,1)=qu(1);
            Vu(i,j,2)=std(v(q(qqu)));
            Vu(i,j,3)=length(qqu)/4;
        end;
        
        
    end;
end;

clear ans i j q qq

subplot(211);
plot([-23 0],[2 2],'k'); 
hold on;
plot([-23 0],[-2 -2],'k'); 
plot([-23 -23],[2 -2],'k');
plot([-10 -10],[2 -2],'k');
plot([0 0],[2 -2],'k');
quiver(mean(Lon'),Lat,squeeze(U(:,:,1)),squeeze(V(:,:,1)))
hold off;
set(gca,'dataaspectratio',[1 1 1]);

% calculate zonal divergence for the two boxes
%  (negative= convergence into box) using values
%  on box edges at 23/10/0W
dU=NaN*ones(2,12);  dV=dU;
for mon=1:12;
  % 23-10W, 2S-2N
    dU(1,mon)=-sum(Umon(:,1,mon).*mean(diff(Lat))*111e3*30) ...
        + sum(Umon(:,4,mon).*mean(diff(Lat))*111e3*30);
  % 10W-0, 2S-2N
    dU(2,mon)=-sum(Umon(:,4,mon).*mean(diff(Lat))*111e3*30) ...
        + sum(Umon(:,7,mon).*mean(diff(Lat))*111e3*30);
end;

% calculate meridional divergence for the two boxes
%  (negative= convergence into box) using values
%  averaged over west and east halves of boxes at 2N and 2S
for mon=1:12;
   % 23-10W, 2S-2N
    dV(1,mon)=(+Vmon(length(Lat),2,mon)+Vmon(length(Lat),3,mon) ...
        -Vmon(1,2,mon)-Vmon(1,3,mon)) ...
        .*diff(Lon(2,:)*111e3*30);
   % 10W-0, 2S-2N
    dV(2,mon)=(+Vmon(length(Lat),5,mon)+Vmon(length(Lat),6,mon) ...
        -Vmon(1,5,mon)-Vmon(1,6,mon)) ...
        .*diff(Lon(5,:)*111e3*30);   
end;

subplot(2,2,3); 
plot(tmon,dU(1,:)/1e6,'--b',tmon,dV(1,:)/1e6,':r',tmon,(dU(1,:)+dV(1,:))/1e6,'k');
title('zonal (b), meridional (r), tot (k), 23-10W')
set(gca,'xtick',tmon,'xticklabel',num2str( (1:12)'))

subplot(2,2,4); 
plot(tmon,dU(2,:)/1e6,'--b',tmon,dV(2,:)/1e6,':r',tmon,(dU(2,:)+dV(2,:))/1e6,'k');
title('zonal (b), meridional (r), tot (k), 10W-0')
set(gca,'xtick',tmon,'xticklabel',num2str( (1:12)'))


            



