function [LWhly, LWdly] = get_LW(lat,lon,elv,stdvc,fidvc,time_in,prcp,temp,rhum,isw,LWclr_i,LWall_i)
% function [LWhly, LWdly] = get_LW(lat,lon,elv,stdvc,fidvc,time_in,prcp,temp,rhum,isw,LWclr_i,LWall_i)
%
% estimating longwave from empirical approaches developed in ~/matlab/LongWaveRadiation
%
% LongWaveRadiation code. script name is "inc_LW_e_all.m".
%
%   [LWdwn,e_all,e_clr] = inc_LW_e_all(M_INPUTS, M_PARAMS, M_OPTIONS)
%
%   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
%   M_PARAMS = structure with the coefficients used for the method.  If you
%   	want to use the defaults, set M_PARAMS.P1 = []; or simply exclude the
%   	variable in the structure.
%   M_OPTIONS = structure with options for the LW method used
%	Method_vpr = 2 = calculate with Clausius-Clapeyron e_sat in mb (hPa) from Murray 1967
%       Method_clr = enter code for clear-sky longwave method, where:
%           10 = Prata (1996)
%           13 = Dilley and O'Brien (1998)
%       Method_all = enter code for all-sky longwave method, where:
%	    7  = Unsworth and Monteith (1975)

% gmauger Mar 2020

% ------------------------------------------------------------------------------------------------
% PREP INPUTS:
path(path,'/home/disk/margaret/mauger/matlab/LongWaveRadiation/lwscript/lwscript'),

M_PARAMS.STA_Elev = elv;

% time arrays:
%M_INPUTS.TIME = time_builder(stdvc(:,1),stdvc(:,2),stdvc(:,3),stdvc(:,4),stdvc(:,5),...
%			     fidvc(:,1),fidvc(:,2),fidvc(:,3),fidvc(:,4),fidvc(:,5),1);
%if rm_lp==1
  % REMOVE LEAP DAYS:
%  M_INPUTS.TIME(find(M_INPUTS.TIME(:,2)==2 & M_INPUTS.TIME(:,3)==29),:) = [];
%end

M_INPUTS.TIME = time_in;	% SWITCH TO USING TIME AS AN INPUT VAR

if size(M_INPUTS.TIME,1) ~= length(temp)
  whos, error('TIME and data time series not of same length'),
end
%M_INPUTS, whos temp prcp isw,

ind = TimeID(M_INPUTS.TIME,'DAY');                                   % Daily indices
ind(:,2) = [];

% Pressure. Assume constant, based on elevation. [kPa]
pressure=((101325*((1-(.0000225577*(3.281*M_PARAMS.STA_Elev)))^5.25588)) *  ...
	ones(size(M_INPUTS.TIME,1),1))./1000;

% convert prcp to cm/hr from mm/hr
prcp = prcp/10;

M_INPUTS.Ta = temp;
M_INPUTS.RH = rhum;
disp(['    ' datestr(now) ' -- done w MINPUTS (met)']),

% --------------------------------------------
% RADIATION:

% Average elevation angle [degrees above horizon]
%M_INPUTS,
El=AVG_EL(M_INPUTS.TIME,lat,lon,8,'BEG');
%whos El,
disp(['    ' datestr(now) ' -- done w AVG_EL']),

% generate a cloud fraction (0-1)
M_INPUTS.c = cloudfactor_Jessica(M_INPUTS.TIME,lat,isw,1);
disp(['    ' datestr(now) ' -- done w cloudfactor_Jessica']),

% Optical [depth?] air mass at 101.3 kPa
m = 35.*sind(El).*((1224.*(sind(El).*sind(El))+1).^-0.5);

% Empirical expression for the product of the first two transmission coeffs (Rayleigh and permanent gases)
%whos m pressure,
tau_r_tau_pg = 1.021-(0.084.*sqrt(m.*(0.00949.*pressure+0.051)));

% transmission coeffecient for absorption by water vapor
tau_w = 1-(0.077.*((prcp.*m).^0.3));

% transmission coeffecient for absorption by aerosols 
tau_a = 0.945.^m;


% ----------------------------------------------------------------------------------------
% calc LW:

% initialize LW arrays:
LWhly = NaN(size(M_INPUTS.TIME,1),length(LWclr_i),length(LWall_i));
LWdly  = NaN(length(ind),length(LWclr_i),length(LWall_i));

M_PARAMS.P1 = [];
M_OPTIONS.Method_vpr = 2;

for iclr = 1:length(LWclr_i)
for iall = 1:length(LWall_i)
  M_OPTIONS.Method_clr = LWclr_i(iclr);
  M_OPTIONS.Method_all = LWall_i(iall);

  [jnk,e_all,e_clr] = inc_LW_e_all(M_INPUTS, M_PARAMS, M_OPTIONS);
  LWhly(:,iclr,iall) = jnk;
  jnkdly	      = TimeAverage(M_INPUTS.TIME,LWhly,'day');
  LWdly(:,iclr,iall) = jnkdly(ind);
  clear jnk jnkdly e_all e_clr,

  disp(['    ' datestr(now) ' -- ' ...
	num2str([iclr length(LWclr_i) iall length(LWall_i)],'%d of %d and %d of %d')]),
end
end

disp(['    ' datestr(now) ' -- DONE.']),
