Contents
versionstring = '1.25';
Configuration and constants
cd('/Users/patrick/Documents/Work/html/Eddycalc');
DATAPATH = '/Users/patrick/Desktop';
FILELENGTH = 30;
AVGTIME = 30;
SAMPLINGRATE = 20;
COORDROT = 'DR';
K = [0.3571 -0.0229 0.9338];
MAXLAG_IRGA = 20;
MEANLAG_CO2 = 6;
MEANLAG_H2O = 6;
DETR = false;
R = 8.31451;
Z = 49;
D = 18;
WINDOFF = 4;
select mat file(s) to load
[files,path] = uigetfile('.mat','Select data file(s)','MultiSelect','on',DATAPATH);
if iscell(files)
files = sort(files);
else
files = num2cell(files,2);
end
n = numel(files);
Eddy covariance data processing
disp(['Eddycalc (',versionstring,') © 2008-2011 Patrick Sturm']);
tic;
fields = {'date';'blocksize';'angle_uw';'wind_dir';'lag_co2';'r_wco2';...
'lag_h2o';'r_wh2o';'mean_Tv';'mean_h2o';'mean_co2';'mean_P';...
'mean_u';'mean_v';'mean_w';'mean_T';'mean_Tin';'mean_Tout';'mean_Tcell';...
'mean_h2o_mr';'mean_co2_mr';'covar_wTv';'covar_wh2o';'covar_wco2';...
'covar_uw';'covar_vw';'covar_wT';'covar_wh2o_mr';'covar_wco2_mr';...
'ustar';'sst_wTv';'sst_wh2o';'sst_wco2';'sst_uw';'var_u2';'var_v2';...
'var_w2';'var_h2o';'var_co2';'var_P';'var_co2_mr';'var_h2o_mr';...
'covar_wP';'stability_parameter';'itc';'qcqa_co2';...
'covar_wco2_wpl';'AGC';'PLL_flag';'Detector_flag';'Chopper_flag';...
'DiffPress_flag';'AuxInput_flag';'Tinlet_flag';'Toutlet_flag'};
num_col = numel(fields);
for i=1:num_col
dataout.(fields{i}) = NaN(n*FILELENGTH/AVGTIME,1);
end
rp = 1;
for i = 1:n
if files{i}==0; continue; end
load(fullfile(path,char(files(i))))
startdate = datenum(strrep(strrep(char(files(i)),'.mat',''),...
'.b',''),'yyyymmddHHMM');
disp([num2str(i),'/',num2str(n),', ',datestr(startdate)])
temp = dec2bin(data.DiagnosticValue,16);
data.AGC = bin2dec(temp(:,13:16))*6.25+6.25;
data.PLL_flag = str2num(temp(:,11));
data.Detector_flag = str2num(temp(:,10));
data.Chopper_flag = str2num(temp(:,9));
data.DiffPress_flag = str2num(temp(:,8));
data.AuxInput_flag = str2num(temp(:,7));
data.Tinlet_flag = str2num(temp(:,6));
data.Toutlet_flag = str2num(temp(:,5));
clear temp
blocksize = SAMPLINGRATE*AVGTIME*60;
to = 0:blocksize:SAMPLINGRATE*FILELENGTH*60;
from = 1:blocksize:SAMPLINGRATE*(FILELENGTH-AVGTIME)*60+1;
data.u = data.AuxiliaryInput1;
data.v = data.AuxiliaryInput2;
data.w = data.AuxiliaryInput3;
data.Tv = data.AuxiliaryInput4;
from(length(from)+1) = length(data.u);
for j=1:length(from)-1
block = structfun(@(x)x(to(j)+1:min(from(j+1)...
,to(j)+blocksize)),data,'UniformOutput',false);
if isempty(block.u); continue; end
dataout.date(rp) = startdate+min(to(j)+blocksize,...
length(data.u))/SAMPLINGRATE/86400;
dataout.blocksize(rp) = length(block.u);
[rotatedWind,theta,phi] = ...
rotateWindVector([block.u block.v block.w],COORDROT,K);
block.u = rotatedWind(:,1);
block.v = rotatedWind(:,2);
block.w = rotatedWind(:,3);
dataout.angle_uw(rp) = phi/pi*180;
if theta<0; theta = theta+2*pi; end
theta = theta/pi*180-WINDOFF;
dataout.wind_dir(rp) = theta;
lag_co2 = findLag(MAXLAG_IRGA,block.CO2dry,block.w);
lag_h2o = findLag(MAXLAG_IRGA,block.H2Odry,block.w);
if lag_co2==MAXLAG_IRGA; lag_co2 = MEANLAG_CO2; end
if lag_h2o==MAXLAG_IRGA; lag_h2o = MEANLAG_H2O; end
dataout.lag_co2(rp) = lag_co2;
dataout.lag_h2o(rp) = lag_h2o;
clear wind
wind.u_co2 = [NaN(lag_co2,1); rotatedWind(1:end-lag_co2,1)];
wind.v_co2 = [NaN(lag_co2,1); rotatedWind(1:end-lag_co2,2)];
wind.w_co2 = [NaN(lag_co2,1); rotatedWind(1:end-lag_co2,3)];
wind.Tv_co2 = [NaN(lag_co2,1); block.Tv(1:end-lag_co2)];
wind.u_h2o = [NaN(lag_h2o,1); rotatedWind(1:end-lag_h2o,1)];
wind.v_h2o = [NaN(lag_h2o,1); rotatedWind(1:end-lag_h2o,2)];
wind.w_h2o = [NaN(lag_h2o,1); rotatedWind(1:end-lag_h2o,3)];
wind.Tv_h2o = [NaN(lag_h2o,1); block.Tv(1:end-lag_h2o)];
T = wind.Tv_h2o;
for l=1:2
h2o_mr = block.H2Odry*1000;
T = wind.Tv_h2o./(1+0.32*h2o_mr/1e6);
end
co2_mr = block.CO2dry;
dataout.mean_Tv(rp) = nanmean(block.Tv);
dataout.mean_h2o(rp) = nanmean(block.H2O);
dataout.mean_co2(rp) = nanmean(block.CO2);
dataout.mean_P(rp) = nanmean(block.TotalPressure);
dataout.mean_u(rp) = nanmean(block.u);
dataout.mean_v(rp) = nanmean(block.v);
dataout.mean_w(rp) = nanmean(block.w);
dataout.mean_T(rp) = nanmean(T);
dataout.mean_h2o_mr(rp) = nanmean(h2o_mr);
dataout.mean_co2_mr(rp) = nanmean(co2_mr);
dataout.mean_Tin(rp) = nanmean(block.TemperatureIn);
dataout.mean_Tout(rp) = nanmean(block.TemperatureOut);
dataout.mean_Tcell(rp) = nanmean(block.CellTemperature);
dataout.var_u2(rp) = nanvar(block.u);
dataout.var_v2(rp) = nanvar(block.v);
dataout.var_w2(rp) = nanvar(block.w);
dataout.var_h2o(rp) = nanvar(block.H2O);
dataout.var_co2(rp) = nanvar(block.CO2);
dataout.var_P(rp) = nanvar(block.TotalPressure);
dataout.var_h2o_mr(rp) = nanvar(h2o_mr);
dataout.var_co2_mr(rp) = nanvar(co2_mr);
dataout.AGC(rp) = nanmean(block.AGC);
dataout.PLL_flag(rp) = min(block.PLL_flag);
dataout.Detector_flag(rp) = min(block.Detector_flag);
dataout.Chopper_flag(rp) = min(block.Chopper_flag);
dataout.DiffPress_flag(rp) = min(block.DiffPress_flag);
dataout.AuxInput_flag(rp) = min(block.AuxInput_flag);
dataout.Tinlet_flag(rp) = min(block.Tinlet_flag);
dataout.Toutlet_flag(rp) = min(block.Toutlet_flag);
if DETR
detrended = nandetrend([block.CO2 block.H2O block.Tv...
block.CO2dry block.H2Odry],1);
block.CO2 = detrended(:,1);
block.H2O = detrended(:,2);
block.Tv = detrended(:,3);
block.CO2dry = detrended(:,4);
block.H2Odry = detrended(:,5);
end
covar1 = cov([block.w block.Tv],1);
dataout.covar_wTv(rp) = covar1(1,2);
covar2 = nancov(block.H2O,wind.w_h2o,1);
dataout.covar_wh2o(rp) = covar2(1,2);
covar3 = nancov(block.CO2,wind.w_co2,1);
dataout.covar_wco2(rp) = covar3(1,2);
covar5 = cov([block.u block.v block.w],1);
dataout.covar_uw(rp) = covar5(1,3);
dataout.covar_vw(rp) = covar5(2,3);
covar = nancov(block.w,T,1);
dataout.covar_wT(rp) = covar(1,2);
covar = nancov([h2o_mr,wind.w_h2o],1);
dataout.covar_wh2o_mr(rp) = covar(1,2);
covar = nancov([co2_mr,wind.w_co2],1);
dataout.covar_wco2_mr(rp) = covar(1,2);
covar = nancov([block.w block.TotalPressure],1);
dataout.covar_wP(rp) = covar(1,2);
dataout.r_wco2(rp) = covar3(1,2)/sqrt(covar3(1,1)*covar3(2,2));
dataout.r_wh2o(rp) = covar2(1,2)/sqrt(covar2(1,1)*covar2(2,2));
dataout.ustar(rp) = (covar5(1,3)^2+covar5(2,3)^2)^(1/4);
n_sub = floor(AVGTIME/5);
covar_sub = NaN(n_sub,5);
for l=1:n_sub
L = floor(length(block.u)/n_sub);
block_sub = structfun(@(x)x(((l-1)*L+1):(l*L)),block,'UniformOutput',false);
wind_sub = structfun(@(x)x(((l-1)*L+1):(l*L)),wind,'UniformOutput',false);
covar = cov([block_sub.w block_sub.Tv],1);
covar_sub(l,1) = covar(1,2);
covar = nancov([block_sub.H2O,wind_sub.w_h2o],1);
covar_sub(l,2) = covar(1,2);
covar = nancov([block_sub.CO2,wind_sub.w_co2],1);
covar_sub(l,3) = covar(1,2);
covar = cov([block_sub.u block_sub.w],1);
covar_sub(l,5) = covar(1,2);
end
dataout.sst_wTv(rp) = (mean(covar_sub(:,1))-covar1(1,2))/covar1(1,2)*100;
dataout.sst_wh2o(rp) = (mean(covar_sub(:,2))-covar2(1,2))/covar2(1,2)*100;
dataout.sst_wco2(rp) = (mean(covar_sub(:,3))-covar3(1,2))/covar3(1,2)*100;
dataout.sst_uw(rp) = (mean(covar_sub(:,5))-covar5(1,3))/covar5(1,3)*100;
rp = rp+1;
end
clear from to block block_sub wind_sub wind T co2_mr h2o_mr detrended rotatedWind
end
L = -abs(dataout.ustar.^3).*dataout.mean_Tv./(0.4*9.81*dataout.covar_wTv);
dataout.stability_parameter = (Z-D)./L;
sigma_uw_theory = 2*abs(dataout.stability_parameter).^0.125;
sigma_uw_theory(abs(dataout.stability_parameter)>1) = 2;
sigma_uw_theory(abs(dataout.stability_parameter)<0.032) = 1.3;
dataout.itc = (sigma_uw_theory-sqrt(dataout.var_w2)./...
dataout.ustar)./sigma_uw_theory*100;
dataout.qcqa_co2 = zeros(length(dataout.itc),1);
dataout.qcqa_co2(abs(dataout.sst_wco2)>30|abs(dataout.itc)>30) = 1;
dataout.qcqa_co2(abs(dataout.sst_wco2)>100|abs(dataout.itc)>100) = 2;
dataout.qcqa_co2(dataout.covar_uw>0) = 2;
[wpl_wh2o wpl_wT,wpl_wP] = calcWPL(dataout.mean_h2o,dataout.covar_wh2o,...
dataout.mean_T,dataout.covar_wT,dataout.mean_co2,dataout.mean_P,...
dataout.covar_wP);
dataout.covar_wco2_wpl = dataout.covar_wco2+wpl_wh2o+wpl_wT+wpl_wP;
dataout_array = reshape(cell2mat(struct2cell(dataout)),length(dataout.date),num_col);
clear blocksize rp i j l startdate fields n_sub L ...
theta phi covar covar1 covar2 covar3 covar4 covar5 covar_sub ...
lag_co2 lag_h2o num_col sigma_uw_theory ...
wpl_wP wpl_wT wpl_wh2o
toc;