% load DL1 climate and flux data for 2004-2018
load('DL1_meas_2004_2017_daily.mat') 
%%
% CLASSIC results
% compare simulations using different PFTs against daily observations
pth_out = 'D:\CLASSIC_output\outputfiles_shrub_2004_2017\';
pth_out2 = 'D:\CLASSIC_output\outputfiles_grass_2004_2017\';
pth_out3 = 'D:\CLASSIC_output\outputfiles_trees_2004_2017\';
pth_out4 = 'D:\CLASSIC_output\outputfiles_shrub_2004_2017_original_beta\';

albs_daily  = squeeze(ncread([pth_out 'albs_daily.nc'],'albs')); %albedo
tas_daily  = squeeze(ncread([pth_out 'tas_daily.nc'],'tas')); %surface T
tsl_daily  = squeeze(ncread([pth_out 'tsl_daily.nc'],'tsl')); %{\itT}_s
mrsl_daily = squeeze(ncread([pth_out 'mrsll_daily.nc'],'mrsll')); %vol water content (liquid)
mrfl_daily = squeeze(ncread([pth_out 'mrsfl_daily.nc'],'mrsfl')); %vol water content (frozen)
rls_daily  = squeeze(ncread([pth_out 'rls_daily.nc'],'rls')); %net longwave
rss_daily  = squeeze(ncread([pth_out 'rss_daily.nc'],'rss')); %net shortwave
hfss_daily  = squeeze(ncread([pth_out 'hfss_daily.nc'],'hfss')); %H
hfls_daily  = squeeze(ncread([pth_out 'hfls_daily.nc'],'hfls')); %LE
nep_daily  = squeeze(ncread([pth_out 'nep_daily.nc'],'nep')); %NEP (kg C m-2 s-1)
npp_daily  = squeeze(ncread([pth_out 'npp_daily.nc'],'npp')); %NPP (kg C m-2 s-1)
gpp_daily  = squeeze(ncread([pth_out 'gpp_daily.nc'],'gpp')); %GPP (kg C m-2 s-1)
ra_daily  = squeeze(ncread([pth_out 'ra_daily.nc'],'ra')); %Ra (kg C m-2 s-1)
rh_daily  = squeeze(ncread([pth_out 'rh_daily.nc'],'rh')); %Rh (kg C m-2 s-1)
nee_daily = -nep_daily; %NEE  (kg C m-2 s-1)

mrro_daily = squeeze(ncread([pth_out 'mrro_daily.nc'],'mrro')); %total runoff from soil
mrros_daily = squeeze(ncread([pth_out 'mrros_daily.nc'],'mrros')); %overland flow from top of soil column
mrrob_daily = squeeze(ncread([pth_out 'mrrob_daily.nc'],'mrrob')); %base flow from bottom
mrroc_daily = squeeze(ncread([pth_out 'mrroc_daily.nc'],'mrroc')); %Liquid/frozen water runoff from vegetation
mrroi_daily = squeeze(ncread([pth_out 'mrroi_daily.nc'],'mrroi')); %interflow
mrron_daily = squeeze(ncread([pth_out 'mrron_daily.nc'],'mrron')); %Liquid runoff from snowpack

zpond_daily = squeeze(ncread([pth_out 'zpond_daily.nc'],'zpond')); %Depth of ponded water [m]
evap_daily = squeeze(ncread([pth_out 'evspsbl_daily.nc'],'evspsbl')); 
groundevap_daily = squeeze(ncread([pth_out 'evspsblsoi_daily.nc'],'evspsblsoi')); 
canopyevap_daily = squeeze(ncread([pth_out 'evspsblveg_daily.nc'],'evspsblveg')); 
tran_daily     = squeeze(ncread([pth_out 'tran_daily.nc'],'tran')); %transpiration (kg m-2 s-1)
pr_daily     = squeeze(ncread([pth_out 'pr_daily.nc'],'pr')); %precipitation (kg m-2 s-1)
snw_daily     = squeeze(ncread([pth_out 'snw_daily.nc'],'snw')); %surface snow amount (kg m-2)
wsnw_daily     = squeeze(ncread([pth_out 'wsnw_daily.nc'],'wsnw')); %Liquid water content of snowpack (kg m-2)
snwdens_daily     = squeeze(ncread([pth_out 'snwdens_daily.nc'],'snwdens')); %density of snow (kg m-3)
tsn_daily  = squeeze(ncread([pth_out 'tsn_daily.nc'],'tsn')); %snowpack T
zsnw_daily = snw_daily./snwdens_daily; % snow depth (m)
gleafmas_daily = squeeze(ncread([pth_out 'gleafmas_daily.nc'],'gleafmas')); %Green leaf mass (kg C m-2)
bleafmas_daily = squeeze(ncread([pth_out 'bleafmas_daily.nc'],'bleafmas')); %Brown leaf mass (kg C m-2)
stemmass_daily = squeeze(ncread([pth_out 'stemmass_daily.nc'],'stemmass')); %Stem mass (kg C m-2)
rootmass_daily = squeeze(ncread([pth_out 'cRoot_daily.nc'],'cRoot')); %Root mass (kg C m-2)
%rootdpth_daily = squeeze(ncread([pth_out 'rootdpth_daily.nc'],'rootdpth')); %Root depth
%lfstatus_daily = squeeze(ncread([pth_out 'lfstatus_daily.nc'],'lfstatus')); %Leaf status (phenology)

tas_daily2  = squeeze(ncread([pth_out2 'tas_daily.nc'],'tas')); %air T
tsl_daily2  = squeeze(ncread([pth_out2 'tsl_daily.nc'],'tsl')); %{\itT}_s
mrsl_daily2 = squeeze(ncread([pth_out2 'mrsll_daily.nc'],'mrsll')); %vol water content (liquid)
mrfl_daily2 = squeeze(ncread([pth_out2 'mrsfl_daily.nc'],'mrsfl')); %vol water content (frozen)
rls_daily2  = squeeze(ncread([pth_out2 'rls_daily.nc'],'rls')); %net longwave
rss_daily2  = squeeze(ncread([pth_out2 'rss_daily.nc'],'rss')); %net shortwave
hfss_daily2  = squeeze(ncread([pth_out2 'hfss_daily.nc'],'hfss')); %H
hfls_daily2  = squeeze(ncread([pth_out2 'hfls_daily.nc'],'hfls')); %LE
albs_daily2  = squeeze(ncread([pth_out2 'albs_daily.nc'],'albs')); %albedo
nep_daily2 = squeeze(ncread([pth_out2 'nep_daily.nc'],'nep')); %NEP
npp_daily2  = squeeze(ncread([pth_out2 'npp_daily.nc'],'npp')); %NPP
gpp_daily2  = squeeze(ncread([pth_out2 'gpp_daily.nc'],'gpp')); %GPP
ra_daily2  = squeeze(ncread([pth_out2 'ra_daily.nc'],'ra')); %Ra
rh_daily2  = squeeze(ncread([pth_out2 'rh_daily.nc'],'rh')); %Rh
nee_daily2 = -nep_daily2; %NEE  (kg C m-2 s-1)

mrro_daily2 = squeeze(ncread([pth_out2 'mrro_daily.nc'],'mrro')); %total runoff from soil
mrros_daily2 = squeeze(ncread([pth_out2 'mrros_daily.nc'],'mrros')); %overland flow from top of soil column
mrrob_daily2 = squeeze(ncread([pth_out2 'mrrob_daily.nc'],'mrrob')); %base flow from bottom
mrroc_daily2 = squeeze(ncread([pth_out2 'mrroc_daily.nc'],'mrroc')); %Liquid/frozen water runoff from vegetation
mrroi_daily2 = squeeze(ncread([pth_out2 'mrroi_daily.nc'],'mrroi')); %interflow
mrron_daily2 = squeeze(ncread([pth_out2 'mrron_daily.nc'],'mrron')); %Liquid runoff from snowpack

zpond_daily2 = squeeze(ncread([pth_out2 'zpond_daily.nc'],'zpond')); %Depth of ponded water [m]
evap_daily2 = squeeze(ncread([pth_out2 'evspsbl_daily.nc'],'evspsbl')); 
groundevap_daily2 = squeeze(ncread([pth_out2 'evspsblsoi_daily.nc'],'evspsblsoi')); 
canopyevap_daily2 = squeeze(ncread([pth_out2 'evspsblveg_daily.nc'],'evspsblveg')); 
tran_daily2     = squeeze(ncread([pth_out2 'tran_daily.nc'],'tran')); %transpiration (kg m-2 s-1)
pr_daily2     = squeeze(ncread([pth_out2 'pr_daily.nc'],'pr')); %precipitation (kg m-2 s-1)
snw_daily2     = squeeze(ncread([pth_out2 'snw_daily.nc'],'snw')); %surface snow amount (kg m-2)
wsnw_daily2     = squeeze(ncread([pth_out2 'wsnw_daily.nc'],'wsnw')); %Liquid water content of snowpack (kg m-2)
snwdens_daily2     = squeeze(ncread([pth_out2 'snwdens_daily.nc'],'snwdens')); %density of snow (kg m-3)
tsn_daily2  = squeeze(ncread([pth_out2 'tsn_daily.nc'],'tsn')); %snowpack T
zsnw_daily2 = snw_daily2./snwdens_daily2; % snow depth (m)
gleafmas_daily2 = squeeze(ncread([pth_out2 'gleafmas_daily.nc'],'gleafmas')); %Green leaf mass (kg C m-2)
bleafmas_daily2 = squeeze(ncread([pth_out2 'bleafmas_daily.nc'],'bleafmas')); %Brown leaf mass (kg C m-2)
stemmass_daily2 = squeeze(ncread([pth_out2 'stemmass_daily.nc'],'stemmass')); %Stem mass (kg C m-2)
rootmass_daily2 = squeeze(ncread([pth_out2 'cRoot_daily.nc'],'cRoot')); %Root mass (kg C m-2)
%rootdpth_daily2 = squeeze(ncread([pth_out2 'rootdpth_daily.nc'],'rootdpth')); %Root depth
t1 = ncread([pth_out 'nep_daily.nc'],'time'); 
t1 = t1 +(tv_av(1)-t1(1));
xData1 = linspace(t1(1),t1(end),180);
[yy, MM, dd, hh, mm] = datevec(t1);

t2 = ncread([pth_out2 'nep_daily.nc'],'time'); 
t2 = t2 +(tv_av(1)-t2(1));
[yy2, MM2, dd2, hh2, mm2] = datevec(t2);
doy2 = doy;


actlyrmax_monthly = squeeze(ncread([pth_out 'actlyrmax_monthly.nc'],'actlyrmax')); %
lai_monthly       = squeeze(ncread([pth_out 'lai_monthly.nc'],'lai')); %
lai_monthly_perpft= squeeze(ncread([pth_out 'lai_monthly_perpft.nc'],'lai')); 
cSoil_monthly     = squeeze(ncread([pth_out 'cSoil_monthly.nc'],'cSoil')); %soil C
cSoil_monthly_perpft= squeeze(ncread([pth_out 'cSoil_monthly_perpft.nc'],'cSoil')); %soil C
t_monthly         = squeeze(ncread([pth_out 'cSoil_monthly.nc'],'time')); %mth
snw_monthly       = squeeze(ncread([pth_out 'snw_monthly.nc'],'snw')); %snow amount
pr_monthly        = squeeze(ncread([pth_out 'pr_monthly.nc'],'pr')); %precipitation (kg m-2 s-1)
evap_monthly = squeeze(ncread([pth_out 'evspsbl_monthly.nc'],'evspsbl'));
groundevap_monthly = squeeze(ncread([pth_out 'evspsblsoi_monthly.nc'],'evspsblsoi')); 
canopyevap_monthly = squeeze(ncread([pth_out 'evspsblveg_monthly.nc'],'evspsblveg')); 
tran_monthly = squeeze(ncread([pth_out 'tran_monthly.nc'],'tran'));
snw_monthly       = squeeze(ncread([pth_out 'snw_monthly.nc'],'snw')); %Mass of snow pack (kg m^{-2} )
tsl_monthly       = squeeze(ncread([pth_out 'tsl_monthly.nc'],'tsl'));
ftable_monthly = squeeze(ncread([pth_out 'ftable_monthly.nc'],'ftable')); % mean depth to frozen water table (m)
cLand_monthly       = squeeze(ncread([pth_out 'cLand_monthly.nc'],'cLand')); % C in all terrestrial C pools (kg C m-2)
cLand_monthly_perpft= squeeze(ncread([pth_out 'cLand_monthly_perpft.nc'],'cLand')); 
cLeaf_monthly     = squeeze(ncread([pth_out 'cLeaf_monthly.nc'],'cLeaf')); 
cLeaf_monthly_perpft= squeeze(ncread([pth_out 'cLeaf_monthly_perpft.nc'],'cLeaf')); 
cLitter_monthly     = squeeze(ncread([pth_out 'cLitter_monthly.nc'],'cLitter')); 
cLitter_monthly_perpft= squeeze(ncread([pth_out 'cLitter_monthly_perpft.nc'],'cLitter')); 
cRoot_monthly     = squeeze(ncread([pth_out 'cRoot_monthly.nc'],'cRoot')); % C in roots
cRoot_monthly_perpft= squeeze(ncread([pth_out 'cRoot_monthly_perpft.nc'],'cRoot'));
cStem_monthly     = squeeze(ncread([pth_out 'cStem_monthly.nc'],'cStem'));
cStem_monthly_perpft = squeeze(ncread([pth_out 'cStem_monthly_perpft.nc'],'cStem'));
cVeg_monthly     = squeeze(ncread([pth_out 'cVeg_monthly.nc'],'cVeg'));
cVeg_monthly_perpft     = squeeze(ncread([pth_out 'cVeg_monthly_perpft.nc'],'cVeg'));
hfg_monthly     = squeeze(ncread([pth_out 'hfg_monthly.nc'],'hfg')); %ground heat flux (W m-2)
gpp_monthly_perpft     = squeeze(ncread([pth_out 'gpp_monthly_perpft.nc'],'gpp'));
rootdpth_monthly_perpft = squeeze(ncread([pth_out 'rootdpth_monthly_perpft.nc'],'rootdpth')); %Root depth

vegHeight_annually = squeeze(ncread([pth_out 'vegHeight_annually.nc'],'vegHeight')); 
vegHeight_annually_perpft = squeeze(ncread([pth_out 'vegHeight_annually_perpft.nc'],['vegHeight'])); 
vegHeight_annually2 = squeeze(ncread([pth_out2 'vegHeight_annually.nc'],'vegHeight')); 
vegHeight_annually_perpft2 = squeeze(ncread([pth_out2 'vegHeight_annually_perpft.nc'],['vegHeight'])); 

actlyrmax_monthly2 = squeeze(ncread([pth_out2 'actlyrmax_monthly.nc'],'actlyrmax')); %
lai_monthly2       = squeeze(ncread([pth_out2 'lai_monthly.nc'],'lai')); %
lai_monthly_perpft2= squeeze(ncread([pth_out2 'lai_monthly_perpft.nc'],'lai')); 
cSoil_monthly2     = squeeze(ncread([pth_out2 'cSoil_monthly.nc'],'cSoil')); %soil C
cSoil_monthly_perpft2 = squeeze(ncread([pth_out2 'cSoil_monthly_perpft.nc'],'cSoil')); %soil C
t_monthly2         = squeeze(ncread([pth_out2 'cSoil_monthly.nc'],'time')); %mth
snw_monthly2       = squeeze(ncread([pth_out2 'snw_monthly.nc'],'snw')); %snow amount
pr_monthly2        = squeeze(ncread([pth_out2 'pr_monthly.nc'],'pr')); %precipitation (kg m-2 s-1)
groundevap_monthly2 = squeeze(ncread([pth_out2 'evspsblsoi_monthly.nc'],'evspsblsoi')); 
canopyevap_monthly2 = squeeze(ncread([pth_out2 'evspsblveg_monthly.nc'],'evspsblveg')); 
snw_monthly2       = squeeze(ncread([pth_out2 'snw_monthly.nc'],'snw')); %Mass of snow pack (kg m^{-2} )
tsl_monthly2       = squeeze(ncread([pth_out2 'tsl_monthly.nc'],'tsl'));
groundevap_monthly2 = squeeze(ncread([pth_out2 'evspsblsoi_monthly.nc'],'evspsblsoi')); 
canopyevap_monthly2 = squeeze(ncread([pth_out2 'evspsblveg_monthly.nc'],'evspsblveg')); 
ftable_monthly2 = squeeze(ncread([pth_out2 'ftable_monthly.nc'],'ftable')); % mean depth to frozen water table (m)
cLand_monthly2       = squeeze(ncread([pth_out2 'cLand_monthly.nc'],'cLand')); % C in all terrestrial C pools (kg C m-2)
cLand_monthly_perpft2= squeeze(ncread([pth_out2 'cLand_monthly_perpft.nc'],'cLand')); 
cLeaf_monthly2     = squeeze(ncread([pth_out2 'cLeaf_monthly.nc'],'cLeaf')); 
cLeaf_monthly_perpft2= squeeze(ncread([pth_out2 'cLeaf_monthly_perpft.nc'],'cLeaf')); 
cLitter_monthly2     = squeeze(ncread([pth_out2 'cLitter_monthly.nc'],'cLitter')); 
cLitter_monthly_perpft2= squeeze(ncread([pth_out2 'cLitter_monthly_perpft.nc'],'cLitter')); 
cRoot_monthly2     = squeeze(ncread([pth_out2 'cRoot_monthly.nc'],'cRoot')); % C in roots
cRoot_monthly_perpft2= squeeze(ncread([pth_out2 'cRoot_monthly_perpft.nc'],'cRoot'));
cStem_monthly2     = squeeze(ncread([pth_out2 'cStem_monthly.nc'],'cStem'));
cStem_monthly_perpft2 = squeeze(ncread([pth_out2 'cStem_monthly_perpft.nc'],'cStem'));
cVeg_monthly2     = squeeze(ncread([pth_out2 'cVeg_monthly.nc'],'cVeg'));
cVeg_monthly_perpft2  = squeeze(ncread([pth_out2 'cVeg_monthly_perpft.nc'],'cVeg'));
hfg_monthly2     = squeeze(ncread([pth_out2 'hfg_monthly.nc'],'hfg')); %ground heat flux (W m-2)
rootdpth_monthly_perpft2 = squeeze(ncread([pth_out2 'rootdpth_monthly_perpft.nc'],'rootdpth')); %Root depth

tsl_daily  = tsl_daily - 273.15; %conversion to degC
tsl_daily2  = tsl_daily2 - 273.15; %conversion to degC

tas_daily3  = squeeze(ncread([pth_out3 'tas_daily.nc'],'tas')); %{\itT}_s
tsl_daily3  = squeeze(ncread([pth_out3 'tsl_daily.nc'],'tsl')); %{\itT}_s
mrsl_daily3 = squeeze(ncread([pth_out3 'mrsll_daily.nc'],'mrsll')); %vol water content (liquid)
mrfl_daily3 = squeeze(ncread([pth_out3 'mrsfl_daily.nc'],'mrsfl')); %vol water content (frozen)
rls_daily3  = squeeze(ncread([pth_out3 'rls_daily.nc'],'rls')); %net longwave
rss_daily3  = squeeze(ncread([pth_out3 'rss_daily.nc'],'rss')); %net shortwave
hfss_daily3  = squeeze(ncread([pth_out3 'hfss_daily.nc'],'hfss')); %H
hfls_daily3  = squeeze(ncread([pth_out3 'hfls_daily.nc'],'hfls')); %LE
albs_daily3  = squeeze(ncread([pth_out3 'albs_daily.nc'],'albs')); %albedo
nep_daily3 = squeeze(ncread([pth_out3 'nep_daily.nc'],'nep')); %NEP
npp_daily3  = squeeze(ncread([pth_out3 'npp_daily.nc'],'npp')); %NPP
gpp_daily3  = squeeze(ncread([pth_out3 'gpp_daily.nc'],'gpp')); %GPP
ra_daily3  = squeeze(ncread([pth_out3 'ra_daily.nc'],'ra')); %Ra
rh_daily3  = squeeze(ncread([pth_out3 'rh_daily.nc'],'rh')); %Rh
nee_daily3 = -nep_daily3; %NEE  (kg C m-2 s-1)

mrro_daily3 = squeeze(ncread([pth_out3 'mrro_daily.nc'],'mrro')); %total runoff from soil
mrros_daily3 = squeeze(ncread([pth_out3 'mrros_daily.nc'],'mrros')); %overland flow from top of soil column
mrrob_daily3 = squeeze(ncread([pth_out3 'mrrob_daily.nc'],'mrrob')); %base flow from bottom
mrroc_daily3 = squeeze(ncread([pth_out3 'mrroc_daily.nc'],'mrroc')); %Liquid/frozen water runoff from vegetation
mrroi_daily3 = squeeze(ncread([pth_out3 'mrroi_daily.nc'],'mrroi')); %interflow
mrron_daily3 = squeeze(ncread([pth_out3 'mrron_daily.nc'],'mrron')); %Liquid runoff from snowpack

zpond_daily3 = squeeze(ncread([pth_out3 'zpond_daily.nc'],'zpond')); %Depth of ponded water [m]
evap_daily3 = squeeze(ncread([pth_out3 'evspsbl_daily.nc'],'evspsbl')); 
groundevap_daily3 = squeeze(ncread([pth_out3 'evspsblsoi_daily.nc'],'evspsblsoi')); 
canopyevap_daily3 = squeeze(ncread([pth_out3 'evspsblveg_daily.nc'],'evspsblveg')); 
tran_daily3     = squeeze(ncread([pth_out3 'tran_daily.nc'],'tran')); %transpiration (kg m-2 s-1)
pr_daily3     = squeeze(ncread([pth_out3 'pr_daily.nc'],'pr')); %precipitation (kg m-2 s-1)
snw_daily3     = squeeze(ncread([pth_out3 'snw_daily.nc'],'snw')); %surface snow amount (kg m-2)
wsnw_daily3     = squeeze(ncread([pth_out3 'wsnw_daily.nc'],'wsnw')); %Liquid water content of snowpack (kg m-2)
snwdens_daily3     = squeeze(ncread([pth_out3 'snwdens_daily.nc'],'snwdens')); %density of snow (kg m-3)
tsn_daily3  = squeeze(ncread([pth_out3 'tsn_daily.nc'],'tsn')); %snowpack T
zsnw_daily3 = snw_daily3./snwdens_daily3; % snow depth (m)
gleafmas_daily3 = squeeze(ncread([pth_out3 'gleafmas_daily.nc'],'gleafmas')); %Green leaf mass (kg C m-2)
bleafmas_daily3 = squeeze(ncread([pth_out3 'bleafmas_daily.nc'],'bleafmas')); %Brown leaf mass (kg C m-2)
stemmass_daily3 = squeeze(ncread([pth_out3 'stemmass_daily.nc'],'stemmass')); %Stem mass (kg C m-2)
rootmass_daily3 = squeeze(ncread([pth_out3 'cRoot_daily.nc'],'cRoot')); %Root mass (kg C m-2)

vegHeight_annually3 = squeeze(ncread([pth_out3 'vegHeight_annually.nc'],'vegHeight')); 
vegHeight_annually_perpft3 = squeeze(ncread([pth_out3 'vegHeight_annually_perpft.nc'],['vegHeight'])); 

actlyrmax_monthly3 = squeeze(ncread([pth_out3 'actlyrmax_monthly.nc'],'actlyrmax')); %
lai_monthly3       = squeeze(ncread([pth_out3 'lai_monthly.nc'],'lai')); %
lai_monthly_perpft3= squeeze(ncread([pth_out3 'lai_monthly_perpft.nc'],'lai')); 
cSoil_monthly3     = squeeze(ncread([pth_out3 'cSoil_monthly.nc'],'cSoil')); %soil C
cSoil_monthly_perpft3 = squeeze(ncread([pth_out3 'cSoil_monthly_perpft.nc'],'cSoil')); %soil C
t_monthly3         = squeeze(ncread([pth_out3 'cSoil_monthly.nc'],'time')); %mth
snw_monthly3       = squeeze(ncread([pth_out3 'snw_monthly.nc'],'snw')); %snow amount
pr_monthly3        = squeeze(ncread([pth_out3 'pr_monthly.nc'],'pr')); %precipitation (kg m-2 s-1)
groundevap_monthly3 = squeeze(ncread([pth_out3 'evspsblsoi_monthly.nc'],'evspsblsoi')); 
canopyevap_monthly3 = squeeze(ncread([pth_out3 'evspsblveg_monthly.nc'],'evspsblveg')); 
snw_monthly3       = squeeze(ncread([pth_out3 'snw_monthly.nc'],'snw')); %Mass of snow pack (kg m^{-2} )
tsl_monthly3       = squeeze(ncread([pth_out3 'tsl_monthly.nc'],'tsl'));
groundevap_monthly3 = squeeze(ncread([pth_out3 'evspsblsoi_monthly.nc'],'evspsblsoi')); 
canopyevap_monthly3 = squeeze(ncread([pth_out3 'evspsblveg_monthly.nc'],'evspsblveg')); 
ftable_monthly3 = squeeze(ncread([pth_out3 'ftable_monthly.nc'],'ftable')); % mean depth to frozen water table (m)
cLand_monthly3       = squeeze(ncread([pth_out3 'cLand_monthly.nc'],'cLand')); % C in all terrestrial C pools (kg C m-2)
cLand_monthly_perpft3= squeeze(ncread([pth_out3 'cLand_monthly_perpft.nc'],'cLand')); 
cLeaf_monthly3     = squeeze(ncread([pth_out3 'cLeaf_monthly.nc'],'cLeaf')); 
cLeaf_monthly_perpft3= squeeze(ncread([pth_out3 'cLeaf_monthly_perpft.nc'],'cLeaf')); 
cLitter_monthly3     = squeeze(ncread([pth_out3 'cLitter_monthly.nc'],'cLitter')); 
cLitter_monthly_perpft3= squeeze(ncread([pth_out3 'cLitter_monthly_perpft.nc'],'cLitter')); 
cRoot_monthly3     = squeeze(ncread([pth_out3 'cRoot_monthly.nc'],'cRoot')); % C in roots
cRoot_monthly_perpft3= squeeze(ncread([pth_out3 'cRoot_monthly_perpft.nc'],'cRoot'));
cStem_monthly3     = squeeze(ncread([pth_out3 'cStem_monthly.nc'],'cStem'));
cStem_monthly_perpft3 = squeeze(ncread([pth_out3 'cStem_monthly_perpft.nc'],'cStem'));
cVeg_monthly3     = squeeze(ncread([pth_out3 'cVeg_monthly.nc'],'cVeg'));
cVeg_monthly_perpft3  = squeeze(ncread([pth_out3 'cVeg_monthly_perpft.nc'],'cVeg'));
hfg_monthly3     = squeeze(ncread([pth_out3 'hfg_monthly.nc'],'hfg')); %ground heat flux (W m-2)
tsl_daily3  = tsl_daily3 - 273.15; %conversion to degC
rootdpth_monthly_perpft3 = squeeze(ncread([pth_out3 'rootdpth_monthly_perpft.nc'],'rootdpth')); %Root depth

rn_daily = rss_daily + rls_daily;   % daily net radiation (W m-2)
rn_daily2 = rss_daily2 + rls_daily2;
rn_daily3 = rss_daily3 + rls_daily3;

albsn_daily  = squeeze(ncread([pth_out 'albsn_daily.nc'],'albsn')); % snow albedo
albsvis_daily  = squeeze(ncread([pth_out 'albsvis_daily.nc'],'albsvis')); % visible albedo
albsir_daily  = squeeze(ncread([pth_out 'albsir_daily.nc'],'albsir')); % infrared albedo

albsn_daily2  = squeeze(ncread([pth_out2 'albsn_daily.nc'],'albsn')); % snow albedo
albsvis_daily2  = squeeze(ncread([pth_out2 'albsvis_daily.nc'],'albsvis')); % visible albedo
albsir_daily2  = squeeze(ncread([pth_out2 'albsir_daily.nc'],'albsir')); % infrared albedo

albsn_daily3  = squeeze(ncread([pth_out3 'albsn_daily.nc'],'albsn')); % snow albedo
albsvis_daily3  = squeeze(ncread([pth_out3 'albsvis_daily.nc'],'albsvis')); % visible albedo
albsir_daily3  = squeeze(ncread([pth_out3 'albsir_daily.nc'],'albsir')); % infrared albedo

albsn_daily4  = squeeze(ncread([pth_out4 'albsn_daily.nc'],'albsn')); % snow albedo
albsvis_daily4  = squeeze(ncread([pth_out4 'albsvis_daily.nc'],'albsvis')); % visible albedo
albsir_daily4  = squeeze(ncread([pth_out4 'albsir_daily.nc'],'albsir')); % infrared albedo

tas_daily4  = squeeze(ncread([pth_out4 'tas_daily.nc'],'tas')); %air T
tsl_daily4  = squeeze(ncread([pth_out4 'tsl_daily.nc'],'tsl')); %{\itT}_s
mrsl_daily4 = squeeze(ncread([pth_out4 'mrsll_daily.nc'],'mrsll')); %vol water content (liquid)
mrfl_daily4 = squeeze(ncread([pth_out4 'mrsfl_daily.nc'],'mrsfl')); %vol water content (frozen)
rls_daily4  = squeeze(ncread([pth_out4 'rls_daily.nc'],'rls')); %net longwave
rss_daily4  = squeeze(ncread([pth_out4 'rss_daily.nc'],'rss')); %net shortwave
hfss_daily4  = squeeze(ncread([pth_out4 'hfss_daily.nc'],'hfss')); %H
hfls_daily4  = squeeze(ncread([pth_out4 'hfls_daily.nc'],'hfls')); %LE
albs_daily4  = squeeze(ncread([pth_out4 'albs_daily.nc'],'albs')); %albedo
nep_daily4 = squeeze(ncread([pth_out4 'nep_daily.nc'],'nep')); %NEP
npp_daily4  = squeeze(ncread([pth_out4 'npp_daily.nc'],'npp')); %NPP
gpp_daily4  = squeeze(ncread([pth_out4 'gpp_daily.nc'],'gpp')); %GPP
ra_daily4  = squeeze(ncread([pth_out4 'ra_daily.nc'],'ra')); %Ra
rh_daily4  = squeeze(ncread([pth_out4 'rh_daily.nc'],'rh')); %Rh
nee_daily4 = -nep_daily4; %NEE  (kg C m-2 s-1)

mrro_daily4 = squeeze(ncread([pth_out4 'mrro_daily.nc'],'mrro')); %total runoff from soil
mrros_daily4 = squeeze(ncread([pth_out4 'mrros_daily.nc'],'mrros')); %overland flow from top of soil column
mrrob_daily4 = squeeze(ncread([pth_out4 'mrrob_daily.nc'],'mrrob')); %base flow from bottom
mrroc_daily4 = squeeze(ncread([pth_out4 'mrroc_daily.nc'],'mrroc')); %Liquid/frozen water runoff from vegetation
mrroi_daily4 = squeeze(ncread([pth_out4 'mrroi_daily.nc'],'mrroi')); %interflow
mrron_daily4 = squeeze(ncread([pth_out4 'mrron_daily.nc'],'mrron')); %Liquid runoff from snowpack

zpond_daily4 = squeeze(ncread([pth_out4 'zpond_daily.nc'],'zpond')); %Depth of ponded water [m]
evap_daily4 = squeeze(ncread([pth_out4 'evspsbl_daily.nc'],'evspsbl')); 
groundevap_daily4 = squeeze(ncread([pth_out4 'evspsblsoi_daily.nc'],'evspsblsoi')); 
canopyevap_daily4 = squeeze(ncread([pth_out4 'evspsblveg_daily.nc'],'evspsblveg')); 
tran_daily4     = squeeze(ncread([pth_out4 'tran_daily.nc'],'tran')); %transpiration (kg m-2 s-1)
pr_daily4     = squeeze(ncread([pth_out4 'pr_daily.nc'],'pr')); %precipitation (kg m-2 s-1)
snw_daily4     = squeeze(ncread([pth_out4 'snw_daily.nc'],'snw')); %surface snow amount (kg m-2)
wsnw_daily4     = squeeze(ncread([pth_out4 'wsnw_daily.nc'],'wsnw')); %Liquid water content of snowpack (kg m-2)
snwdens_daily4     = squeeze(ncread([pth_out4 'snwdens_daily.nc'],'snwdens')); %density of snow (kg m-3)
tsn_daily4  = squeeze(ncread([pth_out4 'tsn_daily.nc'],'tsn')); %snowpack T
zsnw_daily4 = snw_daily4./snwdens_daily4; % snow depth (m)
gleafmas_daily4 = squeeze(ncread([pth_out4 'gleafmas_daily.nc'],'gleafmas')); %Green leaf mass (kg C m-2)
bleafmas_daily4 = squeeze(ncread([pth_out4 'bleafmas_daily.nc'],'bleafmas')); %Brown leaf mass (kg C m-2)
stemmass_daily4 = squeeze(ncread([pth_out4 'stemmass_daily.nc'],'stemmass')); %Stem mass (kg C m-2)
rootmass_daily4 = squeeze(ncread([pth_out4 'cRoot_daily.nc'],'cRoot')); %Root mass (kg C m-2)

vegHeight_annually4 = squeeze(ncread([pth_out4 'vegHeight_annually.nc'],'vegHeight')); 
vegHeight_annually_perpft4 = squeeze(ncread([pth_out4 'vegHeight_annually_perpft.nc'],['vegHeight'])); 

actlyrmax_monthly4 = squeeze(ncread([pth_out4 'actlyrmax_monthly.nc'],'actlyrmax')); %
lai_monthly4       = squeeze(ncread([pth_out4 'lai_monthly.nc'],'lai')); %
lai_monthly_perpft4= squeeze(ncread([pth_out4 'lai_monthly_perpft.nc'],'lai')); 
cSoil_monthly4     = squeeze(ncread([pth_out4 'cSoil_monthly.nc'],'cSoil')); %soil C
cSoil_monthly_perpft4 = squeeze(ncread([pth_out4 'cSoil_monthly_perpft.nc'],'cSoil')); %soil C
t_monthly4         = squeeze(ncread([pth_out4 'cSoil_monthly.nc'],'time')); %mth
snw_monthly4       = squeeze(ncread([pth_out4 'snw_monthly.nc'],'snw')); %snow amount
pr_monthly4        = squeeze(ncread([pth_out4 'pr_monthly.nc'],'pr')); %precipitation (kg m-2 s-1)
groundevap_monthly4 = squeeze(ncread([pth_out4 'evspsblsoi_monthly.nc'],'evspsblsoi')); 
canopyevap_monthly4 = squeeze(ncread([pth_out4 'evspsblveg_monthly.nc'],'evspsblveg')); 
snw_monthly4       = squeeze(ncread([pth_out4 'snw_monthly.nc'],'snw')); %Mass of snow pack (kg m^{-2} )
tsl_monthly4       = squeeze(ncread([pth_out4 'tsl_monthly.nc'],'tsl'));
groundevap_monthly4 = squeeze(ncread([pth_out4 'evspsblsoi_monthly.nc'],'evspsblsoi')); 
canopyevap_monthly4 = squeeze(ncread([pth_out4 'evspsblveg_monthly.nc'],'evspsblveg')); 
ftable_monthly4 = squeeze(ncread([pth_out4 'ftable_monthly.nc'],'ftable')); % mean depth to frozen water table (m)
cLand_monthly4       = squeeze(ncread([pth_out4 'cLand_monthly.nc'],'cLand')); % C in all terrestrial C pools (kg C m-2)
cLand_monthly_perpft4= squeeze(ncread([pth_out4 'cLand_monthly_perpft.nc'],'cLand')); 
cLeaf_monthly4     = squeeze(ncread([pth_out4 'cLeaf_monthly.nc'],'cLeaf')); 
cLeaf_monthly_perpft4= squeeze(ncread([pth_out4 'cLeaf_monthly_perpft.nc'],'cLeaf')); 
cLitter_monthly4     = squeeze(ncread([pth_out4 'cLitter_monthly.nc'],'cLitter')); 
cLitter_monthly_perpft4= squeeze(ncread([pth_out4 'cLitter_monthly_perpft.nc'],'cLitter')); 
cRoot_monthly4     = squeeze(ncread([pth_out4 'cRoot_monthly.nc'],'cRoot')); % C in roots
cRoot_monthly_perpft4= squeeze(ncread([pth_out4 'cRoot_monthly_perpft.nc'],'cRoot'));
cStem_monthly4     = squeeze(ncread([pth_out4 'cStem_monthly.nc'],'cStem'));
cStem_monthly_perpft4 = squeeze(ncread([pth_out4 'cStem_monthly_perpft.nc'],'cStem'));
cVeg_monthly4     = squeeze(ncread([pth_out4 'cVeg_monthly.nc'],'cVeg'));
cVeg_monthly_perpft4  = squeeze(ncread([pth_out4 'cVeg_monthly_perpft.nc'],'cVeg'));
hfg_monthly4     = squeeze(ncread([pth_out4 'hfg_monthly.nc'],'hfg')); %ground heat flux (W m-2)
rootdpth_monthly_perpft4 = squeeze(ncread([pth_out4 'rootdpth_monthly_perpft.nc'],'rootdpth')); %Root depth

rn_daily4 = rss_daily4 + rls_daily4;
tsl_daily4  = tsl_daily4 - 273.15; %conversion to degC
%%
[yy_av, MM_av, dd_av, ~, ~] = datevec(tv_av);

%% calculate mean and standard deviation of climate variables and fluxes grouped by DOY
DL_mod_T = table(t1,yy,MM,dd,doy,nep_daily.*1000.*60*60*24,gpp_daily.*1000.*60*60*24,...
    ra_daily.*1000.*60*60*24,rh_daily.*1000.*60*60*24,hfls_daily,hfss_daily,...
    albs_daily,tas_daily,zsnw_daily,snwdens_daily,evap_daily,tran_daily,groundevap_daily,...
    canopyevap_daily,pr_daily.*86400,zpond_daily,rn_daily,...
    nee_daily.*1000.*60*60*24,mrro_daily,albsvis_daily,albsir_daily);
DL_mod_T.Properties.VariableNames = {'t1','yy','MM','dd','doy','nep','gpp','ra','rh','hfls',...
    'hfss','albs','tas','zsnw','snwdens','evap','tran','groundevap','canopyevap','precip',...
    'zpond','rn','nee','runoff','albsvis','albsir'};
DL_mod_T.tsl1 = tsl_daily(1,:)';
DL_mod_T.tsl3 = tsl_daily(3,:)';
DL_mod_T.tsl6 = tsl_daily(6,:)';
DL_mod_T.mrsl1 = mrsl_daily(1,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod_T.mrsl2 = mrsl_daily(2,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod_T.mrsl3 = mrsl_daily(3,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod_T.mrsl_rel1 = ((mrsl_daily(1,:)-min(mrsl_daily(1,:)))./(max(mrsl_daily(1,:))-min(mrsl_daily(1,:))))';
DL_mod_T.mrsl_rel2 = ((mrsl_daily(2,:)-min(mrsl_daily(2,:)))./(max(mrsl_daily(2,:))-min(mrsl_daily(2,:))))';
DL_mod_T.mrsl_rel3 = ((mrsl_daily(3,:)-min(mrsl_daily(3,:)))./(max(mrsl_daily(3,:))-min(mrsl_daily(3,:))))';
DL_mod_all_T = DL_mod_T;
stats_DL_mod_T = grpstats(DL_mod_T,{'doy'},{'mean','std'},'DataVars',{'nep','gpp',...
    'ra','rh','hfls','hfss','albs','tas','zsnw','snwdens','tsl1','tsl3','tsl6',...
    'mrsl1','mrsl2','mrsl3','mrsl_rel1','mrsl_rel2','mrsl_rel3','evap','tran',...
    'groundevap','canopyevap','precip','zpond','rn','nee','runoff',...
    'albsvis','albsir'});

DL_mod2_T = table(t2,yy2,MM2,dd2,doy2,nep_daily2.*1000.*60*60*24,gpp_daily2.*1000.*60*60*24,...
    ra_daily2.*1000.*60*60*24,rh_daily2.*1000.*60*60*24,hfls_daily2,hfss_daily2,...
    albs_daily2,tas_daily2,zsnw_daily2,snwdens_daily2,evap_daily2,tran_daily2,...
    groundevap_daily2,canopyevap_daily2,pr_daily2.*86400,zpond_daily2,...
    rn_daily2, nee_daily2.*1000.*60*60*24,mrro_daily2,albsvis_daily2,albsir_daily2);
DL_mod2_T.Properties.VariableNames = {'t1','yy','MM','dd','doy','nep','gpp','ra','rh','hfls',...
    'hfss','albs','tas','zsnw','snwdens','evap','tran','groundevap','canopyevap','precip',...
    'zpond','rn','nee','runoff','albsvis','albsir'};
DL_mod2_T.tsl1 = tsl_daily2(1,:)';
DL_mod2_T.tsl3 = tsl_daily2(3,:)';
DL_mod2_T.tsl6 = tsl_daily2(6,:)';
DL_mod2_T.mrsl1 = mrsl_daily2(1,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod2_T.mrsl2 = mrsl_daily2(2,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod2_T.mrsl3 = mrsl_daily2(3,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod2_T.mrsl_rel1 = ((mrsl_daily2(1,:)-min(mrsl_daily2(1,:)))./(max(mrsl_daily2(1,:))-min(mrsl_daily2(1,:))))';
DL_mod2_T.mrsl_rel2 = ((mrsl_daily2(2,:)-min(mrsl_daily2(2,:)))./(max(mrsl_daily2(2,:))-min(mrsl_daily2(2,:))))';
DL_mod2_T.mrsl_rel3 = ((mrsl_daily2(3,:)-min(mrsl_daily2(3,:)))./(max(mrsl_daily2(3,:))-min(mrsl_daily2(3,:))))';
DL_mod2_all_T = DL_mod2_T;
ind_mod2_2018 = find(DL_mod2_all_T.yy >= 2018);
DL_mod2_T(ind_mod2_2018,:) = [];
stats_DL_mod2_T = grpstats(DL_mod2_T,{'doy'},{'mean','std'},'DataVars',{'nep','gpp',...
    'ra','rh','hfls','hfss','albs','tas','zsnw','snwdens','tsl1','tsl3','tsl6',...
    'mrsl1','mrsl2','mrsl3','mrsl_rel1','mrsl_rel2','mrsl_rel3','evap','tran',...
    'groundevap','canopyevap','precip','zpond','rn','nee','runoff',...
    'albsvis','albsir'});

DL_mod3_T = table(t2,yy2,MM2,dd2,doy2,nep_daily3.*1000.*60*60*24,gpp_daily3.*1000.*60*60*24,...
    ra_daily3.*1000.*60*60*24,rh_daily3.*1000.*60*60*24,hfls_daily3,hfss_daily3,...
    albs_daily3,tas_daily3,zsnw_daily3,snwdens_daily3,evap_daily3,tran_daily3,...
    groundevap_daily3,canopyevap_daily3,pr_daily3.*86400,zpond_daily3,...
    rn_daily3,nee_daily3.*1000.*60*60*24,mrro_daily3,albsvis_daily3,albsir_daily3);
DL_mod3_T.Properties.VariableNames = {'t1','yy','MM','dd','doy','nep','gpp','ra','rh','hfls',...
    'hfss','albs','tas','zsnw','snwdens','evap','tran','groundevap','canopyevap','precip',...
    'zpond','rn','nee','runoff','albsvis','albsir'};
DL_mod3_T.tsl1 = tsl_daily3(1,:)';
DL_mod3_T.tsl3 = tsl_daily3(3,:)';
DL_mod3_T.tsl6 = tsl_daily3(6,:)';
DL_mod3_T.mrsl1 = mrsl_daily3(1,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod3_T.mrsl2 = mrsl_daily3(2,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod3_T.mrsl3 = mrsl_daily3(3,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod3_T.mrsl_rel1 = ((mrsl_daily3(1,:)-min(mrsl_daily3(1,:)))./(max(mrsl_daily3(1,:))-min(mrsl_daily3(1,:))))';
DL_mod3_T.mrsl_rel2 = ((mrsl_daily3(2,:)-min(mrsl_daily3(2,:)))./(max(mrsl_daily3(2,:))-min(mrsl_daily3(2,:))))';
DL_mod3_T.mrsl_rel3 = ((mrsl_daily3(3,:)-min(mrsl_daily3(3,:)))./(max(mrsl_daily3(3,:))-min(mrsl_daily3(3,:))))';
DL_mod3_all_T = DL_mod3_T;
DL_mod3_T(ind_mod2_2018,:) = [];
stats_DL_mod3_T = grpstats(DL_mod3_T,{'doy'},{'mean','std'},'DataVars',{'nep','gpp',...
    'ra','rh','hfls','hfss','albs','tas','zsnw','snwdens','tsl1','tsl3','tsl6',...
    'mrsl1','mrsl2','mrsl3','mrsl_rel1','mrsl_rel2','mrsl_rel3','evap','tran',...
    'groundevap','canopyevap','precip','zpond','rn','nee','runoff',...
    'albsvis','albsir'});

DL_meas_T = table(tv_av,yy_av,MM_av,dd_av,doy,NEP,GEP,Re,...
    LE,H,snow,ET,Rn);
DL_meas_T.Properties.VariableNames = {'tv_av','yy','MM','dd','doy','nep','gep','re','le',...
    'hs','snow','et','rn'};
DL_meas_T.Tsoil1 = Tsoil1;
DL_meas_T.Tsoil3 = Tsoil3;
DL_meas_T.Tsoil6 = Tsoil6;
DL_meas_T.sm1 = VWC1;
DL_meas_T.sm3 = VWC3;
DL_meas_T.sm_rel1 = ((VWC1-min(VWC1))./(max(VWC1)-min(VWC1)));
DL_meas_T.sm_rel3 = ((VWC3-min(VWC3))./(max(VWC3)-min(VWC3)));
stats_DL_meas_T = grpstats(DL_meas_T,{'doy'},{'mean','std'},'DataVars',{'nep','gep',...
    're','le','hs','snow','Tsoil1','Tsoil3','Tsoil6','sm1','sm3',...
    'sm_rel1','sm_rel3','et','rn'});

DL_mod4_T = table(t2,yy2,MM2,dd2,doy2,nep_daily4.*1000.*60*60*24,gpp_daily4.*1000.*60*60*24,...
    ra_daily4.*1000.*60*60*24,rh_daily4.*1000.*60*60*24,hfls_daily4,hfss_daily4,...
    albs_daily4,tas_daily4,zsnw_daily4,snwdens_daily4,evap_daily4,tran_daily4,...
    groundevap_daily4,canopyevap_daily4,pr_daily4.*86400,zpond_daily4,...
    rn_daily4, nee_daily4.*1000.*60*60*24,mrro_daily4,albsvis_daily4,albsir_daily4);
DL_mod4_T.Properties.VariableNames = {'t1','yy','MM','dd','doy','nep','gpp','ra','rh','hfls',...
    'hfss','albs','tas','zsnw','snwdens','evap','tran','groundevap','canopyevap','precip',...
    'zpond','rn','nee','runoff','albsvis','albsir'};
DL_mod4_T.tsl1 = tsl_daily4(1,:)';
DL_mod4_T.tsl3 = tsl_daily4(3,:)';
DL_mod4_T.tsl6 = tsl_daily4(6,:)';
DL_mod4_T.mrsl1 = mrsl_daily4(1,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod4_T.mrsl2 = mrsl_daily4(2,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod4_T.mrsl3 = mrsl_daily4(3,:)'./100; %conversion from kg/m2 to m3/m3
DL_mod4_T.mrsl_rel1 = ((mrsl_daily4(1,:)-min(mrsl_daily4(1,:)))./(max(mrsl_daily4(1,:))-min(mrsl_daily4(1,:))))';
DL_mod4_T.mrsl_rel2 = ((mrsl_daily4(2,:)-min(mrsl_daily4(2,:)))./(max(mrsl_daily4(2,:))-min(mrsl_daily4(2,:))))';
DL_mod4_T.mrsl_rel3 = ((mrsl_daily4(3,:)-min(mrsl_daily4(3,:)))./(max(mrsl_daily4(3,:))-min(mrsl_daily4(3,:))))';
DL_mod4_all_T = DL_mod4_T;
ind_mod4_2018 = find(DL_mod4_all_T.yy >= 2018);
DL_mod4_T(ind_mod4_2018,:) = [];
stats_DL_mod4_T = grpstats(DL_mod4_T,{'doy'},{'mean','std'},'DataVars',{'nep','gpp',...
    'ra','rh','hfls','hfss','albs','tas','zsnw','snwdens','tsl1','tsl3','tsl6',...
    'mrsl1','mrsl2','mrsl3','mrsl_rel1','mrsl_rel2','mrsl_rel3','evap','tran',...
    'groundevap','canopyevap','precip','zpond','rn','nee','runoff',...
    'albsvis','albsir'});
%%
% 5-day averages
% mod
[doy_5av, nep_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_nep, [], [], 5, 1);
[~, nep_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_nep, [], [], 5, 1);
[~, gpp_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_gpp, [], [], 5, 1);
[~, gpp_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_gpp, [], [], 5, 1);
[~, ra_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_ra, [], [], 5, 1);
[~, ra_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_ra, [], [], 5, 1);
[~, rh_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_rh, [], [], 5, 1);
[~, rh_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_rh, [], [], 5, 1);
[~, hfls_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_hfls, [], [], 5, 1);
[~, hfls_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_hfls, [], [], 5, 1);
[~, hfss_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_hfss, [], [], 5, 1);
[~, hfss_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_hfss, [], [], 5, 1);
[~, albs_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_albs, [], [], 5, 1);
[~, albs_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_albs, [], [], 5, 1);
[~, albsvis_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_albsvis, [], [], 5, 1);
[~, albsir_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_albsir, [], [], 5, 1);
[~, albsvis_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_albsvis, [], [], 5, 1);
[~, albsir_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_albsir, [], [], 5, 1);
[~, tas_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_tas, [], [], 5, 1);
[~, tas_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_tas, [], [], 5, 1);
[~, zsnw_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_zsnw, [], [], 5, 1);
[~, zsnw_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_zsnw, [], [], 5, 1);
[~, snwdens_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_snwdens, [], [], 5, 1);
[~, snwdens_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_snwdens, [], [], 5, 1);
[~, tsl1_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_tsl1, [], [], 5, 1);
[~, tsl1_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_tsl1, [], [], 5, 1);
[~, tsl3_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_tsl3, [], [], 5, 1);
[~, tsl3_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_tsl3, [], [], 5, 1);
[~, tsl6_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_tsl6, [], [], 5, 1);
[~, tsl6_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_tsl6, [], [], 5, 1);
[~, mrsl1_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_mrsl1, [], [], 5, 1);
[~, mrsl1_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_mrsl1, [], [], 5, 1);
[~, mrsl2_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_mrsl2, [], [], 5, 1);
[~, mrsl2_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_mrsl2, [], [], 5, 1);
[~, mrsl3_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_mrsl3, [], [], 5, 1);
[~, mrsl3_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_mrsl3, [], [], 5, 1);
[~, mrsl_rel1_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_mrsl_rel1, [], [], 5, 1);
[~, mrsl_rel1_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_mrsl_rel1, [], [], 5, 1);
[~, mrsl_rel2_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_mrsl_rel2, [], [], 5, 1);
[~, mrsl_rel2_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_mrsl_rel2, [], [], 5, 1);
[~, mrsl_rel3_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_mrsl_rel3, [], [], 5, 1);
[~, mrsl_rel3_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_mrsl_rel3, [], [], 5, 1);
[~, evap_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_evap, [], [], 5, 1);
[~, evap_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_evap, [], [], 5, 1);
[~, tran_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_tran, [], [], 5, 1);
[~, tran_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_tran, [], [], 5, 1);
[~, groundevap_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_groundevap, [], [], 5, 1);
[~, groundevap_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_groundevap, [], [], 5, 1);
[~, canopyevap_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_canopyevap, [], [], 5, 1);
[~, canopyevap_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_canopyevap, [], [], 5, 1);
[~, precip_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_precip, [], [], 5, 1);
[~, precip_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_precip, [], [], 5, 1);
[~, zpond_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_zpond, [], [], 5, 1);
[~, zpond_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_zpond, [], [], 5, 1);
[~, rn_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_rn, [], [], 5, 1);
[~, rn_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_rn, [], [], 5, 1);
[~, nee_mod_5av] = dailystats_number(doy,stats_DL_mod_T.mean_nee, [], [], 5, 1);
[~, nee_std_mod_5av] = dailystats_number(doy,stats_DL_mod_T.std_nee, [], [], 5, 1);
% mod2
[~, nep_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_nep, [], [], 5, 1);
[~, nep_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_nep, [], [], 5, 1);
[~, gpp_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_gpp, [], [], 5, 1);
[~, gpp_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_gpp, [], [], 5, 1);
[~, ra_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_ra, [], [], 5, 1);
[~, ra_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_ra, [], [], 5, 1);
[~, rh_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_rh, [], [], 5, 1);
[~, rh_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_rh, [], [], 5, 1);
[~, hfls_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_hfls, [], [], 5, 1);
[~, hfls_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_hfls, [], [], 5, 1);
[~, hfss_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_hfss, [], [], 5, 1);
[~, hfss_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_hfss, [], [], 5, 1);
[~, albs_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_albs, [], [], 5, 1);
[~, albsvis_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_albsvis, [], [], 5, 1);
[~, albsir_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_albsir, [], [], 5, 1);
[~, albsvis_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_albsvis, [], [], 5, 1);
[~, albsir_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_albsir, [], [], 5, 1);
[~, albs_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_albs, [], [], 5, 1);
[~, tas_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_tas, [], [], 5, 1);
[~, tas_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_tas, [], [], 5, 1);
[~, zsnw_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_zsnw, [], [], 5, 1);
[~, zsnw_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_zsnw, [], [], 5, 1);
[~, snwdens_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_snwdens, [], [], 5, 1);
[~, snwdens_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_snwdens, [], [], 5, 1);
[~, tsl1_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_tsl1, [], [], 5, 1);
[~, tsl1_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_tsl1, [], [], 5, 1);
[~, tsl3_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_tsl3, [], [], 5, 1);
[~, tsl3_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_tsl3, [], [], 5, 1);
[~, tsl6_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_tsl6, [], [], 5, 1);
[~, tsl6_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_tsl6, [], [], 5, 1);
[~, mrsl1_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_mrsl1, [], [], 5, 1);
[~, mrsl1_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_mrsl1, [], [], 5, 1);
[~, mrsl2_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_mrsl2, [], [], 5, 1);
[~, mrsl2_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_mrsl2, [], [], 5, 1);
[~, mrsl3_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_mrsl3, [], [], 5, 1);
[~, mrsl3_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_mrsl3, [], [], 5, 1);
[~, mrsl_rel1_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_mrsl_rel1, [], [], 5, 1);
[~, mrsl_rel1_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_mrsl_rel1, [], [], 5, 1);
[~, mrsl_rel2_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_mrsl_rel2, [], [], 5, 1);
[~, mrsl_rel2_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_mrsl_rel2, [], [], 5, 1);
[~, mrsl_rel3_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_mrsl_rel3, [], [], 5, 1);
[~, mrsl_rel3_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_mrsl_rel3, [], [], 5, 1);
[~, evap_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_evap, [], [], 5, 1);
[~, evap_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_evap, [], [], 5, 1);
[~, tran_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_tran, [], [], 5, 1);
[~, tran_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_tran, [], [], 5, 1);
[~, groundevap_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_groundevap, [], [], 5, 1);
[~, groundevap_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_groundevap, [], [], 5, 1);
[~, canopyevap_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_canopyevap, [], [], 5, 1);
[~, canopyevap_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_canopyevap, [], [], 5, 1);
[~, precip_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_precip, [], [], 5, 1);
[~, precip_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_precip, [], [], 5, 1);
[~, zpond_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_zpond, [], [], 5, 1);
[~, zpond_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_zpond, [], [], 5, 1);
[~, rn_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_rn, [], [], 5, 1);
[~, rn_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_rn, [], [], 5, 1);
[~, nee_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.mean_nee, [], [], 5, 1);
[~, nee_std_mod2_5av] = dailystats_number(doy,stats_DL_mod2_T.std_nee, [], [], 5, 1);
% mod3
[~, nep_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_nep, [], [], 5, 1);
[~, nep_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_nep, [], [], 5, 1);
[~, gpp_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_gpp, [], [], 5, 1);
[~, gpp_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_gpp, [], [], 5, 1);
[~, ra_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_ra, [], [], 5, 1);
[~, ra_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_ra, [], [], 5, 1);
[~, rh_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_rh, [], [], 5, 1);
[~, rh_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_rh, [], [], 5, 1);
[~, hfls_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_hfls, [], [], 5, 1);
[~, hfls_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_hfls, [], [], 5, 1);
[~, hfss_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_hfss, [], [], 5, 1);
[~, hfss_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_hfss, [], [], 5, 1);
[~, albs_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_albs, [], [], 5, 1);
[~, albs_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_albs, [], [], 5, 1);
[~, albsvis_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_albsvis, [], [], 5, 1);
[~, albsir_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_albsir, [], [], 5, 1);
[~, albsvis_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_albsvis, [], [], 5, 1);
[~, albsir_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_albsir, [], [], 5, 1);
[~, tas_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_tas, [], [], 5, 1);
[~, tas_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_tas, [], [], 5, 1);
[~, zsnw_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_zsnw, [], [], 5, 1);
[~, zsnw_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_zsnw, [], [], 5, 1);
[~, snwdens_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_snwdens, [], [], 5, 1);
[~, snwdens_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_snwdens, [], [], 5, 1);
[~, tsl1_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_tsl1, [], [], 5, 1);
[~, tsl1_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_tsl1, [], [], 5, 1);
[~, tsl3_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_tsl3, [], [], 5, 1);
[~, tsl3_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_tsl3, [], [], 5, 1);
[~, tsl6_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_tsl6, [], [], 5, 1);
[~, tsl6_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_tsl6, [], [], 5, 1);
[~, mrsl1_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_mrsl1, [], [], 5, 1);
[~, mrsl1_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_mrsl1, [], [], 5, 1);
[~, mrsl2_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_mrsl2, [], [], 5, 1);
[~, mrsl2_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_mrsl2, [], [], 5, 1);
[~, mrsl3_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_mrsl3, [], [], 5, 1);
[~, mrsl3_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_mrsl3, [], [], 5, 1);
[~, mrsl_rel1_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_mrsl_rel1, [], [], 5, 1);
[~, mrsl_rel1_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_mrsl_rel1, [], [], 5, 1);
[~, mrsl_rel2_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_mrsl_rel2, [], [], 5, 1);
[~, mrsl_rel2_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_mrsl_rel2, [], [], 5, 1);
[~, mrsl_rel3_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_mrsl_rel3, [], [], 5, 1);
[~, mrsl_rel3_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_mrsl_rel3, [], [], 5, 1);
[~, evap_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_evap, [], [], 5, 1);
[~, evap_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_evap, [], [], 5, 1);
[~, tran_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_tran, [], [], 5, 1);
[~, tran_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_tran, [], [], 5, 1);
[~, groundevap_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_groundevap, [], [], 5, 1);
[~, groundevap_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_groundevap, [], [], 5, 1);
[~, canopyevap_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_canopyevap, [], [], 5, 1);
[~, canopyevap_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_canopyevap, [], [], 5, 1);
[~, precip_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_precip, [], [], 5, 1);
[~, precip_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_precip, [], [], 5, 1);
[~, zpond_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_zpond, [], [], 5, 1);
[~, zpond_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_zpond, [], [], 5, 1);
[~, rn_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_rn, [], [], 5, 1);
[~, rn_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_rn, [], [], 5, 1);
[~, nee_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.mean_nee, [], [], 5, 1);
[~, nee_std_mod3_5av] = dailystats_number(doy,stats_DL_mod3_T.std_nee, [], [], 5, 1);
% mod4
[~, nep_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_nep, [], [], 5, 1);
[~, nep_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_nep, [], [], 5, 1);
[~, gpp_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_gpp, [], [], 5, 1);
[~, gpp_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_gpp, [], [], 5, 1);
[~, ra_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_ra, [], [], 5, 1);
[~, ra_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_ra, [], [], 5, 1);
[~, rh_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_rh, [], [], 5, 1);
[~, rh_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_rh, [], [], 5, 1);
[~, hfls_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_hfls, [], [], 5, 1);
[~, hfls_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_hfls, [], [], 5, 1);
[~, hfss_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_hfss, [], [], 5, 1);
[~, hfss_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_hfss, [], [], 5, 1);
[~, albs_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_albs, [], [], 5, 1);
[~, albsvis_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_albsvis, [], [], 5, 1);
[~, albsir_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_albsir, [], [], 5, 1);
[~, albsvis_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_albsvis, [], [], 5, 1);
[~, albsir_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_albsir, [], [], 5, 1);
[~, albs_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_albs, [], [], 5, 1);
[~, tas_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_tas, [], [], 5, 1);
[~, tas_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_tas, [], [], 5, 1);
[~, zsnw_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_zsnw, [], [], 5, 1);
[~, zsnw_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_zsnw, [], [], 5, 1);
[~, snwdens_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_snwdens, [], [], 5, 1);
[~, snwdens_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_snwdens, [], [], 5, 1);
[~, tsl1_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_tsl1, [], [], 5, 1);
[~, tsl1_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_tsl1, [], [], 5, 1);
[~, tsl3_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_tsl3, [], [], 5, 1);
[~, tsl3_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_tsl3, [], [], 5, 1);
[~, tsl6_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_tsl6, [], [], 5, 1);
[~, tsl6_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_tsl6, [], [], 5, 1);
[~, mrsl1_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_mrsl1, [], [], 5, 1);
[~, mrsl1_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_mrsl1, [], [], 5, 1);
[~, mrsl2_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_mrsl2, [], [], 5, 1);
[~, mrsl2_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_mrsl2, [], [], 5, 1);
[~, mrsl3_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_mrsl3, [], [], 5, 1);
[~, mrsl3_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_mrsl3, [], [], 5, 1);
[~, mrsl_rel1_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_mrsl_rel1, [], [], 5, 1);
[~, mrsl_rel1_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_mrsl_rel1, [], [], 5, 1);
[~, mrsl_rel2_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_mrsl_rel2, [], [], 5, 1);
[~, mrsl_rel2_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_mrsl_rel2, [], [], 5, 1);
[~, mrsl_rel3_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_mrsl_rel3, [], [], 5, 1);
[~, mrsl_rel3_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_mrsl_rel3, [], [], 5, 1);
[~, evap_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_evap, [], [], 5, 1);
[~, evap_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_evap, [], [], 5, 1);
[~, tran_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_tran, [], [], 5, 1);
[~, tran_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_tran, [], [], 5, 1);
[~, groundevap_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_groundevap, [], [], 5, 1);
[~, groundevap_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_groundevap, [], [], 5, 1);
[~, canopyevap_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_canopyevap, [], [], 5, 1);
[~, canopyevap_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_canopyevap, [], [], 5, 1);
[~, precip_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_precip, [], [], 5, 1);
[~, precip_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_precip, [], [], 5, 1);
[~, zpond_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_zpond, [], [], 5, 1);
[~, zpond_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_zpond, [], [], 5, 1);
[~, rn_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_rn, [], [], 5, 1);
[~, rn_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_rn, [], [], 5, 1);
[~, nee_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.mean_nee, [], [], 5, 1);
[~, nee_std_mod4_5av] = dailystats_number(doy,stats_DL_mod4_T.std_nee, [], [], 5, 1);
% meas
[~, nep_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_nep, [], [], 5, 1);
[~, nep_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_nep, [], [], 5, 1);
[~, gep_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_gep, [], [], 5, 1);
[~, gep_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_gep, [], [], 5, 1);
[~, re_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_re, [], [], 5, 1);
[~, re_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_re, [], [], 5, 1);
%[~, le_meas_5av] = dailystats_number(doy,stats_DL_meas_T.sum_le, [], [], 5, 1);
[~, le_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_le, [], [], 5, 1);
[~, le_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_le, [], [], 5, 1);
%[~, hs_meas_5av] = dailystats_number(doy,stats_DL_meas_T.sum_hs, [], [], 5, 1);
[~, hs_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_hs, [], [], 5, 1);
[~, hs_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_hs, [], [], 5, 1);
[~, snow_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_snow, [], [], 5, 1);
[~, snow_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_snow, [], [], 5, 1);
[~, Tsoil1_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_Tsoil1, [], [], 5, 1);
[~, Tsoil1_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_Tsoil1, [], [], 5, 1);
[~, Tsoil3_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_Tsoil3, [], [], 5, 1);
[~, Tsoil3_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_Tsoil3, [], [], 5, 1);
[~, Tsoil6_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_Tsoil6, [], [], 5, 1);
[~, Tsoil6_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_Tsoil6, [], [], 5, 1);
[~, sm1_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_sm1, [], [], 5, 1);
[~, sm1_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_sm1, [], [], 5, 1);
[~, sm3_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_sm3, [], [], 5, 1);
[~, sm3_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_sm3, [], [], 5, 1);
[~, sm_rel1_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_sm_rel1, [], [], 5, 1);
[~, sm_rel1_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_sm_rel1, [], [], 5, 1);
[~, sm_rel3_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_sm_rel3, [], [], 5, 1);
[~, sm_rel3_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_sm_rel3, [], [], 5, 1);
[~, et_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_et, [], [], 5, 1);
[~, et_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_et, [], [], 5, 1);
[~, rn_meas_5av] = dailystats_number(doy,stats_DL_meas_T.mean_rn, [], [], 5, 1);
[~, rn_std_meas_5av] = dailystats_number(doy,stats_DL_meas_T.std_rn, [], [], 5, 1);
[~, chambers_mean_year_5av] = dailystats_number(doy(1:366),Re_chambers(1:366), [], [], 5, 1);

[~, chambers_std_year_5av] = dailystats_number(doy(1:366),Re_chambers(1:366), [], [], 5, 1);
nep_chambers_mean_year_5av = chambers_mean_year_5av*(-1);
%%
fig = 0;
%%
%plots for average year  
%%
      fig = eh_set_fig(fig);
  [hl, hp] = boundedline(doy_5av', zsnw_mod_5av,zsnw_std_mod_5av,'-b',...
      doy_5av',snow_meas_5av,snow_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[0 0.5]);
  xlabel('DOY');
  ylabel('Snow depth (m)'); legend('Shrubs','Obs')
%%
fig = eh_set_fig(fig);
set(gcf, 'Units', 'Centimeters', 'Position', [0 0 24 18]);
subplot(2,2,1)
 [hl, hp] = boundedline(doy_5av', hfls_mod_5av,hfls_std_mod_5av,'-b',...
      doy_5av',hfls_mod2_5av,hfls_std_mod2_5av,'-y',...
      doy_5av',hfls_mod3_5av,hfls_std_mod3_5av,'-r',...
      doy_5av',le_meas_5av,le_std_meas_5av,'-.k','alpha','nan','gap');
  set(gca,'XLim',[0 366],'YLim',[-10 135]);
  ylabel('LE (W m^{-2})'); legend('Shrubs','Grass','Trees','Obs')
  text(7,130,'a)','FontSize',12);

 subplot(2,2,2)
 [hl, hp] = boundedline(doy_5av', hfss_mod_5av,hfss_std_mod_5av,'-b',...
      doy_5av',hfss_mod2_5av,hfss_std_mod2_5av,'-y',...
      doy_5av',hfss_mod3_5av,hfss_std_mod3_5av,'-r',...
      doy_5av',hs_meas_5av,hs_std_meas_5av,'-.k','alpha','nan','gap');
  set(gca,'XLim',[0 366],'YLim',[-20 135]);
  ylabel('{\itH} (W m^{-2})'); legend('Shrubs','Grass','Trees','Obs')
  text(7,130,'b)','FontSize',12);
    
subplot(2,2,3)
 [hl, hp] = boundedline(doy_5av', hfls_mod_5av,hfls_std_mod_5av,'-b',...
      doy_5av',hfls_mod4_5av,hfls_std_mod4_5av,'-r',...
      doy_5av',le_meas_5av,le_std_meas_5av,'-.k','alpha','nan','gap');
  set(gca,'XLim',[0 366],'YLim',[-10 135]);
  xlabel('DOY');
  ylabel('LE (W m^{-2})'); legend('Shrubs','Shrubs - original \beta','Obs')
  text(7,130,'c)','FontSize',12);

 subplot(2,2,4)
 [hl, hp] = boundedline(doy_5av', hfss_mod_5av,hfss_std_mod_5av,'-b',...
      doy_5av',hfss_mod4_5av,hfss_std_mod4_5av,'-r',...
      doy_5av',hs_meas_5av,hs_std_meas_5av,'-.k','alpha','nan','gap');
  set(gca,'XLim',[0 366],'YLim',[-20 135]);
  xlabel('DOY');
  ylabel('{\itH} (W m^{-2})'); legend('Shrubs','Shrubs - original \beta','Obs')
  text(7,130,'d)','FontSize',12);  
%%
fig = eh_set_fig(fig);
  set(gcf, 'Units', 'Centimeters', 'Position', [0 0 24 12]);
subplot(1,3,1)
  [hl, hp] = boundedline(doy_5av', nep_mod_5av, nep_std_mod_5av,'-b',...
      doy_5av', nep_mod4_5av, nep_std_mod4_5av,'-r',...
      doy_5av',nep_meas_5av,nep_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[-1 3]);
  xlabel('DOY');
  ylabel('NEP (g C m^{-2})'); legend('Shrubs','Shrubs - original \beta','Obs')
  %title('Shrubs')
  text(7,2.8,'a)','FontSize',12);

 subplot(1,3,2)
  [hl, hp] = boundedline(doy_5av', gpp_mod_5av, gpp_std_mod_5av,'-b',...
      doy_5av',gpp_mod4_5av,gpp_std_mod4_5av,'-r',...
      doy_5av',gep_meas_5av,gep_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[-0.5 5]);
  xlabel('DOY');
  ylabel('GPP (g C m^{-2})'); legend('Shrubs','Shrubs - original \beta','Obs')
  text(7,4.7,'b)','FontSize',12);
  
 subplot(1,3,3)
 [hl, hp] = boundedline(doy_5av', ra_mod_5av + rh_mod_5av,ra_std_mod_5av + rh_std_mod_5av,'-b',...
      doy_5av',ra_mod4_5av + rh_mod4_5av,ra_std_mod4_5av + rh_std_mod4_5av,'-r',...
      doy_5av',re_meas_5av,re_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[-0.5 5]);
  xlabel('DOY');
  ylabel('{\itR}_e (g C m^{-2})'); legend('Shrubs','Shrubs - original \beta','Obs')
  text(7,4.7,'c)','FontSize',12);
%%
% VWC including old CEVAP simulation
fig = eh_set_fig(fig);
set(gcf, 'Units', 'Centimeters', 'Position', [0 0 24 18]);
subplot(2,3,1)
  [hl, hp] = boundedline(doy_5av', mrsl_rel1_mod_5av, mrsl_rel1_std_mod_5av,'-b',...
      doy_5av', mrsl_rel1_mod4_5av, mrsl_rel1_std_mod4_5av,'-r',...
     doy_5av',sm_rel1_meas_5av,sm_rel1_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[0 1]);
  set(gca,'XTick',(0:50:366),'Fontsize',12);
  xticklabels({'0','','100','','200','','300',''});
  title('Shrubs');
  ylabel('Relative VWC Layer 1 (0-10cm)'); 
  legend('New \beta','Original \beta','Obs')
  text(5,1.1,'a)','Fontsize',12);
  
  subplot(2,3,2)
  [hl, hp] = boundedline(doy_5av', mrsl_rel1_mod2_5av, mrsl_rel1_std_mod2_5av,'-b',...
     doy_5av',sm_rel1_meas_5av,sm_rel1_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[0 1]);
  set(gca,'XTick',(0:50:366),'Fontsize',12);
  xticklabels({'0','','100','','200','','300',''});
  title('Grass');
subplot(2,3,3)
  [hl, hp] = boundedline(doy_5av', mrsl_rel1_mod3_5av, mrsl_rel1_std_mod3_5av,'-b',...
     doy_5av',sm_rel1_meas_5av,sm_rel1_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[0 1]);
  set(gca,'XTick',(0:50:366),'Fontsize',12);
  xticklabels({'0','','100','','200','','300',''});
  title('Trees');
subplot(2,3,4)
  [hl, hp] = boundedline(doy_5av', mrsl_rel2_mod_5av, mrsl_rel2_std_mod_5av,'-b',...
      doy_5av', mrsl_rel2_mod4_5av, mrsl_rel2_std_mod4_5av,'-r',...
     doy_5av',sm_rel3_meas_5av,sm_rel3_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[0 1]);
  set(gca,'XTick',(0:50:366),'Fontsize',12);
  xticklabels({'0','','100','','200','','300',''});
  xlabel('DOY');
  ylabel('Relative VWC Layer 2 (10-20cm)'); 
  text(5,1.1,'b)','Fontsize',12);
  
  subplot(2,3,5)
  [hl, hp] = boundedline(doy_5av', mrsl_rel2_mod2_5av, mrsl_rel2_std_mod2_5av,'-b',...
     doy_5av',sm_rel3_meas_5av,sm_rel3_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[0 1]);
  set(gca,'XTick',(0:50:366),'Fontsize',12);
  xticklabels({'0','','100','','200','','300',''});
  xlabel('DOY');
subplot(2,3,6)
  [hl, hp] = boundedline(doy_5av', mrsl_rel2_mod3_5av, mrsl_rel2_std_mod3_5av,'-b',...
     doy_5av',sm_rel3_meas_5av,sm_rel3_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[0 1]);
  set(gca,'XTick',(0:50:366),'Fontsize',12);
  xticklabels({'0','','100','','200','','300',''});
  xlabel('DOY');
%%
fig = eh_set_fig(fig);
set(gcf, 'Units', 'Centimeters', 'Position', [0 0 24 10]);
ha = tight_subplot(1, 3, [.08 .08], [.15 .01], [.08 .02]); %[verDist horDist] [bottom top] [left right]
subplot(1,3,1)
  [hl, hp] = boundedline(doy_5av', evap_mod_5av.*86400, evap_std_mod_5av.*86400,'-b',...
     doy_5av',tran_mod_5av.*86400,tran_std_mod_5av.*86400,'-y',...
     doy_5av',groundevap_mod_5av.*86400,groundevap_std_mod_5av.*86400,'-r',...
     doy_5av',canopyevap_mod_5av.*86400,canopyevap_std_mod_5av.*86400,'-m',...
     doy_5av',et_meas_5av,et_std_meas_5av,'-.k','alpha','nan','gap');
set(gca,'XLim',[0 366],'YLim',[-0.5 4]);
title('Shrubs');
xlabel('DOY');
ylabel('Water fluxes (mm d^{-1})'); 
legend('ET','{\itT}','Ground {\itE}','Canopy {\itE}','Obs ET')
subplot(1,3,2)
  [hl, hp] = boundedline(doy_5av', evap_mod2_5av.*86400, evap_std_mod2_5av.*86400,...
     '-b',doy_5av',tran_mod2_5av.*86400,tran_std_mod2_5av.*86400,'-y',...
     doy_5av',groundevap_mod2_5av.*86400,groundevap_std_mod2_5av.*86400,'-r',...
     doy_5av',canopyevap_mod2_5av.*86400,canopyevap_std_mod2_5av.*86400,'-m',...
     doy_5av',et_meas_5av,et_std_meas_5av,'-.k','alpha','nan','gap');
set(gca,'XLim',[0 366],'YLim',[-0.5 4]);
title('Grass');
xlabel('DOY');
legend('ET','{\itT}','Ground {\itE}','Canopy {\itE}','Obs ET')

 subplot(1,3,3)
  [hl, hp] = boundedline(doy_5av', evap_mod3_5av.*86400, evap_std_mod3_5av.*86400,...
     '-b',doy_5av',tran_mod3_5av.*86400,tran_std_mod3_5av.*86400,'-y',...
     doy_5av',groundevap_mod3_5av.*86400,groundevap_std_mod3_5av.*86400,'-r',...
     doy_5av',canopyevap_mod3_5av.*86400,canopyevap_std_mod3_5av.*86400,'-m',...
     doy_5av',et_meas_5av,et_std_meas_5av,'-.k','alpha','nan','gap');
 %doy_5av',et_meas_5av.*86400,et_std_meas_5av.*86400,'-.k','alpha','nan','gap');
set(gca,'XLim',[0 366],'YLim',[-0.5 4]);
title('Trees');
xlabel('DOY');
legend('ET','{\itT}','Ground {\itE}','Canopy {\itE}','Obs ET')    
%%
% soil temperatures for shrub and old cevap run

fig = eh_set_fig(fig);
set(gcf, 'Units', 'Centimeters', 'Position', [0 0 24 16]);
subplot(1,3,1)
  [hl, hp] = boundedline(doy_5av', tsl1_mod_5av, tsl1_std_mod_5av,'-b',...
      doy_5av', tsl1_mod4_5av, tsl1_std_mod4_5av,'-r',...
     doy_5av',Tsoil1_meas_5av,Tsoil1_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[-22 15]);
  xlabel('DOY');
  title('Layer 1 (0-10cm)')
  ylabel('{\itT}_s (^oC)');
  legend('Shrubs','Obs')
subplot(1,3,2)
  [hl, hp] = boundedline(doy_5av', tsl3_mod_5av, tsl3_std_mod_5av,'-b',...
      doy_5av', tsl3_mod4_5av, tsl3_std_mod4_5av,'-r',...
     doy_5av',Tsoil3_meas_5av,Tsoil3_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[-22 15]);
  xlabel('DOY');
  title('Layer 3 (20-30cm)')
  legend('Shrubs','Obs')
subplot(1,3,3)
  [hl, hp] = boundedline(doy_5av', tsl6_mod_5av, tsl6_std_mod_5av,'-b',...
      doy_5av', tsl6_mod4_5av, tsl6_std_mod4_5av,'-r',...
     doy_5av',Tsoil6_meas_5av,Tsoil6_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[-22 15]);
  xlabel('DOY');
  title('Layer 6 (50-60cm)') 
  legend('Shrubs','Shrubs - original \beta','Obs')
%%
% soil temperatures at different depths for shrub run

fig = eh_set_fig(fig);
set(gcf, 'Units', 'Centimeters', 'Position', [0 0 24 14]);
subplot(1,3,1)
  [hl, hp] = boundedline(doy_5av', tsl1_mod_5av, tsl1_std_mod_5av,'-b',...
      doy_5av', tsl1_mod2_5av, tsl1_std_mod2_5av,'-y',... 
      doy_5av', tsl1_mod3_5av, tsl1_std_mod3_5av,'-r',...
      doy_5av',Tsoil1_meas_5av,Tsoil1_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[-22 15]);
  xlabel('DOY');
  title('Layer 1 (0-10cm)')
  ylabel('{\itT}_s (^oC)');
  legend('Shrubs','Grass','Trees','Obs')
subplot(1,3,2)
  [hl, hp] = boundedline(doy_5av', tsl3_mod_5av, tsl3_std_mod_5av,'-b',...
      doy_5av', tsl3_mod2_5av, tsl3_std_mod2_5av,'-y',... 
      doy_5av', tsl3_mod3_5av, tsl3_std_mod3_5av,'-r',...
      doy_5av',Tsoil3_meas_5av,Tsoil3_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[-22 15]);
  xlabel('DOY');
  title('Layer 3 (20-30cm)')
  legend('Shrubs','Grass','Trees','Obs')
subplot(1,3,3)
  [hl, hp] = boundedline(doy_5av', tsl6_mod_5av, tsl6_std_mod_5av,'-b',...
      doy_5av', tsl6_mod2_5av, tsl6_std_mod2_5av,'-y',... 
      doy_5av', tsl6_mod3_5av, tsl6_std_mod3_5av,'-r',...
      doy_5av',Tsoil6_meas_5av,Tsoil6_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[-22 15]);
  xlabel('DOY');
  title('Layer 6 (50-60cm)')
  legend('Shrubs','Grass','Trees','Obs')
%%
% mod only lines 

fig = eh_set_fig(fig);
set(gcf,'color','w');
h1 = plot(t1, cumsum(nep_daily.*1000.*60*60*24),'Color','b','LineWidth',1.1);hold on;
h2 = plot(t2, cumsum(nep_daily2.*1000.*60*60*24),'Color',[1 215/255 0],'LineWidth',1.1);
h3 = plot(t2, cumsum(nep_daily3.*1000.*60*60*24),'r','LineWidth',1.1);
datetick('x','YYYY');xlabel('Year');
yline(0,'--');
xline(t1(2192),':');
ylabel('Cumulative NEP (g C m^{-2})'); 
legend([h1 h2 h3],'Shrubs','Grass','Trees','Location','Southwest')
%%
fig = eh_set_fig(fig);  
[hl, hp] = boundedline(doy_5av', rn_mod_5av,rn_std_mod_5av,'-b',...
      doy_5av', rn_mod2_5av,rn_std_mod2_5av,'-y',...
      doy_5av', rn_mod3_5av,rn_std_mod3_5av,'-r',...
      doy_5av',rn_meas_5av,rn_std_meas_5av,'-.k','alpha');
set(gca,'XLim',[0 366]);
xlabel('DOY');
ylabel('{\itR}_n (W m^{-2})'); 
legend('Shrubs','Grass','Trees','Obs','Location','Northeast')
%%
  fig = eh_set_fig(fig);
set(gcf, 'Units', 'Centimeters', 'Position', [0 0 24 14]);
subplot(1,2,1)
  [hl, hp] = boundedline(doy_5av', ra_mod_5av,ra_std_mod_5av,'-b',...
      doy_5av', ra_mod2_5av,ra_std_mod2_5av,'-y',...
      doy_5av', ra_mod3_5av,ra_std_mod3_5av,'-r','alpha');
  set(gca,'XLim',[0 366],'YLim',[0 1.7]);
  xlabel('DOY');
  %title('Shrubs');
  ylabel('{\itR}_a (g C m^{-2})'); legend('Shrubs','Grass','Trees')
  subplot(1,2,2)
  [hl, hp] = boundedline(doy_5av', rh_mod_5av,rh_std_mod_5av,'-b',...
      doy_5av', rh_mod2_5av,rh_std_mod2_5av,'-y',...
      doy_5av', rh_mod3_5av,rh_std_mod3_5av,'-r','alpha');
  set(gca,'XLim',[0 366],'YLim',[0 1.7]);
  xlabel('DOY');
  ylabel('{\itR}_h (g C m^{-2})'); legend('Shrubs','Grass','Trees')
%%
fig = eh_set_fig(fig);
  set(gcf, 'Units', 'Centimeters', 'Position', [0 0 24 18]);
subplot(2,2,1)
   [hl, hp] = boundedline(doy_5av', nep_mod_5av,nep_std_mod_5av,'-b',...
      'alpha');
 [hl2, hp2] = boundedline(doy_5av',nep_mod2_5av,nep_std_mod2_5av,'-y',...
      'alpha');
  [hl3, hp3] = boundedline(doy_5av',nep_mod3_5av,nep_std_mod3_5av,'-r',...
      'alpha');
  [hl4, hp4] = boundedline(doy_5av',nep_meas_5av,nep_std_meas_5av,'-.k',...
      'alpha');
  hold on;hl5 = plot(doy_5av,nep_chambers_mean_year_5av,'-','Color',[0.5 0.5 0.5]);
  set(gca,'XLim',[0 366],'YLim',[-1 5]);
  xlabel('DOY');
  ylabel('NEP (g C m^{-2})'); legend([hp hp2 hp3 hp4 hl5],{'Shrubs','Grass','Trees','Obs EC','Chambers'})
  text(7,4.7,'a)','FontSize',12);

 subplot(2,2,2)
  [hl, hp] = boundedline(doy_5av', gpp_mod_5av, gpp_std_mod_5av,'-b',...
      doy_5av',gpp_mod2_5av,gpp_std_mod2_5av,'-y',...
      doy_5av',gpp_mod3_5av,gpp_std_mod3_5av,'-r',...
      doy_5av',gep_meas_5av,gep_std_meas_5av,'-.k','alpha');
  set(gca,'XLim',[0 366],'YLim',[-1 5]);
  xlabel('DOY');
  ylabel('GPP (g C m^{-2})'); legend('Shrubs','Grass','Trees','Obs')
  text(7,4.7,'b)','FontSize',12);
  
 subplot(2,2,3)
 [hl, hp] = boundedline(doy_5av', ra_mod_5av + rh_mod_5av,ra_std_mod_5av + rh_std_mod_5av,'-b',...
      'alpha');
 [hl2, hp2] = boundedline(doy_5av',ra_mod2_5av + rh_mod2_5av,ra_std_mod2_5av + rh_std_mod2_5av,'-y',...
      'alpha');
  [hl3, hp3] = boundedline(doy_5av',ra_mod3_5av + rh_mod3_5av,ra_std_mod3_5av + rh_std_mod3_5av,'-r',...
      'alpha');
  [hl4, hp4] = boundedline(doy_5av',re_meas_5av,re_std_meas_5av,'-.k',...
      'alpha');
  hold on;hl5 = plot(doy_5av,chambers_mean_year_5av,'-','Color',[0.5 0.5 0.5]);
  set(gca,'XLim',[0 366],'YLim',[-1 5]);
  xlabel('DOY');
  ylabel('{\itR}_e (g C m^{-2})'); 
  legend([hp hp2 hp3 hp4 hl5],{'Shrubs','Grass','Trees','Obs EC','Chambers'})
  text(7,4.7,'c)','FontSize',12);
  
  subplot(2,2,4)
 [hl, hp] = boundedline(doy_5av', ra_mod_5av + rh_mod_5av,ra_std_mod_5av + rh_std_mod_5av,'-b',...
      'alpha');
 [hl2, hp2] = boundedline(doy_5av',ra_mod_5av,ra_std_mod_5av,'-r',...
      'alpha');
  [hl3, hp3] = boundedline(doy_5av',rh_mod_5av,rh_std_mod_5av,'-m',...
      'alpha');
  [hl4, hp4] = boundedline(doy_5av',re_meas_5av,re_std_meas_5av,'-.k',...
      'alpha');
  hold on;hl5 = plot(doy_5av,chambers_mean_year_5av,'-','Color',[0.5 0.5 0.5]);
  set(gca,'XLim',[0 366],'YLim',[-1 5]);
  xlabel('DOY');
  ylabel('{\itR} (g C m^{-2})'); 
  legend([hp hp2 hp3 hp4 hl5],'{\itR}_e shrubs','{\itR}_a shrubs','{\itR}_h shrubs',...
      '{\itR}_e obs EC','{\itR}_e Chambers')
  text(7,4.7,'d)','FontSize',12);