%% Functions needed

    % TimeID.m
    % TimeAverage.m
    % time_builder.m
    % AVG_EL
    % inc_LW_e_all.m
    % inc_LW_e_clr.m
    % in_to_m
    % greg2julian
    % isarray



%% Inputs

%select WRF grid cell (0 to 54)
gn = 54
%Enter begin and end time of PNNL data
time_hourly=time_builder(1981,01,01,0,0,2015,12,31,23,0,1);

%load data for all WRF grid cells
cd 'D:\UW_PhD_files\LongWaveRadiation\lwscript\input'
P = load('Pmat.mat'); %precipitation [m/hr]
RH = load('RHmat.mat'); %relative humidity [%]
SW = load('SWmat.mat'); %short wave radiation [W/m2]
T = load('Tmat.mat'); %temperature [deg C]
LW = load('LWmat.mat'); %long wave radiation from WRF [W/m2]
LD = readtable('Sauk_P_mappingfile_i.csv', 'HeaderLines',1); 


%for each pnnl grid cell (right now manually iterate by changing gn and
%rerunninng script)

%parameters
field = ['key' num2str(gn)]
Lat = LD(gn+1,4).Var4 ;           % [decimal degrees] 
Lon = LD(gn+1,5).Var5 ;           % [decimal degrees] 
Elev = LD(gn+1,6).Var6 *3.281;           % [feet]
%variables
SW_TR_derived = SW.(field).' ; % [W m-2] derived SW data %short wave radiation ->from pnnl wrf
Tadj = T.(field).';           % temperature [deg C] ->from pnnl wrf
RH   = RH.(field).'/100;           % [decimal] [0-1] ->estimated from pnnl wrf using Feld et al., 2013 
Ppt_hourly = P.(field).'*1000;     % [mm/hr] from pnnl wrf
LW_WRF = LW.(field).'; %[W m-2]


%% 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]
El=AVG_EL(time_hourly,Lat(1,1),-Lon(1,1),8,'MID');        %j Average elevation angle [degrees above horizon]
nh = length(time_hourly(:,7)); %number of hours in time series
LWdown       = nan(nh,130);                  % Hourly LongWave Generation
LWday        = nan(nh/24,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]
        pressure=((101325*((1-(.0000225577*(3.281*Elev)))^5.25588))*ones(length(El),1))./1000; %j 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.STA_Elev    = 3.281*(Elev); %j
                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
                
                disp('clrMethod')
                disp('allMethod')
            end
        end
        



%% Check mean longwave radiation of all model runs:
mnLWdown = nanmean(LWdown,1)
%figure size
x0=10;
y0=10;
width=1200;
height=400;

cd D:\UW_PhD_files\LongWaveRadiation\lwscript\output

figure
set(gcf,'position',[x0,y0,width,height])
plot(linspace(1,130,130),mnLWdown,'r')
hold on
plot([0,140],[mean(LW_WRF),mean(LW_WRF)],'k--')
ylabel('mean LW [W/m2]')
xlabel('run #')
legend('model run','PNNL WRF')
grid on
grid minor
saveas(gcf,['MeanLW_g' num2str(gn) '.jpg'])

%% View results of specific model runs:

%model run #127
%Method_clr
%13 = Dilley and O'Brien (1998)
%Method_all
%7 = Unsworth and Monteith (1975)

%model run #67
%Method_clr
%10 = Prata (1996)
%Method_all
%7 = Unsworth and Monteith (1975)

%model run #82
%Method_clr
%Niemela et al. (2001)
%Method_all
%Iziomon et al. (2003)
%compute mean of each run



%figure size
x0=10;
y0=10;
width=800;
height=400;

figure
set(gcf,'position',[x0,y0,width,height])
plot(time_hourly(:,7),LWdown(:,127),'r')
hold on
plot(time_hourly(:,7),LWdown(:,67),'b')
plot(time_hourly(:,7),LWdown(:,82),'g')
plot(time_hourly(:,7),LW_WRF,'k')
ylabel('LW (W/m2)')
legend('D&O1998 clr, U&M1975 allLW','Prata1996 clr, U&M1975 allLW' ,'Niemela2001 clr, Iziomon2003 allLW','PNNL WRF')
axis([datenum(2008,10,1) datenum(2009,3,30) -inf inf])
datetick('x','yy/mm/dd','keeplimits')
saveas(gcf,['LWCompare_g' num2str(gn) '.jpg'])