% General Settings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;                % Refresh command window
clear;              % Clear all variables
close all;          % Close all figures
format long g       % 14 decimal places
sec = 1; m = 60*sec; h = 60*m; d = 24*h; y = 365*d;      % time units (in seconds)
addpath(genpath('../../Functions')) % Load functions
physicalConstants

monsPlot  = 0;
printPlot = 0;

seasID = {'summer' 'winter'};


mainPath  = '';
inputDir  = './Input/';
inputFile = 'FLX_DE-Tha_FLUXNET2015_SUBSET_HH_1996-2014_1-3.csv';

disp('Read raw data..')
DATA = importdata([mainPath inputDir inputFile]);

disp('Process raw data..')
varIndex = [3 5 6 8 10 12 14 16 18 19 20 21 22 24 26 27 30 42 46 48 54];
varNames = DATA.textdata(varIndex);
varNum   = length(varNames);
Data     = DATA.data(:,varIndex);

startStr   = num2str(DATA.data(1,1));
firstStr   = num2str(DATA.data(1,2));
finalStr   = num2str(DATA.data(end,1));
dateFormIn = 'yyyymmddHHMM';

startNum   = datenum(startStr,dateFormIn);
firstNum   = datenum(firstStr,dateFormIn);
finalNum   = datenum(finalStr,dateFormIn);

startDate  = datestr(startNum);
firstDate  = datestr(firstNum);
finalDate  = datestr(finalNum);

dt         = round((firstNum-startNum)*d);  % time step in seconds
dtd        = dt/d;
nextNum    = finalNum + dtd;
dayNum     = nextNum-startNum;
ndt        = dayNum/dtd;
dtTotal    = size(Data,1);
timeVec    = linspace(startNum,finalNum,ndt)';

startRun    = datenum('01-Jan-2006'); 
finalRun    = datenum('31-Dec-2014 23:30:00'); 
startRunIdx = find(abs(timeVec - startRun)<1e-7);
finalRunIdx = find(abs(timeVec - finalRun)<1e-7);

rowTimes = datetime(timeVec(startRunIdx:finalRunIdx),'ConvertFrom','datenum');
Data     = DATA.data(startRunIdx:finalRunIdx,varIndex);
Data(Data==-9999) = NaN;
raw      = array2timetable(Data,'RowTimes',rowTimes,'VariableNames',varNames);



raw.T_SFC = (raw.LW_OUT./(stbo *emiss)).^0.25;
raw.P_F   = raw.P_F/dt;     % Precipitation   unit mm   --> kg/(m^2s) 
raw.PA_F  = raw.PA_F*1000;  % Air pressure    unit kPa  --> Pa
raw.TA_F  = raw.TA_F + T0;  % Air temperature unit degC --> K
raw.RH    = raw.RH/100;     % Relative humid  unit %    --> -
raw.Q_AIR = RH2q(raw.TA_F,raw.PA_F,raw.RH);
raw.VP    = raw.VPD_F./(1./raw.RH-1)*100;
[~,raw.E] = RH2q(raw.TA_F,raw.PA_F,raw.RH);
raw.Q_SAT = qsatNew(raw.TA_F,raw.PA_F);
[~,raw.E_SAT] = qsatNew(raw.TA_F,raw.PA_F);

raw.T_SFC = round(raw.T_SFC*100)/100;
raw.Q_AIR = round(raw.Q_AIR*100000)/100000;
raw.VP    = round(raw.VP*100)/100;
raw.E     = round(raw.E*100)/100;
raw.Q_SAT = round(raw.Q_SAT*100000)/100000;
raw.E_SAT = round(raw.E_SAT*100)/100;


detNaN = mean([raw.NETRAD raw.T_SFC raw.TA_F raw.Q_AIR raw.G_F_MDS raw.H_F_MDS raw.LE_F_MDS raw.SW_OUT raw.LW_OUT],2);
idxNaN = isnan(detNaN);

raw.NETRAD(idxNaN)   = NaN;
raw.SW_IN_F(idxNaN)  = NaN;
raw.SW_OUT(idxNaN)   = NaN;
raw.LW_OUT(idxNaN)   = NaN;
raw.LW_IN_F(idxNaN)  = NaN;
raw.T_SFC(idxNaN)    = NaN;
raw.TA_F(idxNaN)     = NaN;
raw.Q_AIR(idxNaN)    = NaN;
raw.G_F_MDS(idxNaN)  = NaN;
raw.H_F_MDS(idxNaN)  = NaN;
raw.LE_F_MDS(idxNaN) = NaN;

posIdx = raw.NETRAD>0;
negIdx = raw.NETRAD<=0;

raw.G = raw.G_F_MDS;
raw.G(posIdx) = 0.1*raw.NETRAD(posIdx);
raw.G(negIdx) = 0.3*raw.NETRAD(negIdx);

imbal = raw.NETRAD - raw.G - raw.LE_F_MDS - raw.H_F_MDS;

raw.LE = raw.LE_F_MDS;
raw.H  = raw.H_F_MDS;


raw.LE(posIdx) = raw.LE(posIdx) + imbal(posIdx);
raw.H(negIdx)  = raw.H(negIdx) + imbal(negIdx);

imbalTest = raw.NETRAD - raw.G - raw.LE - raw.H;


disp('Save raw data...')
save('./Output/thar.ObsNew.raw.mat','raw')

