function y = scs (timerain, CNII, evaprate, data4)
%
% Accepts a 2-column array, timerain, with columns for time and ppt
% Sd, the storage depth, and evaprate, the rate of evaporation.
% It goes through the precipitation table and takes off the storage
% depth for the entire array.
% The evaporation rate must be in compatible units with Sd and timerain
% (e.g., if Sd is in inches and the time step in the timerain table is 
% hours, then evaprate must be in inches per hour).
%
Pi=0; Wi=0; Wo=0;
antRain = 0;
datapts = size(timerain,1); % Determines the number of data points
ARCIII = [0 17 26 33 39 45 50 55 60 65 70 75 79 83 87 91 94 97 98 99 100];
ARCII = [0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100];
ARCI = [0 2 4 7 9 12 15 19 23 27 31 35 40 45 51 57 63 70 78 87 100];
CNIII = interp1(ARCII,ARCIII,CNII);
CNI=interp1(ARCII,ARCI,CNII);
for i = 1:datapts
    newtimerain(i,1) = timerain(i,1); % time
    if i < 481
        antRain = antRain + timerain(i,1);
    else
        antRain = antRain + timerain(i,1) - timerain(i-480,1);
    end
    if data4(i,3)==0
        if antRain > 5.33
            CN = CNIII;
        elseif antRain < 3.55
            CN = CNI;
        else
            CN = ((antRain-3.55)/(1.78)*(CNIII-CNII))+CNI;
        end
    end
    Smax = 2.54*(1000/CN-10);
    if data4(i,3)==1 
        Wo=Wi;
        Pi=Pi+timerain(i,2);  %Cumm. precip counter
        Wi=((Pi-0.2*Smax)^2)/(Pi+0.8*Smax);  %SCS Cumm. Runoff Formula
        if Pi<0.8*Smax
            Wi=0;
        end
        newtimerain(i,2)=Wi-Wo;  %incr. Runoff
    end
    if data4(i,3)==0 
        Pi=0; Wi=0; Wo=0;
        newtimerain(i,2) = timerain(i,2);% reset to precip counters to zero
    end  
end
y = newtimerain;