%% Functions needed

    % TimeID.m
    % TimeAverage.m
    % time_builder.m
    % AVG_EL
    % inc_LW_e_all.m
    % inc_LW_e_clr.m


%% Inputs

time_hourly=time_builder(1980,12,31,2015,12,31,1);
Lat = 48.05700;            % decimal degrees
Lon = -122.03726;
%48.05700,-122.03726,103.3,71.6,48.05700,-122.04000,data_48.05700_-122.03730

Elev = 103.3;           % feet
SW_TR_derived = XX ; % [W m-2] derived SW data
Tadj = XX;           % temperature [deg C]
RH   = XX;           % decimal [0-1]
Ppt_hourly = XX;     % mm/hr

%% Code

ind = TimeID(time_hourly,'DAY');                                   % Daily indices
ind = ind(1:end,1);

El(:,jj,1)=AVG_EL(time_hourly,Lat(1,1),-Lon(1,1),8,'MID');        % Average elevation angle [degrees above horizon]

LWdown       = nan(length(time_hourly(:,7)),130);                  % Hourly LongWave Generation
LWday        = nan(length(Met.site_names(1).time(:,7)),130);       % Daily LW


        % SYNTAX
        %   [LWdwn,e_all,e_clr] = inc_LW_e_all(M_INPUTS, M_PARAMS, M_OPTIONS)
        %
        % INPUTS
        %   M_INPUTS = structure with the following variables (must have exact name)
        %       TIME = time matrix (per time_builder.m format), same number of rows as other inputs
        %       Ta = air temperature, K or C (will assume C if mean is less than 40 C)
        %       eo = vapor pressure (kPa) - optional if RH is provided and option Method_vpr>0
        %       RH = fractional relative humidity
        %       c = time series of cloud fraction (0-1 range) for Method_all from 1-8
        %       s = solar index (0-1 range) for Method_all from 9-10
        %       tau = transmissiivty (0-1 range) for method 16

        [c]=cloudfactor_Jessica(time_hourly(:,7),Lat,SW_TR_derived,1);                                % generate a cloud fraction (0-1)
        M_INPUTS.c=c;
        
        % Compute solar index
        pressure=((101325*((1-(.0000225577*(in_to_m(Elev.*12))))^5.25588))*ones(length(El),1))./1000; % assume a constant pressure, not measured at sites, [kPa]
        
        
        % As described in the Craford and Duchon 1999 method
        % s_clr=S_0 (sin(alpha)) T_r T_pg T_w T_a
        % S_0   is the solar constant taken as 1360 W m-2
        % alpha is the angle of the sun above the horizon
        % T_r   is the transmission coeffeicnet for rayleigh scattering
        % T_pg  is the transmission coefffeicent for water vapor
        % T_a   is the transmission coeffeicnet for aerosol scattering
        
        m            = 35.*sind(El).*((1224.*(sind(El).*sind(El))+1).^-0.5);         % Pptical air mass at 101.3 kPa
        tau_r_tau_pg = 1.021-(0.084.*sqrt(m.*(0.00949.*pressure+0.051)));            % Empirical expression for the product of the first two transmission coeffecients (Rayleigh and permanent gases)
        w            = (Ppt_hourly.*3600);                                      % Convert to cm/hr from kg/m2 s
        tau_w        = 1-(0.077.*((w.*m).^0.3));                                     % transmission coeffecient for absorption by water vapor 
        tau_a        = 0.945.^m;                                                     % transmission coeffecient for absorption by aerosols 
        s_clr        = 1360.*sind(El).*tau_r_tau_pg.*tau_w.*tau_a;                   % as defined in Appendix A of Flechinger
        
        count=0;                    % Keeps track of the longwave combination 
        for clrMethod=1:13
            for allMethod=1:10
                
                count=count+1;
                M_INPUTS.TIME        = time_hourly;
                M_INPUTS.Ta          = Tadj;
                M_INPUTS.RH          = RH./100;
                M_INPUTS.s           = SW_TR_derived./s_clr;
                M_INPUTS.tau         = SW_TR_derived./1360;     % Atmospheric transmissivity: a common index of cloud-cover taken from Sicart et al. 2006
                
                M_PARAMS.STA_Elev    = in_to_m(Elev.*12);
                M_PARAMS.P1          = [];
                
                M_OPTIONS.Method_vpr = 2;
                M_OPTIONS.Method_clr = clrMethod;
                M_OPTIONS.Method_all = allMethod;
                
                [LWdwn,e_all,e_clr]  = inc_LW_e_all(M_INPUTS, M_PARAMS, M_OPTIONS);  % Derive LW
                
                LWdown(:,count)= LWdwn;                                        % Hourly LW
                LWdayTmp       = TimeAverage(time_hourly,LWdwn,'day');                    
                LWday(:,count) = LWdayTmp(ind);                                % Mean daily LW
                
                
            end
        end
