function x=adjout(y,thr,tflag)

%    -- Adjust for outliers using fraction of IQR
%
%       y = Data series
%       thr = threshold in multiples of IQR
%       tflag = 0  == replace with missing value
%               1  == replace with maximum value
%               2  == replace with median value
%               3  == replace with local median (obs + or - 3 on each side)
%               4  == replace with one-sided median (5 preceding obs)

small = 1.0e-06;

% -- Compute IQR -- %
z=y(~isnan(y));
z=sort(z,1);
nzrow=size(z,1);
zm=z(ceil(0.5*nzrow));
iqr=z(ceil(.75*nzrow))-z(ceil(.25*nzrow));

if iqr < small
    x=nan(size(y));
    return
end

ya=abs(y-zm);

x=y;
iya = ya > (thr*iqr);
iyb = ya <= (thr*iqr);

switch tflag
    case 0
        x=(iyb .* y) + (iya .* nan);
    case 1
        isign = y > 0;
        jsign = -(y < 0);
        isign=isign+jsign;
        yt=(zm.*ones(length(y),1)) + isign .* (thr .* ones(length(y),1));
        x=(iyb .* y) + (iya .* yt);
    case 2
        x=(iyb .* y) + (iya .* zm);
    case 3
        % Compute rolling median %
        iwin=3;  % Window on either side %
        ymvec=nan(length(y),1);
        for i=1:length(y)
            j1=max(1,(i-iwin));
            j2=min(length(y),(i+iwin));
            ymvec(i)=nanmedian(y(j1:j2));
        end
        x=(iyb .* y) + (iya .* ymvec);
    case 4
        % Compute rolling median %
        iwin=5;  % Window on ones side %
        ymvec=miss(zeros(rows(y),1),0);
        for i=1:size(y,1)
            j1=maxc(1,(i-iwin));
            j2=i;
            ymvec(i)=nanmedian(y(j1:j2));
        end
        x=(iyb .* y) + (iya .* ymvec);
    otherwise
        error('choose tflag in {0,1,2,3,4}!')
end

end
