% test approach to estimating longwave from empirical approaches developed in ~/matlab/LongWaveRadiation
%
% forcing file format:
% Date Time, Air Temp (C), Wind (m/s), Rel. Hum. (%), Shortwave (W/m^2), Longwave (W/m^2), Precip (mm)
% 
% e.g.:
% 01/01/1981-16:00:00 9.49566650 1.78458345 98.73087373 0.00000000 351.39556885 0.00001283
% 01/01/1981-17:00:00 9.74948120 1.67217803 98.66921515 55.56172562 347.69848633 0.00000828
% 01/01/1981-18:00:00 9.93814087 2.08545280 98.62618649 150.49931335 340.99572754 0.00000625
% 01/01/1981-19:00:00 9.91506958 2.12686610 98.62954544 279.03155518 311.97860718 0.00000321
%
% 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:
%	    ????****????

% gmauger Mar 2020

disp([datestr(now) ' -- Starting...']),
clear,

% ------------------------------------------------------------------------------------------------
% *** TEMPORARY: LATER, THIS WILL BE A FUNCTION AND THE VARS WILL BE PROVIDED AS INPUTS ***

frcfl = '/home/disk/rocinante/DATA/temp/pnnl/westWAOR/4_output/forcs_dhsvm/data_47.80490_-121.75440';
%frcfl = '/home/disk/rocinante/DATA/temp/TNC_stormwater/bias_correction/scripts/bias_correction_wrf/lw/testdata';
iu=findstr(frcfl,'_');
lat = str2num(frcfl(iu(end-1)+1:iu(end)-1));
lon = str2num(frcfl(iu(end)+1:end));

% read in raw forcing file:
fid = fopen(frcfl,'r');
indat = textscan(fid,'%s %f %f %f %f %f %f','Delimiter',' ');
dtstr = indat{1};
temp  = indat{2};
wind  = indat{3};
rhum  = indat{4};
isw   = indat{5};
olw   = indat{6};
prcp  = indat{7};
fclose(fid);
clear indat fid,

stdvc = datevec(dtstr{1});
fidvc = datevec(dtstr{end});
disp([datestr(now) ' -- done reading in forcing file']),

% ------------------------------------------------------------------------------------------------
% PREP INPUTS:

M_PARAMS.STA_Elev = 1148.3;

% 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);

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]
El=AVG_EL(M_INPUTS.TIME,lat,lon,8,'MID');

% generate a cloud fraction (0-1)
M_INPUTS.c = cloudfactor_Jessica(M_INPUTS.TIME(:,7),lat,isw,1);

% 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)
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;

% as defined in Appendix A of Flechinger
s_clr = 1360.*sind(El).*tau_r_tau_pg.*tau_w.*tau_a;

M_INPUTS.s = isw./s_clr;
M_INPUTS.tau = isw/1360; % Atmos. Transmissivity: a common index of cloud-cover from Sicart et al. 2006
disp([datestr(now) ' -- done w MINPUTS (radiation)']),

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

% initialize LW arrays:
%LWclr_i = [10 13];
%LWall_i = 1:20;
LWclr_i = [13];
LWall_i = [7];
LWdown = NaN(size(M_INPUTS.TIME,1),length(LWclr_i),length(LWall_i));
LWday  = 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);
  LWdown(:,iclr,iall) = jnk;
  jnkdly	      = TimeAverage(M_INPUTS.TIME,LWdown,'day');
  LWday(:,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

% check corr/diffs
dhr = NaN(length(LWclr_i),length(LWall_i),3);
ddy = NaN(length(LWclr_i),length(LWall_i),3);
rhr = NaN(length(LWclr_i),length(LWall_i),3);
rdy = NaN(length(LWclr_i),length(LWall_i),3);
for iclr = 1:length(LWclr_i)
for iall = 2:length(LWall_i)
  dhr(iclr,iall,:) = prctile(LWdown(:,iclr,1)./LWdown(:,iclr,iall),[2.5 50 97.5]);
  ddy(iclr,iall,:) = prctile(LWday(:,iclr,1) ./LWday(:,iclr,iall), [2.5 50 97.5]);

  rhr(iclr,iall,:) = corr_err_fisher(LWdown(:,iclr,1),LWdown(:,iclr,iall),0.2/24);
  rdy(iclr,iall,:) = corr_err_fisher(LWday(:,iclr,1), LWday(:,iclr,iall), 0.2);
end
end
% all all-sky methods are identical on a daily time step, but can be very different at hourly time steps
% 	some are even negatively correlated with each other

rclrhr = NaN(length(LWall_i),3);
rclrdy = NaN(length(LWall_i),3);
for iall = 1:length(LWall_i)
  rclrhr(iall,:) = corr_err_fisher(LWdown(:,1,iall),LWdown(:,2,iall),0.2/24);
  rclrdy(iall,:) = corr_err_fisher(LWday(:,1,iall), LWday(:,2,iall), 0.2);
end
clear iall,
% same for the clear-sky methods -- differences are only at hourly time-step

rall1hr = NaN(length(LWall_i),length(LWall_i),3);
rall2hr = NaN(length(LWall_i),length(LWall_i),3);
for iall1 = 1:length(LWall_i)
for iall2 = 1:length(LWall_i)
  if iall1 ~= iall2
    rall1hr(iall1,iall2,:) = corr_err_fisher(LWdown(:,1,iall1),LWdown(:,1,iall2),0.2/24);
    rall2hr(iall1,iall2,:) = corr_err_fisher(LWdown(:,2,iall1),LWdown(:,2,iall2),0.2/24);
  end
end
end

squeeze(rall1hr(:,:,1)),

figure(1), clf,
hp=plot(squeeze(LWdown(1:24,1,:))); hl=legend(hp,cellstr(num2str([1:20]')));
set(hp(9),'linewidth',4),
set(hp(10),'linewidth',2),
set(hp(18),'linewidth',2),
set(hp(18),'linewidth',4),
set(hp(19),'linewidth',2),

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