function [time, SnDepth, SoilTemp] = read_soil(path,period)
z = dir([path 'gradobs*']);
Archive = [];
for i = 1:length(z)
    x = read_GMD([path z(i).name]);
    hrmn = x(:,4);
    hour = fix(hrmn/100);
    mnt = rem(hrmn,100);
    time = x(:,3) + datenum([2018 1 0]) + hour/24 + mnt/24/60;
    Time = datevec(time);
    Archive = [Archive; time x(:,24:34)];
end
Archive(Archive(:,1) < datenum(period(1,:)) | Archive(:,1) > datenum(period(2,:)) + 1 - 1/24/60,:) = [];
Archive(Archive == -9999.9999) = NaN;

time = Archive(:,1);
SnDepth = Archive(:,2);
SoilTemp = Archive(:,3:end);

% QC
load([path 'SnDepth_qc.mat'])
load([path 'SoilTemp_qc.mat'])
qc_SD(qc_SD(:,1) < datenum(period(1,:)) | qc_SD(:,1) > datenum(period(2,:))+ 1 - 1/24/60,:) = [];
qc_soil(qc_soil(:,1) < datenum(period(1,:)) | qc_soil(:,1) > datenum(period(2,:))+ 1 - 1/24/60,:) = [];

SoilTemp(qc_soil(:,2:11) > 0) = NaN;
SnDepth(qc_SD(:,2) > 0) = NaN;

% convert
SnDepth = 3.6576 - SnDepth/1000;
SoilTemp = SoilTemp' + 273.15;
