% Julia Casagrande
% BEE 1510: Intoduction to Computer Programming
% Laboratory Exercise #6
% Title: Rainfall Excess and Surface Runoff: Files and Arrays
% Purpose: determine snow depth and snow melt
%          for each day of the year, based on the weather data
% Laboratory TA: Maureen Murage
%                Graham Nekut
% Last Modified: 10/13/10

% Variable Dictionary
    % precip_data = precipitation data file
    % precip_read = function to read precipitation data
    % A = array of precipitation data
    % temp_data = temperature data file
    % tp = function to read temperature
    % B = array of temperature data
    % C = matrix of temperature data when precipitation is rain
    % D = matrix of temperature data when precipitation is rain
    % isFrozen = days when ground is frozen
    % isDormant = days when ground is dormant
    % day = day of the year
    % K = snow melta parameter; .45
    % H2O_shed_data = water shed data file
    % H2Oread = function to read water shed data
    % E = array of water shed data
    % Q_ti = depth of water runoff
    % S_ti = potential retention of precipitation (max soil water
    %        retention parameter)
days = 1:365;

precip_data = '151PRECP.DAT';
precip_read = readPrecipitation(precip_data);

temp_data = '151TEMP.DAT';
tp = readTempFile(temp_data);

rain = find(tp >= 0);
snow = find(tp < 0);


rainfall = zeros(365,1);
rainfall(rain) = precip_read(rain);

snowfall = zeros(365,1);
snowfall(snow) = precip_read(snow);

isFrozen = false(365,1);
isDormant = false(365,1);

days_frozen = 120:305;
days_dormant = 121:304;

isFrozen(days_frozen) = true;
isDormant(days_dormant) = true;


isFrozen = days <= 120 | days >= 305;
isDormant = days >= 121 & days <= 304;

K = 0.45;
snow = 0;
snow_base = zeros(365,1);
snowmelt = zeros(365,1);
snowmelt_actual = zeros(365,1);

for days=[1:365]
    snowmelt(days) = max((K*tp(days)),0);
    snow = snow + snowfall(days) - snowmelt(days);
    snow_base(days) = snow;
    snowmelt_actual(days) = min(snowmelt(days), snow_base(days));
end

H2O_shed_data = 'H2OSHED.DAT';
[land_use_area, reference_curve_num] = readH2O(H2O_shed_data);
land_use_area_m = land_use_area.*((1000)^2);

am = antecedentMoisture(rainfall, snowmelt_actual);

% Calculate antecedent moisture
for days = 1:365
    runningsum = 0;
    for J = (days-1):(days-5)
       am = rainfall(J) + snowmelt_actual(J);
       runningsum = runningsum + am;
    end
end




for land_type = 1:length(reference_curve_num)

    cn = curveNumber(reference_curve_num(land_type), am, snowmelt_actual, isDormant, isFrozen);

    S(land_type, days) = (2540/cn(days))-25.4;

    for days = 1:365
        % Calculate daily curve numbers




        % Calculate depth of water runoff
        if precip(days)> 0.2*S(days)
            Q(days, land_type) = ((precip(days)- 0.2*S(days, land_type))^2)/(precip(days)- 0.8*S(days))*(land_use_area_m);
        else
            Q(days, land_type) = 0;
        end
    end
end

plot(Q)