% bias-correct WRF forcing files for stilly grid cells 
%
% 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
%

% gmauger May 2020
% gmauger Jun 2020 -- added elevation correction for WRF-DEM elev difference

disp([datestr(now) ' -- Starting...']),
clear,
ist = 1;
ifi = 9999;

% ------------------------------------------------------------------------------------------------
lapse_rate = 0.0/1000;         % K/m

bc_T = 0.2;	% additive
bc_P = 1.15;    % multiplicative
bc_W = 0.6;	% multiplicative
bc_SW = 0.9;	% multiplicative

bc_T = 0;		% additive
bc_P = 1;    % multiplicative
bc_W = 1;		% multiplicative
bc_SW = 1;	% multiplicative



scrdir  = '/home/disk/margaret/mauger/2020_12_SnohoCounty_Flooding/scripts';
outbase = '/home/disk/rocinante/DATA/temp/WRF/scripts/regrid_c2c/longwave/output/';
rawdir  = '/home/disk/rocinante/DATA/temp/WRF/scripts/regrid_c2c/longwave/data/';
%rawdir  = '/home/disk/rocinante/DATA/temp/kcp3/data/prismWRFhr/';

%outdir_uni = [outbase '/bc_forcs/lapse_rate_' num2str(abs(lapse_rate)*1000,'%.1f') '/pnnl_wrf/unifm_bc'];
outdir_uni = [outbase '/bc_forcs/lapse_rate_' num2str(abs(lapse_rate)*1000,'%.1f') '/pnnl_wrf'];
%outdir_uni = [outbase '/prismWRF/'];

% ------------------------------------------------------------------------------------------------
% get lat/lon list:
%fid = fopen('/home/disk/rocinante/DATA/temp/WRF/scripts/regrid_c2c/longwave/data/filter_elev.csv', 'r');
fid = fopen('/home/disk/rocinante/DATA/temp/WRF/scripts/regrid_c2c/longwave/data/elev.csv', 'r');


jnk = textscan(fid,'%f %f %f %f %f %f %s','Delimiter',',','HeaderLines',1);
frclat = jnk{1};
frclon = jnk{2};
pnlelv = jnk{3};
demelv = jnk{4};
jsnfnm = jnk{7};
fclose(fid);
clear jnk fid,

% ------------------------------------------------------------------------------------------------
if ~exist(outdir_uni,'dir')
  mkdir(outdir_uni),
end

% ------------------------------------------------------------------------------------------------
for fi=ist:min(ifi,length(frclat))
  % ------------------------------------------------
  frc_fnm = ['data_' num2str(frclat(fi),'%.5f') '_' num2str(frclon(fi),'%.5f')];

  % read in raw forcings:
  frcfl = [rawdir '/' jsnfnm{fi}];
  if ~exist(frcfl,'file')
    frcfl, error('frcfl does not seem to exist. naming problem?'),
  end

%  if ~exist([outdir_uni '/' frc_fnm],'file')

    % ------------------------------------------------
    % get raw forcing data:
    fid = fopen(frcfl,'r');
    indat = textscan(fid,'%s %f %f %f %f %f %f','Delimiter',' ');
    dnm = datenum(indat{1}); %-8/24;       % FIX TIME ZONE
    dvc = datevec(dnm);

    % check for leap days:
    jnk = datevec(dnm+8/24);                    % start with GMT dates, since those follow the model cal
    if sum(jnk(:,2)==2 & jnk(:,3)==29) == 0     % remove 2/29 in time-adjusted dvc, leave alone if isleap
      disp('NO LEAP DAYS IN FORCINGS???'),
      ilp=find(dvc(:,2)==2 & dvc(:,3)==29);
      dvc(ilp,3)=28;
      clear ilp,
    else
      disp('CONFIRMING: PNNL forcings have leap days'),
    end
    clear jnk,

    dtstr = cellstr(datestr(dnm,'mm/dd/yyyy-HH:MM:SS'));
    temp  = indat{2};
    wind  = indat{3};
    rhum  = indat{4};
    isw   = indat{5};
    LWwrf = indat{6};
    prcp  = indat{7};
    fclose(fid);
    clear indat fid,

    % round up:
    i60=find(dvc(:,5)>55);
    dvc(i60,4) = dvc(i60,4) + 1;
    dvc(i60,5:6) = 0;
    clear i60,

    % round down:
    dvc(find(dvc(:,5)>0 & dvc(:,5)<5),5) = 0;
    dvc(find(dvc(:,6)>0),6) = 0;

    stdvc = datevec(dtstr{1});
    fidvc = datevec(dtstr{end});
    fidvc(:,5) = 1;     % to ensure final hour is counted. If zero it's sometimes dropped by time_builder

    % ------------------------------------------------
    % Elevation correction:
	elvcorr = (demelv(fi)-pnlelv(fi))*lapse_rate;	% adjusting temp to match demelv(fi)
    disp(num2str([demelv(fi) pnlelv(fi) elvcorr],'dem=%.0f, wrf=%.0f (m), elvcorr=%+.2f (C)')),

    % ------------------------------------------------
    % get LW (Dilley and O'Brien + Unsworth and Monteith (1975)):
    %   (needs to be done with the temp correction applied)
    [LWdob, jnk] = get_LW(frclat(fi),frclon(fi),demelv(fi),stdvc,fidvc,[dvc dnm],...
                                prcp,temp+bc_T+elvcorr,rhum,isw,13,7);
    clear jnk,
	% ------------------------------------------------
    % write to files:

    % Uniform deltas:

    fid = fopen([outdir_uni '/' frc_fnm],'w');
    for ti=1:length(dtstr)
      fprintf(fid,[dtstr{ti} ' ' num2str(temp(ti)+bc_T+elvcorr,'%.8f') ' ' ...
	num2str(wind(ti)*bc_W,'%.8f') ' ' num2str(rhum(ti),'%.8f') ' ' ...
	num2str(isw(ti)*bc_SW,'%.8f') ' ' num2str(LWdob(ti),'%.8f') ' ' ...
	num2str(prcp(ti)*bc_P,'%.8f') '\n']);
    end
	fclose(fid);

%  end

  clear fid dtstr temp wind rhum isw LWwrf LWprt LWdob prcp ti elvcorr,

  disp([datestr(now) ' -- done w: ' frc_fnm]),
end
clear fi,

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

