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

clc;                % Refresh command window
%clear;              % Clear all variables
clearvars -except DATA
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






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

DATA = importdata([mainPath inputDir inputFile]);

varIndex = [3 5 6 8 10 12 14 16 18 19 20 21 22 24 26 27 30 42 46 48 54];
varID    = DATA.textdata(varIndex);
varNum   = length(varID);
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);

D.time = timeVec(startRunIdx:finalRunIdx);
for i = 1:varNum
    currVar = DATA.data(startRunIdx:finalRunIdx,varIndex(i));
    currVar(currVar==-9999) = NaN;
    D.(varID{i}) = currVar;
end

varNames  = fieldnames(D);
nVar      = length(varNames);

% for iVar = 1:nVar
%     figure
%     plot(D.(varNames{iVar}))
%     title(varNames{iVar},'interpreter','none')
% end

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

thar = D;
save('tharFluxnet.mat','thar')






