%Version 4 - Allsites_flexible_fillerV4
%Uncomment the site you want to run
%You can either hardwire the energy budget closure, or have the program
%calculate it by setting Energyclose = 0

%Version 2 - In general, a simple (conservative) approach)
% R calculated by extrapolating u*>0.3 periods to darkness using nonlinear
% fit.  Includes cleaning to remove outliers and uses regression info from
% previous period if results are unreasonable (really small or large R)
% 
% GEE calcuted from same nonlinear fit, based on sobtracting R from NEE or 
% calculating the GEE fit by dropping outthe first (R) term
%
% NEE calculated similar to GEE
%
% Et claculated from linear regression between Et and K
% Et set to zero at night - this was needed since regrsssion yielded some
% high night Ets in summer




% Read in data
% 
% load('C:\Analysis\California\Tower data\Sage_v3_6.mat');
% Startyear = 2006;
% %Energyclose = 0.90;
% Energyclose = 0; 
% Good = DATA(:,1) >= 200;
% DATA = DATA(Good,:);
% 
% load('C:\Analysis\California\Tower data\Grass_v3_6.mat');
% Startyear = 2006;
% %Energyclose = 0.90;
% Energyclose = 0; 
% Good = DATA(:,1) >= 200;
% DATA = DATA(Good,:);
% 
% load('C:\Analysis\California\Tower data\JamesRes_v3_2.mat');
% Startyear = 2006;
% %Energyclose = 0.73;
% Energyclose = 0; 

% load('C:\Analysis\California\Tower data\Pinyon_v3_6.mat');
% Startyear = 2006;
% %Energyclose = 0.93;
% Energyclose = 0; 
% 
% load('C:\Analysis\California\Tower data\PinyonBurn_v3_6.mat');
% Startyear = 2006;
% %Energyclose = 0.75;
% Energyclose = 0; 
% 
% load('C:\Analysis\California\Tower data\LowDesert_v3_2.mat');
% Startyear = 2006;
% Energyclose = 0.67; 
% 
% load('C:\Analysis\California\Tower data\SJER_v3_6.mat');
% Startyear = 2010;
% %Energyclose = 0.80;
% Energyclose = 0; 
% % 
% load('C:\Analysis\California\Tower data\Soaproot_v3_6.mat');
% Startyear = 2010;
% %Energyclose = 0.77;
% Energyclose = 0; 
% % 
% load('C:\Analysis\California\Tower data\P301_v3_6.mat');
% Startyear = 2008;
% %Energyclose = 0.78;  
% Energyclose = 0; 
% 
load('C:\Analysis\California\Tower data\Shorthair_v3_6.mat');
Startyear = 2010;
%Energyclose = 0.78;  %Estimated
Energyclose = 0.887;  %Estimated using Rnet vs G regresion from P301

% Set some filters
Ustarfilter = 0.2;
dimthresh = 200;
brightthresh = 400;

%Calculate energy budget closure if needed
if Energyclose == 0;
    H = DATA(:,17);
    LE = DATA(:,18);
    Rnet = DATA(:,15);
    G = DATA(:,16);
    good = isfinite(Rnet + G);
    Gfit = polyfit(Rnet(good),G(good),1);
    Gfit = Rnet * Gfit(1) + Gfit(2); 
    G(~good) = Gfit(~good);
    binwid = 50;
    starter = floor(min(DATA(:,1)));
    lngth = floor(max(DATA(:,1) - min(DATA(:,1)))/50) + 2;
    Closureoutbin = NaN * zeros(lngth,2);
    good = isfinite(H + LE + Rnet + G);
    for i = 1:lngth;
        beginner = starter + (i-1)*binwid;
        ender = beginner + binwid;
        focus = DATA(:,1) >= beginner & DATA(:,1) < ender;
        xxx = Rnet(good & focus) - G(good & focus);
        yyy = H(good & focus) + LE(good & focus);
        Closureoutbin(i,1) = beginner + binwid/2;
        Closureoutbin(i,2) = xxx\yyy;
    end
    bad = Closureoutbin(:,2) == 0;
    threshbad = median(Closureoutbin(~bad,2));
    bad = Closureoutbin(:,2) < threshbad;
    %bad = Closureoutbin(:,2) < 0.55;
    Closureoutbin = Closureoutbin(~bad,:);
    Closures = interp1(Closureoutbin(:,1), Closureoutbin(:,2), DATA(:,1), 'linear', 'extrap');
else
    Closures = zeros(size(DATA(:,1))) + Energyclose;
end


%Creat continous meteorology fields first
%use two filling approaches
%for short missing periods (<=4 missing intervals) just linear extrapolate
%for longer missing intervals fill in the mean of 25 day blocks

%create filled T first
%start by calculating length of missing period
Tconsensus = DATA(:,5);
bad = isnan(Tconsensus);
Tconsensus(bad) = DATA(bad,6);
bad = isnan(Tconsensus);
Tunfill = Tconsensus;
Tfill = Tconsensus;
missing = isnan(Tunfill);
counter = 0; 

for i = 1:length(missing);
    if missing(i) == 0;
       counter = 0 ;
    else;
       counter = counter + 1;
   end;
cummissing(i) = counter;
end;

%now create liner extrapolate
Tinterp = interp1(DATA(~missing,1),Tunfill(~missing),DATA(:,1),'linear');

%now creat series of mean diel curves for 25-day chnunks
binwid = 25;
start = min(DATA(:,1));
binnum = floor((max(DATA(:,1)- start)/25));
bincenter = (0:(binnum-1))* binwid + start;
Toutbin = zeros(binnum,48);
Hour = 24 * (DATA(:,1) - floor(DATA(:,1)));
good = isfinite(Tunfill);

for i = 1:binnum;
    focusday = (DATA(:,1) > (bincenter(i) - 12.5)) & (DATA(:,1) <= (bincenter(i) + 12.5));
    for j = 1:48;
        focushour = Hour > j/2 - 0.75 & Hour < j/2 - 0.25;
        Toutbin(i,j) = mean(Tunfill(good&focusday&focushour));
    end
end
Toutbin = [bincenter' Toutbin];

%Now do the filling
for i = 1:length(missing);
    if missing(i);
        if cummissing(i) < 5;
            Tfill(i) = Tinterp(i);
        else
            rownumb = floor((max(DATA(i,1)- start+12.5)/25)) + 1;
            colnumb = round(Hour(i)*2) + 2;
            if rownumb > binnum;
                rownumb = binnum;
            end
            Tfill(i) = Toutbin(rownumb,colnumb);
        end
    end
end


%Now create filled light
%start by calculating length of missing period
Kunfill = DATA(:,13);
missing = isnan(Kunfill);
Kunfill(missing) = DATA(missing,11).*0.52647;
Kfill = Kunfill;
missing = isnan(Kunfill);
counter = 0; 

for i = 1:length(missing);
    if missing(i) == 0;
       counter = 0 ;
    else;
       counter = counter + 1;
   end;
cummissing(i) = counter;
end;

%now create liner extrapolate
Kinterp = interp1(DATA(~missing,1),Kunfill(~missing),DATA(:,1),'linear');

%now creat series of mean diel curves for 25-day chnunks
binwid = 25;
start = min(DATA(:,1));
binnum = floor((max(DATA(:,1)- start)/25));
bincenter = (0:(binnum-1))* binwid + start;
Koutbin = zeros(binnum,48);
Hour = 24 * (DATA(:,1) - floor(DATA(:,1)));
good = isfinite(Kunfill);

for i = 1:binnum;
    focusday = (DATA(:,1) > (bincenter(i) - 12.5)) & (DATA(:,1) <= (bincenter(i) + 12.5));
    for j = 1:48;
        focushour = Hour > j/2 - 0.75 & Hour < j/2 - 0.25;
        Koutbin(i,j) = mean(Kunfill(good&focusday&focushour));
    end
end
Koutbin = [bincenter' Koutbin];

%Now do the filling
for i = 1:length(missing);
    if missing(i);
        if cummissing(i) < 5;
            Kfill(i) = Kinterp(i);
        else
            rownumb = floor((max(DATA(i,1)- start+12.5)/25)) + 1;
            colnumb = round(Hour(i)*2) + 2;
            if rownumb > binnum;
                rownumb = binnum;
            end
            Kfill(i) = Koutbin(rownumb,colnumb);
        end
    end
end



%Now create filled Rnet
%start by calculating length of missing period
Rnetunfill = DATA(:,15);
missing = isnan(Rnetunfill);
Rnetunfill(missing) = Kfill(missing)*0.70582-43.776;
Rnetfill = Rnetunfill;


% %Now do CO2 flux filling 
% %%% Now calculate Resp using light curves
NEEunfill = DATA(:,19);
calm = DATA(:,2) < Ustarfilter & isfinite(DATA(:,2));
NEEunfill(calm) = NaN;
bad = DATA(:,19) > 10 | DATA(:,19) < -35;
NEEunfill(bad) = NaN;

NEEfit = DATA(:,19) * NaN;
Respfit = DATA(:,19) * NaN;
GEEfit = DATA(:,19) * NaN;

%now creat series of light curves for 10-day chnunks and calculate the flux
%at each irradaince

binwid = 10;
start = binwid/2+min(DATA(:,1));
binnum = floor((max(DATA(:,1)- start)/binwid));
bincenter = (0:(binnum-1))* binwid + start + binwid/2;
Toutbin = zeros(binnum,48);
Hour = 24 * (DATA(:,1) - floor(DATA(:,1)));
good = isfinite(Kunfill);
GEEfit = NEEunfill * NaN;

parfitdimold = [-0.001  3];
parfitmediumold = [-0.001  3];
parfitbrightold = [-0.001  3];

for i = 1:binnum;
    focusday = (DATA(:,1) > (bincenter(i) - 5)) & (DATA(:,1) <= (bincenter(i) + 5));
    par_day = Kfill(focusday & isfinite(NEEunfill));
    nee_day = NEEunfill(focusday & isfinite(NEEunfill));
  
    dim = Kfill < dimthresh;
    medium = Kfill >= dimthresh & Kfill < brightthresh;
    bright = Kfill >= brightthresh;

    par_day = Kfill(focusday & isfinite(NEEunfill) & dim);
    nee_day = NEEunfill(focusday & isfinite(NEEunfill) & dim);
    if length(nee_day) > 5
        parfitdim = polyfit(par_day, nee_day, 1);
        parfitdimold = parfitdim;
    else 
        parfitdim = parfitdimold;   
    end
        
    par_day = Kfill(focusday & isfinite(NEEunfill) & medium);
    nee_day = NEEunfill(focusday & isfinite(NEEunfill) & medium);
    if length(nee_day) > 5
        parfitmedium = polyfit(par_day, nee_day, 1);
        parfitmediumold = parfitmedium;
    else 
        parfitmedium = parfitmediumold;   
    end
    
    par_day = Kfill(focusday & isfinite(NEEunfill) & bright);
    nee_day = NEEunfill(focusday & isfinite(NEEunfill) & bright);
    if length(nee_day) > 5
        parfitbright = polyfit(par_day, nee_day, 1);
        parfitbrightold = parfitbright;
    else 
        parfitbright = parfitmediumold;   
    end
    
    Respfit(focusday) = parfitdim(2);        
    NEEfit(focusday & dim) = parfitdim(2) + (parfitdim(1)*Kfill(focusday & dim));
    NEEfit(focusday & medium) = parfitmedium(2) + (parfitmedium(1)*Kfill(focusday & medium));
    NEEfit(focusday & bright) = parfitbright(2) + (parfitbright(1)*Kfill(focusday & bright));
    GEEfit(focusday) =  NEEfit(focusday) - Respfit(focusday);
    
    par_day = Kfill(focusday);
    nee_day = NEEunfill(focusday);
    
    %     plot(par_day,nee_day,'.')
    %     i
    %     b
    %     pause
    
end

% 
% 
% %Now calculate annual respiration
missing = isnan(Respfit);
Respfit(missing) = 0;
Respfill = Respfit;

% %%% Now calculate GEE
GEEunfill = DATA(:,19) - Respfill;
calm = DATA(:,2) < Ustarfilter & isfinite(DATA(:,2));
GEEunfill(calm) = NaN;
light = Kfill > 2;
GEEunfill(~light) = 0;
missing = isnan(GEEunfill);
GEEfill = GEEunfill;
GEEfill(missing) = GEEfit(missing);
missing = isnan(GEEfill);
GEEfill(missing) = 0;

% %Finally calculate filled NEE and NEP
NEEunfill = DATA(:,19);
calm = DATA(:,2) < Ustarfilter & isfinite(DATA(:,2));
NEEunfill(calm) = NaN;
NEEfill = NEEunfill;
missing = isnan(NEEfill);
NEEfill(missing) = NEEfit(missing);
missing = isnan(NEEfill);
NEEfill(missing) = 0;



%%%%%%
%%% Now do water and energy
%%%%%%
%%% Now calculate E
Eunfill = DATA(:,20);
calm = DATA(:,2) < Ustarfilter & isfinite(DATA(:,2));
Eunfill(calm) = NaN;
Efit = Eunfill * NaN;

%now creat series of light curves for 10-day chnunks and calculate the flux
%at each irradaince

binwid = 10;
start = binwid/2+min(DATA(:,1));
binnum = floor((max(DATA(:,1)- start)/binwid));
bincenter = (0:(binnum-1))* binwid + start + binwid/2;
Toutbin = zeros(binnum,48);
Hour = 24 * (DATA(:,1) - floor(DATA(:,1)));
good = isfinite(Kunfill);
Efit = Eunfill * NaN;
parfitold = [0.8 0.006];

for i = 1:binnum;
    focusday = (DATA(:,1) > (bincenter(i) - 5)) & (DATA(:,1) <= (bincenter(i) + 5));
    par_day = Kfill(focusday & isfinite(Eunfill));
    nee_day = Eunfill(focusday & isfinite(Eunfill));
    
    X = [ones(size(par_day)) par_day];
    y = nee_day;
    
    if length(nee_day) > 20
        b = regress(y,X);
    else
        b = parfitold;
    end
    parfitold = b;
    
    Efit(focusday) = b(1) + b(2)*Kfill(focusday);
    
    %     plot(par_day,nee_day,'.')
    %     i
    %     b
    %     pause
    
end

%Now do the E filling
missing = isnan(Eunfill);
Efill = Eunfill;
Efill(missing) = Efit(missing);
missing = isnan(Efill);
Efill(missing) = 0;


%%%%%%
%%%%%%
%%% Now calculate H
Hunfill = DATA(:,17);

calm = DATA(:,2) < Ustarfilter & isfinite(DATA(:,2));
Hunfill(calm) = NaN;

Hfit = Hunfill * NaN;

%now creat series of light curves for 10-day chnunks and calculate the flux
%at each irradaince

binwid = 10;
start = binwid/2+min(DATA(:,1));
binnum = floor((max(DATA(:,1)- start)/binwid));
bincenter = (0:(binnum-1))* binwid + start + binwid/2;
Toutbin = zeros(binnum,48);
Hour = 24 * (DATA(:,1) - floor(DATA(:,1)));
good = isfinite(Kunfill);
Hfit = Hunfill * NaN;
parfitold = [-20 0.6];

for i = 1:binnum;
    focusday = (DATA(:,1) > (bincenter(i) - 5)) & (DATA(:,1) <= (bincenter(i) + 5));
    par_day = Kfill(focusday & isfinite(Hunfill));
    nee_day = Hunfill(focusday & isfinite(Hunfill));
    
    X = [ones(size(par_day)) par_day];
    y = nee_day;
    
    if length(nee_day) > 10
        b = regress(y,X);
    else
        b = parfitold;
    end
    parfitold = b;
    
    Hfit(focusday) = b(1) + b(2)*Kfill(focusday);
    
%         plot(par_day,nee_day,'.')
%         i
%         b
%         pause
    
end

%Now do the H filling
missing = isnan(Hunfill);
Hfill = Hunfill;
Hfill(missing) = Hfit(missing);
missing = isnan(Hfill);
Hfill(missing) = 0;

%%%%%%%%
%%%%%%%%
%%%%%%%%
%Now calculate calculate all of the sums and annual fluxes
cumH = cumsum(Hfill./Closures)*(1800)/(1000*1000);
cumE = cumsum(Efill./Closures)*(18.02*1800)/(1000*1000);
cumNEE = cumsum(NEEfill./Closures)/(2.31*2);
cumGEE = cumsum(GEEfill./Closures)/(2.31*2);
cumResp = cumsum(Respfill./Closures)/(2.31*2);

%Now calculate annual sums
wateryrstart = datenum(2006,10,1) - datenum(2006,1,0);
%wateryrstart = datenum(2007,1,0) - datenum(2006,1,0);   %Use this for
%calander years

if Startyear <= 2006;
is2006 = DATA(:,1) > wateryrstart - 0.01 & DATA(:,1) < wateryrstart + 0.01;
Resp2006 = cumResp(is2006);
GEE2006 = cumGEE(is2006);
NEE2006 = cumNEE(is2006);
E2006 = cumE(is2006);
H2006 = cumH(is2006);

else
Resp2006 = 0;
GEE2006 = 0;
NEE2006 = 0;
E2006 = 0;
H2006 = 0;
end

if Startyear <= 2007;
is2007 = DATA(:,1) > wateryrstart - 0.01 + 1*365.25 & DATA(:,1) < wateryrstart + 0.01 + 1*365.25;
Resp2007 = cumResp(is2007) - Resp2006;
GEE2007 = cumGEE(is2007) - GEE2006;
NEE2007 = cumNEE(is2007) - NEE2006;
E2007 = cumE(is2007) - E2006;
H2007 = cumH(is2007) - H2006;

else
Resp2007 = 0;
GEE2007 = 0;
NEE2007 = 0;
E2007 = 0;
H2007 = 0;
end

if Startyear <= 2008;
is2008 = DATA(:,1) > wateryrstart - 0.01 + 2*365.25 & DATA(:,1) < wateryrstart + 0.01 + 2*365.25;
Resp2008 = cumResp(is2008) - Resp2006 - Resp2007;
GEE2008 = cumGEE(is2008) - GEE2006 - GEE2007;
NEE2008 = cumNEE(is2008) - NEE2006 - NEE2007;
E2008 = cumE(is2008) - E2006 - E2007;
H2008 = cumH(is2008) - H2006 - H2007;

else
Resp2008 = 0;
GEE2008 = 0;
NEE2008 = 0;
E2008 = 0;
H2008 = 0;
end

if Startyear <= 2009;
is2009 = DATA(:,1) > wateryrstart - 0.01 + 3*365.25 & DATA(:,1) < wateryrstart + 0.01 + 3*365.25;
Resp2009 = cumResp(is2009) - Resp2006 - Resp2007 - Resp2008;
GEE2009 = cumGEE(is2009) - GEE2006 - GEE2007 - GEE2008;
NEE2009 = cumNEE(is2009) - NEE2006 - NEE2007 - NEE2008;
E2009 = cumE(is2009) - E2006 - E2007 - E2008;
H2009 = cumH(is2009) - H2006 - H2007 - H2008;

else
Resp2009 = 0;
GEE2009 = 0;
NEE2009 = 0;
E2009 = 0;
H2009 = 0;
end

if Startyear <= 2010;
is2010 = DATA(:,1) >  wateryrstart - 0.01 + 4*365.25 & DATA(:,1) < wateryrstart + 0.01 + 4*365.25;
Resp2010 = cumResp(is2010) - Resp2006 - Resp2007 - Resp2008 - Resp2009;
GEE2010 = cumGEE(is2010) - GEE2006 - GEE2007 - GEE2008 - GEE2009;
NEE2010 = cumNEE(is2010) - NEE2006 - NEE2007 - NEE2008 - NEE2009;
E2010 = cumE(is2010) - E2006 - E2007 - E2008 - E2009;
H2010 = cumH(is2010) - H2006 - H2007 - H2008 - H2009;

else
Resp2010 = 0;
GEE2010 = 0;
NEE2010 = 0;
E2010 = 0;
H2010 = 0;
end

if Startyear <= 2011;
is2011 = DATA(:,1) >  wateryrstart - 0.01 + 5*365.25 & DATA(:,1) < wateryrstart + 0.01 + 5*365.25;
Resp2011 = cumResp(is2011) - Resp2006 - Resp2007 - Resp2008 - Resp2009 - Resp2010;
GEE2011 = cumGEE(is2011) - GEE2006 - GEE2007 - GEE2008 - GEE2009 - GEE2010;
NEE2011 = cumNEE(is2011) - NEE2006 - NEE2007 - NEE2008 - NEE2009 - NEE2010;
E2011 = cumE(is2011) - E2006 - E2007 - E2008 - E2009 - E2010;
H2011 = cumH(is2011) - H2006 - H2007 - H2008 - H2009 - H2010;

else
Resp2011 = 0;
GEE2011 = 0;
NEE2011 = 0;
E2011 = 0;
H2011 = 0;
end

if Startyear <= 2012;
is2012 = DATA(:,1) >  wateryrstart - 0.01 + 6*365.25 & DATA(:,1) < wateryrstart + 0.01 + 6*365.25;
Resp2012 = cumResp(is2012) - Resp2006 - Resp2007 - Resp2008 - Resp2009 - Resp2010 - Resp2011;
GEE2012 = cumGEE(is2012) - GEE2006 - GEE2007 - GEE2008 - GEE2009 - GEE2010 - GEE2011;
NEE2012 = cumNEE(is2012) - NEE2006 - NEE2007 - NEE2008 - NEE2009 - NEE2010 - NEE2011;
E2012 = cumE(is2012) - E2006 - E2007 - E2008 - E2009 - E2010 - E2011;
H2012 = cumH(is2012) - H2006 - H2007 - H2008 - H2009 - H2010 - H2011;

else
Resp2012 = 0;
GEE2012 = 0;
NEE2012 = 0;
E2012 = 0;
H2012 = 0;
end

if Startyear <= 2013;
is2013 = DATA(:,1) >  wateryrstart - 0.01 + 7*365.25 & DATA(:,1) < wateryrstart + 0.01 + 7*365.25;
Resp2013 = cumResp(is2013) - Resp2006 - Resp2007 - Resp2008 - Resp2009 - Resp2010 - Resp2011 - Resp2012;
GEE2013 = cumGEE(is2013) - GEE2006 - GEE2007 - GEE2008 - GEE2009 - GEE2010 - GEE2011 - GEE2012;
NEE2013 = cumNEE(is2013) - NEE2006 - NEE2007 - NEE2008 - NEE2009 - NEE2010 - NEE2011 - NEE2012;
E2013 = cumE(is2013) - E2006 - E2007 - E2008 - E2009 - E2010 - E2011 - E2012;
H2013 = cumH(is2013) - H2006 - H2007 - H2008 - H2009 - H2010 - H2011 - H2012;

else
Resp2013 = 0;
GEE2013 = 0;
NEE2013 = 0;
E2013 = 0;
H2013 = 0;
end

if Startyear <= 2014;
is2014 = DATA(:,1) >  wateryrstart - 0.01 + 8*365.25 & DATA(:,1) < wateryrstart + 0.01 + 8*365.25;
Resp2014 = cumResp(is2014) - Resp2006 - Resp2007 - Resp2008 - Resp2009 - Resp2010 - Resp2011 - Resp2012 - Resp2013;
GEE2014 = cumGEE(is2014) - GEE2006 - GEE2007 - GEE2008 - GEE2009 - GEE2010 - GEE2011 - GEE2012 - GEE2013;
NEE2014 = cumNEE(is2014) - NEE2006 - NEE2007 - NEE2008 - NEE2009 - NEE2010 - NEE2011 - NEE2012 - NEE2013;
E2014 = cumE(is2014) - E2006 - E2007 - E2008 - E2009 - E2010 - E2011 - E2012 - E2013;
H2014 = cumH(is2014) - H2006 - H2007 - H2008 - H2009 - H2010 - H2011 - H2012 - H2013;

else
Resp2014 = 0;
GEE2014 = 0;
NEE2014 = 0;
E2014 = 0;
H2014 = 0;
end


if Startyear <= 2015;
is2015 = DATA(:,1) >  wateryrstart - 0.01 + 9*365.25 & DATA(:,1) < wateryrstart + 0.01 + 9*365.25;
Resp2015 = cumResp(is2015) - Resp2006 - Resp2007 - Resp2008 - Resp2009 - Resp2010 - Resp2011 - Resp2012 - Resp2013 - Resp2014;
GEE2015 = cumGEE(is2015) - GEE2006 - GEE2007 - GEE2008 - GEE2009 - GEE2010 - GEE2011 - GEE2012 - GEE2013 - GEE2014;
NEE2015 = cumNEE(is2015) - NEE2006 - NEE2007 - NEE2008 - NEE2009 - NEE2010 - NEE2011 - NEE2012 - NEE2013 - NEE2014;
E2015 = cumE(is2015) - E2006 - E2007 - E2008 - E2009 - E2010 - E2011 - E2012 - E2013 - E2014;
H2015 = cumH(is2015) - H2006 - H2007 - H2008 - H2009 - H2010 - H2011 - H2012 - H2013 - H2014;

else
Resp2015 = 0;
GEE2015 = 0;
NEE2015 = 0;
E2015 = 0;
H2015 = 0;
end


if Startyear <= 2016;
is2016 = DATA(:,1) >  wateryrstart - 0.01 + 10*365.25 & DATA(:,1) < wateryrstart + 0.01 + 10*365.25;
Resp2016 = cumResp(is2016) - Resp2006 - Resp2007 - Resp2008 - Resp2009 - Resp2010 - Resp2011 - Resp2012 - Resp2013 - Resp2014 - Resp2015;
GEE2016 = cumGEE(is2016) - GEE2006 - GEE2007 - GEE2008 - GEE2009 - GEE2010 - GEE2011 - GEE2012 - GEE2013 - GEE2014 - GEE2015;
NEE2016 = cumNEE(is2016) - NEE2006 - NEE2007 - NEE2008 - NEE2009 - NEE2010 - NEE2011 - NEE2012 - NEE2013 - NEE2014 - NEE2015;
E2016 = cumE(is2016) - E2006 - E2007 - E2008 - E2009 - E2010 - E2011 - E2012 - E2013 - E2014 - E2015;
H2016 = cumH(is2016) - H2006 - H2007 - H2008 - H2009 - H2010 - H2011 - H2012 - H2013 - H2014 - H2015;

else
Resp2015 = 0;
GEE2015 = 0;
NEE2015 = 0;
E2015 = 0;
H2015 = 0;
end


R = [Resp2007 Resp2008 Resp2009 Resp2010 Resp2011 Resp2012 Resp2013 Resp2014 Resp2015 Resp2016];
GEE = [GEE2007 GEE2008 GEE2009 GEE2010 GEE2011 GEE2012 GEE2013 GEE2014 GEE2015 GEE2016];
NEE = [NEE2007 NEE2008 NEE2009 NEE2010 NEE2011 NEE2012 NEE2013 NEE2014 NEE2015 NEE2016];
E = [E2007 E2008 E2009 E2010 E2011 E2012 E2013 E2014 E2015 E2016];
H = [H2007 H2008 H2009 H2010 H2011 H2012 H2013 H2014 H2015 H2016];
AnnC = [R' GEE' NEE'];
AnnFluxes = [R' GEE' NEE' E' H']


%Use this code if you want to write filled time series
Filledup = [DATA(:,1) Kfill Rnetfill Tfill NEEfill GEEfill Respfill Efill Hfill];
% save('C:\Analysis\California\Tower data\Sage_Filled.mat','Filledup');
% save('C:\Analysis\California\Tower data\Grass_Filled.mat','Filledup');
% save('C:\Analysis\California\Tower data\Pinyon_Filled.mat','Filledup');
% save('C:\Analysis\California\Tower data\PFburn_Filled.mat','Filledup');
% save('C:\Analysis\California\Tower data\Desert_Filled.mat','Filledup');
% save('C:\Analysis\California\Tower data\James_Filled.mat','Filledup');
%save('C:\Analysis\California\Tower data\SJER_Filled.mat','Filledup');
%save('C:\Analysis\California\Tower data\Soaproot_Filled.mat','Filledup');
%save('C:\Analysis\California\Tower data\P301_Filled.mat','Filledup');
%save('C:\Analysis\California\Tower data\Shorthair_Filled.mat','Filledup');

% 
% %Use this code to calculate daily sums
hoursrunning = round(48 * (DATA(:,1) - floor(DATA(:,1))));
cumulativeET = cumE;
cumulativeGEE = cumGEE;
midnight = hoursrunning == 14;
Date = DATA(midnight,1);
Date = Date(2:length(Date));
Date = Date + datenum(2006,1,0) - datenum(1900,1,0);
%Date = datestr(Date);
ET = cumulativeET(midnight);
ET = ET(2:length(ET)) - ET(1:(length(ET)-1));
DailyET = [Date ET];
GEE = cumulativeGEE(midnight);
GEE = GEE(2:length(GEE)) - GEE(1:(length(GEE)-1));
DailyGEE = [Date GEE];
% 
% dlmwrite('C:\Analysis\California\Tower data\P301_DailyET2.csv',DailyET,'delimiter',',','precision', 7);
% dlmwrite('C:\Analysis\California\Tower data\Soaproot_DailyET2.csv',DailyET,'delimiter',',','precision', 7);
 dlmwrite('C:\Analysis\California\Tower data\SJER_DailyET2.csv',DailyET,'delimiter',',','precision', 7);
% dlmwrite('C:\Analysis\California\Tower data\P301_DailyGEE2.csv',DailyGEE,'delimiter',',','precision', 7);
% dlmwrite('C:\Analysis\California\Tower data\Soaproot_DailyGEE2.csv',DailyGEE,'delimiter',',','precision', 7);
 dlmwrite('C:\Analysis\California\Tower data\SJER_DailyGEE2.csv',DailyGEE,'delimiter',',','precision', 7);
%     

% 
% %Use this code to calculate monthly sums
% hoursrunning = round(48 * (DATA(:,1) - floor(DATA(:,1))));
% midnight = hoursrunning == 14;
% dayofmonth = day(DATA(:,1));
% pulldate = dayofmonth == 1 & midnight;
% 
% cumulativeET = cumE;
% cumulativeGEE = cumGEE;
% Date = DATA(pulldate,1);
% numdays = Date(2:length(Date)) - Date(1:(length(Date)-1));
% Date = Date(2:length(Date));
% Date = Date + datenum(2006,1,0) - datenum(1900,1,0);
% ET = cumulativeET(pulldate);
% ET = ET(2:length(ET)) - ET(1:(length(ET)-1));
% ET = ET./numdays;
% MonthlyET = [Date ET];
% GEE = cumulativeGEE(pulldate);
% GEE = GEE(2:length(GEE)) - GEE(1:(length(GEE)-1));
% GEE = GEE./numdays;
% MonthlyGEE = [Date GEE];
% % 
% % dlmwrite('C:\Analysis\California\Tower data\P301_MonthlyET2.csv',MonthlyET,'delimiter',',','precision', 7);
% % dlmwrite('C:\Analysis\California\Tower data\Soaproot_MonthlyET2.csv',MonthlyET,'delimiter',',','precision', 7);
% % dlmwrite('C:\Analysis\California\Tower data\SJER_MonthlyET2.csv',MonthlyET,'delimiter',',','precision', 7);
% % dlmwrite('C:\Analysis\California\Tower data\P301_MonthlyGEE2.csv',MonthlyGEE,'delimiter',',','precision', 7);
% % dlmwrite('C:\Analysis\California\Tower data\Soaproot_MonthlyGEE2.csv',MonthlyGEE,'delimiter',',','precision', 7);
% % dlmwrite('C:\Analysis\California\Tower data\SJER_MonthlyGEE2.csv',MonthlyGEE,'delimiter',',','precision', 7);
% %     
% 
% 
% 



% %Use this code if you want to create coarser time series
% % 
% %Summing
% %binwid = 10; %length of averaging period in days
% start = binwid/2+min(DATA(:,1));
% binnum = floor((max(DATA(:,1)- start)/binwid));
% bincenter = (0:(binnum-1))* binwid + start + binwid/2;
% Sageoutbin = zeros(binnum,6);
% Sageoutbin(:,1) = bincenter;
% for i = 1:binnum;
%     focusday = (DATA(:,1) > (bincenter(i) - 0.5)) & (DATA(:,1) <= (bincenter(i) + 0.5));
%     Sageoutbin(i,2) = sum(FilledSage(focusday,5))./(2.31*2);%*0.90 accoutns for energy budget closure
%     Sageoutbin(i,3) = sum(FilledSage(focusday,6))./(2.31*2);%*0.90 accoutns for energy budget closure
%     Sageoutbin(i,4) = sum(FilledSage(focusday,7))./(2.31*2);%*0.90 accoutns for energy budget closure
%     Sageoutbin(i,5) = sum(FilledSage(focusday,8))*(18.02*1800)/(1000*1000);%*0.90 accoutns for energy budget closure
%     Sageoutbin(i,6) = sum(FilledSage(focusday,9))*(1800)/(1000*1000)./0.90;%*0.90 accoutns for energy budget closure
% end
% 
% % 

%Use this to create and write time series of daily soil water
hour = 24*(DATA(:,1) - floor(DATA(:,1)));
focus = hour == 0;
soilwater = horzcat(DATA(focus,1),DATA(focus,33:38));

dlmwrite('C:\Analysis\California\Tower data\SJER_dailysoilmatric.csv',soilwater,'delimiter',',','precision', 7);
%  dlmwrite('C:\Analysis\California\Tower data\Soaproot_dailysoilmatric.csv',soilwater,'delimiter',',','precision', 7);
%  dlmwrite('C:\Analysis\California\Tower data\P301_dailysoilmatric.csv',soilwater,'delimiter',',','precision', 7);

