function qc_sst(filename);

% function qc_sst(filename);
% experimental code to quality-control
% the sst observations made by the drifters.
%   IDs are stored in C5 (category 5), C3 and C1.
%   1=odd-looking, 3=suspect, 5=major problems.

% calcreynoldssst must be run first, to interpolate
%  weekly Reynolds SST to the drifter positions
%  (output is "sst_rey")

% load data
eval(['load ',filename,' lon lat t id sst sst_rey']);

del_thresh=.75;  % threshold for "bad" sst difference
del_t=median(diff(t(1:100)));

% calculate D, the difference between 
%  drifter SST and Reynolds SST
D=sst-sst_rey;
% calculate dtD, the rate of change of D (used
% to identify spikes; units are delta-degC/s.
q=find(diff(id));
q=[0;q;length(id)];
dtD=NaN*ones(size(D));
for i=2:length(q);
    a=q(i-1)+1:q(i);
    q2=find(isfinite(sst(a)-sst_rey(a)));
    if (length(q2)>4);
        sstlp=NaN*ones(size(a));
        sstlp(q2)=sst(a(q2))-sst_rey(a(q2));
        dtD(min(a)+2:max(a)-2)= ...
            (sstlp(5:length(a))-sstlp(1:length(a)-4));
        dtD(min(a)+2:max(a)-2)= ...
            dtD(min(a)+2:max(a)-2) ./ ...
            (t(a(5:length(a)))- ...
                t(a(1:length(a)-4))) /3600/24;
    end;
end;

% create vectors corresponding to the list of unique
% IDs to store qc metrics:
%   md: the fraction of lifetime for which 
%       abs(D)>del_thresh.
%   sd: the lifetime-mean standard deviation of D
%   T:  the lifetime (days)
%   TD: lifetime multiplied by lifetime-mean abs(D)
ID=unique(id);
md=NaN*ones(size(ID));
sd=md;T=md;TD=md;

for i=1:length(ID);
    fprintf('%d of %d\n',i,length(ID));
    q=find(id==ID(i));
    q2=find(isfinite(D(q)));
    if (length(q2)>2)    
        T(i)=length(q2)*del_t;  % number of days of sst obs
        md(i)=length(find(abs(D(q))>del_thresh))/length(q2);
        sd(i)=std(D(q(q2)));
        TD(i)=T(i)*abs(mean(D(q(q2))));
    end;
end;

% Now identify the problem drifters
% CAT 5: SST (drifter) and SST (Reynolds) differ
% by more than 10 degrees at any time,
%    OR  TD>200
%    OR  d_sst>.5e-4 dgC/s
%    OR  std.dev>10 dgC
%    OR  D always exceeds threshold and T>30d
C5=find(abs(D)>10 | abs(dtD)>.5e-4);
C5=unique(id(C5));
C5=unique([C5;ID(find(TD>200)); ...
    ID(find(md==1 & T>30));ID(find(sd>10))]);

% CAT 3: mean exceeds threshold >99% of the time,
%  OR std.dev. is greater than 3degC,
%  OR TD>150,
%     AND (in all cases) T>10 days.
C3=find( T>10 & ...
    (md>.99  | sd>3 | TD>150) );
C3=ID(C3);
% remove redundancies (IDs in C5)
for i=1:length(C5);
    q=find(C3==C5(i));
    C3(q)=NaN;
end;
C3=C3(find(isfinite(C3)));

% CAT 1: mean exceeds threshold >50% of the time,
%  AND std.dev. is greater than 0.5 degC.
C1=find(md>.5 & sd>0.5 );
C1=ID(C1);
% remove redundancies (IDs in C5 and C3)
for i=1:length(C5);
    q=find(C1==C5(i));
    C1(q)=NaN;
end;
for i=1:length(C3);
    q=find(C1==C3(i));
    C1(q)=NaN;
end;
C1=C1(find(isfinite(C1)));

fprintf('CAT 5 (probably bad SST): [');
for i=1:length(C5);
    fprintf(' %d',C5(i));
    if (i/5==round(i/5)) fprintf('\n     '); end;
end;
fprintf(']\n');
fprintf('CAT 3 (questionable SST): [');
for i=1:length(C3);
    fprintf(' %d',C3(i));
    if (i/5==round(i/5)) fprintf('\n     '); end;
end;
fprintf(']\n');
fprintf('CAT 1 (odd/unusual SST): [');
for i=1:length(C1);
    fprintf(' %d',C1(i));
    if (i/5==round(i/5)) fprintf('\n     '); end;
end;
fprintf(']\n');
