Contents

% EDDYCALC Eddy covariance data processing toolbox.
%
% Scripts for analysing LI-7200 eddy covariance data.
%
% Raw data files are first saved as MATLAB binary files (*.mat, conversion
% with txt2mat.m) to speed up data processing.
%
% Currently available features are:
% - different time averaging intervals
% - coordinate system rotation (double rotation, triple rotation or
%   planar-fit method)
% - determination of lags between wind and concentration data
% - point-by-point conversion to mixing ratios
% - linear detrending, filtering, despiking
% - QC/QA: stationarity test, intergral turbulence characteristics test
% - calculation of WPL terms (not really needed for LI-7200)
% - spectral analysis
%
% A better documentation is not available yet, sorry.
%
% These scripts are also available in R (translated by Matthias Zeeman and
% developed into a R package (Reddy) by Rebecca Hiller).
%
% Required MATLAB toolboxes: Statistics Toolbox
%                            Signal Processing Toolbox
%
% AUTHORS: Patrick Sturm <patrick.sturm@gmx.ch>
%          Albin Hammerle <albin.hammerle@uibk.ac.at>
% Version: 1.25
% Date: 07-Oct-2011
%
% REFERENCE:
% Sturm P, Eugster W, Knohl A (2012) Eddy covariance measurements of CO2
% isotopologues with a quantum cascade laser absorption spectrometer,
% Agricultural and Forest Meteorology, 152, 73-82,
% doi:10.1016/j.agrformet.2011.09.007

versionstring = '1.25';

% COPYRIGHT 2008-2011 Patrick Sturm
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% For a copy of the GNU General Public License, see
% <http://www.gnu.org/licenses/>.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Configuration and constants

% Paths
cd('/Users/patrick/Documents/Work/html/Eddycalc'); % m-files
DATAPATH = '/Users/patrick/Desktop'; % data files (yyyymmddHHMMSS.mat)

% New file interval (min)
FILELENGTH = 30;
% Averaging interval (min)
AVGTIME = 30; % set AVGTIME <= FILELENGTH !
% sampling frequency (data points per second)
SAMPLINGRATE = 20; % (Hz)

% Coordinate rotation (PF: planar fit, DR: double roation, TR: triple
% rotation:
COORDROT = 'DR';
% Planar fit coordinate rotation: unit vector K of new z-axis
% calculated with getPlanarFitCoeffs.m:
K = [0.3571 -0.0229 0.9338];

% Maximum lag (in data points) for findLag.m:
MAXLAG_IRGA = 20; % Maximum lag for IRGA data
% Set lag to MEANLAG if no maximum is found:
MEANLAG_CO2 = 6; % Mean lag for IRGA CO2 data (system specific)
MEANLAG_H2O = 6; % Mean lag for IRGA H2O data (system specific)

% Detrending (true: linear detrending, false: no detrending):
DETR = false;

% Universal gas constant (J/K/mol):
R = 8.31451;

% Site specific constants
% Laegeren:
Z = 49; % Measuring height (m)
D = 18; % Displacement height (m)
WINDOFF = 4; % wind direction offset of anemometer in degrees

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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']);

% start timer
tic;

% Initialize output variables
% fields of output data
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'}; % fields of output data
num_col   = numel(fields); % number of output fields
for i=1:num_col
    dataout.(fields{i}) = NaN(n*FILELENGTH/AVGTIME,1);
end
rp = 1; % row pointer

% process files sequentially
for i = 1:n

    % load file
    if files{i}==0; continue; end
    load(fullfile(path,char(files(i))))
    startdate = datenum(strrep(strrep(char(files(i)),'.mat',''),...
        '.b',''),'yyyymmddHHMM'); % timestamp of file
    disp([num2str(i),'/',num2str(n),', ',datestr(startdate)])

    % translate dianostic value
    temp = dec2bin(data.DiagnosticValue,16);
    data.AGC = bin2dec(temp(:,13:16))*6.25+6.25; % Automatic Gain Control, < 62 = ok
    data.PLL_flag = str2num(temp(:,11)); % Lock bit, indicates that optical wheel is rotating at the correct rate, 1 = ok
    data.Detector_flag = str2num(temp(:,10)); % Detector temperature, 1 = ok
    data.Chopper_flag = str2num(temp(:,9)); % Chopper wheel temperature, 1 = ok
    data.DiffPress_flag = str2num(temp(:,8)); % Differential pressure sensor, 1 = ok
    data.AuxInput_flag = str2num(temp(:,7)); % Aux input, 1 = ok
    data.Tinlet_flag = str2num(temp(:,6)); % Thermocouple inlet, 1 = ok
    data.Toutlet_flag = str2num(temp(:,5)); % Thermocouple outlet, 1 = ok
    clear temp

    blocksize = SAMPLINGRATE*AVGTIME*60;
    to = 0:blocksize:SAMPLINGRATE*FILELENGTH*60;
    from = 1:blocksize:SAMPLINGRATE*(FILELENGTH-AVGTIME)*60+1;

    % rename auxiliary input variables
    data.u = data.AuxiliaryInput1;
    data.v = data.AuxiliaryInput2;
    data.w = data.AuxiliaryInput3;
    data.Tv = data.AuxiliaryInput4;

    % loop through blocks
    from(length(from)+1) = length(data.u);
    for j=1:length(from)-1

        % get block of data
        block = structfun(@(x)x(to(j)+1:min(from(j+1)...
            ,to(j)+blocksize)),data,'UniformOutput',false);

        if isempty(block.u); continue; end

        % end time of data block
        dataout.date(rp) = startdate+min(to(j)+blocksize,...
            length(data.u))/SAMPLINGRATE/86400;

        % number of data per block
        dataout.blocksize(rp) = length(block.u);

        % rotate wind coordinate system
        [rotatedWind,theta,phi] = ...
            rotateWindVector([block.u block.v block.w],COORDROT,K);
        block.u = rotatedWind(:,1);
        block.v = rotatedWind(:,2);
        block.w = rotatedWind(:,3);

        % Rotation angle of vertical wind
        dataout.angle_uw(rp) = phi/pi*180;

        % wind direction
        if theta<0; theta = theta+2*pi; end
        % in degrees minus wind direction offset of sonic
        theta = theta/pi*180-WINDOFF;
        dataout.wind_dir(rp) = theta;

        % Find lag between co2, h2o, and w using findLag.m
        lag_co2 = findLag(MAXLAG_IRGA,block.CO2dry,block.w);
        lag_h2o = findLag(MAXLAG_IRGA,block.H2Odry,block.w);
        % confine lag to realistic values
        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;

        % shift wind data to match concentration signals
        % Tv is also shifted to match co2 and h2o for conversion to T
        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)];

        % Convert from sonic virtual temperature to air temperature and
        % from h2o molar density to volume mixing ratio
        % (air temperature is then used for the sensible heat flux instead
        % of applying a sonic heat flux correction)
        T = wind.Tv_h2o;
        for l=1:2 % iterative approach to find h2o_mr and T
            % h2o[mmol/m3]*1000*R*T/100*P = h2o[ppm] with P in hPa
%             h2o_mr = block.H2O*R./block.TotalPressure.*T*10; % mmol/m3 to ppm
            h2o_mr = block.H2Odry*1000;
            T      = wind.Tv_h2o./(1+0.32*h2o_mr/1e6); % (Kaimal and Gaynor 1991)
        end

        % Mixing ratios (umol/mol) instead of molar density (mol/m3) for IRGA
        % co2. This should eliminate the need for WPL terms.
        % Not needed for LI-7200!
%         co2_mr = block.CO2/1000./(block.TotalPressure*100./...
%             (R.*wind.Tv_co2)-block.H2O/1000)*1e6;
        co2_mr = block.CO2dry;
%         co2_mr = block.CO2*1000*R.*wind.Tv_co2./(block.TotalPressure*100)./(1-h2o_mr/1e6)

        % calculate mean of blocks
        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);
        % variance of blocks
        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);

        % calculate several diagnostic values (LI7200 manual page 4-21)
        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);

        % detrend
        if DETR
            detrended = nandetrend([block.CO2 block.H2O block.Tv...
                block.CO2dry block.H2Odry],1);
%             temp2 = [block.CO2 block.H2O block.Tv block.CO2dry block.H2Odry];
%             detrended = temp2 - filterLowpass1(temp2,1/SAMPLINGRATE,600);
%             clear temp2
            block.CO2 = detrended(:,1);
            block.H2O = detrended(:,2);
            block.Tv = detrended(:,3);
            block.CO2dry = detrended(:,4);
            block.H2Odry = detrended(:,5);
        end

        % calculate covariances (nancov removes any rows with NaNs,
        % 'pairwise' option is very slow -> do it seperately)
        covar1 = cov([block.w block.Tv],1); % wTv
        dataout.covar_wTv(rp)      = covar1(1,2);
        covar2 = nancov(block.H2O,wind.w_h2o,1); % wh2o
        dataout.covar_wh2o(rp)     = covar2(1,2);
        covar3 = nancov(block.CO2,wind.w_co2,1); % wco2
        dataout.covar_wco2(rp)     = covar3(1,2);
        covar5 = cov([block.u block.v block.w],1); % momentum
        dataout.covar_uw(rp)       = covar5(1,3);
        dataout.covar_vw(rp)       = covar5(2,3);
        covar = nancov(block.w,T,1);    % wT
        dataout.covar_wT(rp)       = covar(1,2);
        covar = nancov([h2o_mr,wind.w_h2o],1); % wh2o_mr
        dataout.covar_wh2o_mr(rp)  = covar(1,2);
        covar = nancov([co2_mr,wind.w_co2],1); % wco2_mr
        dataout.covar_wco2_mr(rp)  = covar(1,2);
        covar = nancov([block.w block.TotalPressure],1); % wP
        dataout.covar_wP(rp)       = covar(1,2);

        % correlation coefficients between w and co2,h2o
        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));

        % friction velocity
        %dataout.ustar(rp) = -1*sign(covar5(1,3))*sqrt(abs(covar5(1,3))); % W. Eugster's definition
        dataout.ustar(rp) = (covar5(1,3)^2+covar5(2,3)^2)^(1/4);

        % quality control/quality assurance (CarboEurope scheme)
        % stationarity test (steady-state turbulence test)
        n_sub = floor(AVGTIME/5); % number of sub-intervals
        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);       % wT_v
            covar_sub(l,1) = covar(1,2);
            covar = nancov([block_sub.H2O,wind_sub.w_h2o],1); % wh2o
            covar_sub(l,2) = covar(1,2);
            covar = nancov([block_sub.CO2,wind_sub.w_co2],1); % wco2
            covar_sub(l,3) = covar(1,2);
            covar = cov([block_sub.u block_sub.w],1);     % momentum
            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

% stability parameter (with virtual instead of potential temperature)
L = -abs(dataout.ustar.^3).*dataout.mean_Tv./(0.4*9.81*dataout.covar_wTv);  % Obukhov length
dataout.stability_parameter = (Z-D)./L;

% Integral turbulence characteristics test (parametrisation after Foken/CarboEurope)
sigma_uw_theory = 2*abs(dataout.stability_parameter).^0.125; % 0.032 <= |z/L| <= 1
sigma_uw_theory(abs(dataout.stability_parameter)>1) = 2; % range |z/L| > 1
sigma_uw_theory(abs(dataout.stability_parameter)<0.032) = 1.3; % range |z/L| < 0.032

dataout.itc = (sigma_uw_theory-sqrt(dataout.var_w2)./...
    dataout.ustar)./sigma_uw_theory*100;

% quality control/quality assurance (CarboEurope scheme)
% Data quality flag
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; % W. Eugster's special criterion

% WPL correction
[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;

% Output data in matrix form
dataout_array = reshape(cell2mat(struct2cell(dataout)),length(dataout.date),num_col);

% clear temporary variables
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;