%-----------------------%
%   Read JSBACH Ouput   %
%-----------------------%

clc;                % Refresh Command Window
clear;              % clear Workspace
close all;          % close all Figures
format long g       % 14 decimal places
addpath(genpath('../../Functions')) % Load functions
physicalConstants


siteName    = 'thar';             % choose between trop thar and dice
selSeas     = 'summer';

switch siteName
    case 'trop'
        siteFolder  = 'Tropics/'; 
    case 'thar'
        siteFolder  = 'Tharandt/';
    case 'dice'
        siteFolder  = 'DICE/';
    otherwise
        error('Site not avaiblable!')
end


obs = load('./Output/thar.ObsNew.proc.mat');

Tsobs = obs.proc.Cyc.(selSeas).T_SFC;
Hobs = obs.proc.Cyc.(selSeas).H;
LEobs = obs.proc.Cyc.(selSeas).LE;
Robs = obs.proc.Cyc.(selSeas).NETRAD;
Gobs = obs.proc.Cyc.(selSeas).G;
Nobs = Robs - Hobs - LEobs - Gobs;

Luobs = obs.proc.Cyc.(selSeas).LW_OUT;
Ldobs = obs.proc.Cyc.(selSeas).LW_IN_F;
Suobs = obs.proc.Cyc.(selSeas).SW_OUT;
Sdobs = obs.proc.Cyc.(selSeas).SW_IN_F;
Snobs = Sdobs - Suobs + Ldobs - Luobs;


inputFolder = './Output/';

surfScheme  = 'CEBa';
inputFile   = [siteName '.' surfScheme '.Ref.proc.mat'];
ceba(2) = load([inputFolder inputFile]);

surfScheme  = 'Clas';
inputFile   = [siteName '.' surfScheme '.Ref.proc.mat'];
ceba(1) = load([inputFolder inputFile]);


expNum = 2;
H = ceba(expNum).proc.Cyc.(selSeas).Hflx;
LE = ceba(expNum).proc.Cyc.(selSeas).Lflx;
R = ceba(expNum).proc.Cyc.(selSeas).Rnet;
G = ceba(expNum).proc.Cyc.(selSeas).Gflx;
S = ceba(expNum).proc.Cyc.(selSeas).Stot;
N = R + H + LE + G - S;

Lu = ceba(expNum).proc.Cyc.(selSeas).LWup;
Ld = ceba(expNum).proc.Cyc.(selSeas).LWdo;
Sd = ceba(expNum).proc.Cyc.(selSeas).SWdo;
Su = Sd+Ld-Lu-R;



dt    = 7.5*m;
dtd   = d/dt;
% Create time x-axis for all plots        
xStart = '01-Jan-2015';           % arbitary startDate
xFinal = '02-Jan-2015';           % arbitary finalDate
xS = datenum(xStart);
xF = datenum(xFinal)-1/dtd;
nDay = round(xF-xS);
x    = linspace(xS,xF,nDay*dtd);
xdt     = 4;                        % xtick step in Hours
xInt    = xdt*hr/d;
xTicks  = xS:xInt:xF;
xTicks  = [xTicks xTicks(1)+nDay];

xLabels = cellstr(datestr(xTicks,'HH:MM'));
xLabels = [xLabels;xLabels{1}];


dt    = 30*m;
dtd   = d/dt;
% Create time x-axis for all plots        
xStart = '01-Jan-2015';           % arbitary startDate
xFinal = '02-Jan-2015';           % arbitary finalDate
xS = datenum(xStart);
xF = datenum(xFinal)-1/dtd;
nDay = round(xF-xS);
xobs       = linspace(xS,xF,nDay*dtd);


colNum = [2 1];



figure
plot(xobs,Tsobs-T0,'--','Color','k','LineWidth',1.5)
hold on
for i = 1:length(ceba)
    plot(x,(ceba(i).proc.Cyc.(selSeas).LWup/emiss/stbo).^0.25-T0,'Color',newColor(colNum(i)),'LineWidth',2.5)
end

set(gca,'FontSize',22,'FontName','Verdana','xtick',xTicks,'xlim',[xS xTicks(1)+nDay],'xticklabel',xLabels)
ylabel('T_{rad} [°C]')
box off
skipXTickLabel(gca,2)

drucken('./Plots/tharT.png')


figure
plot(xobs,Hobs,'--','Color','k','LineWidth',1.5)
hold on
for i = 1:length(ceba)
    plot(x,-ceba(i).proc.Cyc.(selSeas).Hflx,'Color',newColor(colNum(i)),'LineWidth',2.5)
end

set(gca,'FontSize',22,'FontName','Verdana','xtick',xTicks,'xlim',[xS xTicks(1)+nDay],'xticklabel',xLabels)
ylabel('H [W/m^2]')
box off
skipXTickLabel(gca,2)
drucken('./Plots/tharH.png')

figure
plot(xobs,LEobs,'--','Color','k','LineWidth',1.5)
hold on
for i = 1:length(ceba)
    plot(x,-ceba(i).proc.Cyc.(selSeas).Lflx,'Color',newColor(colNum(i)),'LineWidth',2.5)
end

set(gca,'FontSize',22,'FontName','Verdana','xtick',xTicks,'xlim',[xS xTicks(1)+nDay],'xticklabel',xLabels)
ylabel('LE [W/m^2]')
box off
skipXTickLabel(gca,2)
drucken('./Plots/tharLE.png')


