function [time, ghf] = read_ghf(path, period)
warning off
z = dir([path 'nsasebsE10.b1.*']);
Archive = [];
for i = 1:length(z)
    x = rd_netcdf([path z(i).name]);
    y = attribute_list([path z(i).name]);
    time = double(x.base_time)/60/60/24 + datenum([1970 1 1 0 0 0]) + double(x.time/60/60/24);
    % QC
    ghf = x.surface_soil_heat_flux_avg;
    ghf(x.qc_surface_soil_heat_flux_avg > 0) = NaN;
    % accumulate
    Archive = [Archive; time double(ghf)];
end

Archive(Archive(:,1) < datenum(period(1,:)) | Archive(:,1) > datenum(period(2,:)) + 1 - 30/24/60,:) = [];
clear mex

% interpret
time = Archive(:,1);
ghf = Archive(:,2);
