%read 10-min SHEBA "stats" data
clear all
load D:\SHEBA\data\proc_stat_clb_rdo10_10min.txt
indata=proc_stat_clb_rdo10_10min;
nz=length(indata);
clear proc_stat_clb_rdo10_10min
jd10=indata(:,1);%decimal JD
jd1097=indata(:,2);
for k=1:5
  z(1:nz,k)=indata(:,k+2);%instrument heights
  ws10(1:nz,k)=indata(:,k+7);%streamwise winds
  wd10(1:nz,k)=indata(:,k+12);%wind direction
end
tisn10=indata(:,18);%snow-ice interface T
tsnw10=indata(:,19);%snow surface T (thermistor)
tsr10=indata(:,20);%PIR radiative surface T
tbr10=indata(:,21);%Barnes surface T in Cent, uncorrected for emissivity/reflection
Vbr10=indata(:,22);%Barnes voltage
tge10=indata(:,23);%General Eastern air T
tdge10=indata(:,24);%General Eastern air Td
for k=1:5
  t10(1:nz,k)=indata(:,k+24);%Vaisala T
  tsn10(1:nz,k)=indata(:,k+29);%sonic temperatures (uncalibrated)
  rhwv10(1:nz,k)=indata(:,k+34);%Vaisala RH (% wrt water)
  ww10(1:nz,k)=indata(:,k+51);;%sonic vertical velocities, m/s
end
Rldx10=indata(:,40);%thermopile longwave down (incoming)
Rlux10=indata(:,41);%thermopile longwave up (outgoing)
Tcdx10=indata(:,42);%Tcase longwave down 
Tcux10=indata(:,43);%Tcase longwave up
Tddx10=indata(:,44);%Tdome longwave down
Tdux10=indata(:,45);%Tdome longwave up
lwd10=indata(:,46);%longwave down (incoming)
lwu10=indata(:,47);%longwave up (outgoing)
swd10=indata(:,48);%solar down (incoming)
swu10=indata(:,49);%solar up (outgoing)
org10=indata(:,50);%optical raingauge, mm/h
ncr10=indata(:,51);%ncar tipping raingauge, mm/h

%edit data
t10(find(abs(t10)>100))=NaN;
tisn10(find(abs(tisn10)>100))=NaN;
tsnw10(find(abs(tsnw10)>100))=NaN;
tsr10(find(abs(tsr10)>100))=NaN;
tbr10(find(abs(tbr10)>100))=NaN;
tge10(find(abs(tge10)>100))=NaN;
tdge10(find(abs(tdge10)>100))=NaN;
ws10(find(abs(ws10)>100))=NaN;
wd10(find(abs(wd10)>400))=NaN;
rhwv10(find(abs(rhwv10)>200))=NaN;
ww10(find(abs(ww10)>10))=NaN;
z(find(abs(z)>20))=NaN;
lwd10(abs(lwd10)>900)=NaN;
lwu10(abs(lwu10)>900)=NaN;
swd10(abs(swd10)>900)=NaN;
swu10(abs(swu10)>900)=NaN;
org10(abs(org10)>900)=NaN;
ncr10(abs(ncr10)>900)=NaN;

for k=1:5
  zave(k)=nanomean(z(:,k));
end

%plot data
fn=0;
fn=fn+1;
figure(fn)
subplot(3,1,1),plot(jd1097,t10(:,1),'r-',jd1097,t10(:,2),'b-',jd1097,t10(:,3),'m-',jd1097,t10(:,4),'c-',jd1097,t10(:,5),'g-',[jd1097(1) jd1097(end)],[0 0],'k--')
legend([num2str(round(zave(1)*10)/10) ' m'],[num2str(round(zave(2)*10)/10) ' m'],[num2str(round(zave(3)*10)/10) ' m'],[num2str(round(zave(4)*10)/10) ' m'],[num2str(round(zave(5)*10)/10) ' m'])
xlabel('Year Day wrt Jan 1 1997');ylabel('Air Temperature (deg C)')
ylim([-45 5])

subplot(3,1,2),plot(jd1097,tsr10,'r-',jd1097,tsnw10,'b-',jd1097,tbr10,'m-',jd1097,tisn10,'c-',[jd1097(1) jd1097(end)],[0 0],'k--')
legend(['T_s_f_c PIR'],['T_s_f_c thrmstr'],['T_s_f_c Barnes'],['T snow-ice intrfc'])
xlabel('Year Day wrt Jan 1 1997');ylabel('Temperature (deg C)')
ylim([-45 5])

subplot(3,1,3),plot(jd1097,rhwv10(:,1),'r-',jd1097,rhwv10(:,2),'b-',jd1097,rhwv10(:,3),'m-',jd1097,rhwv10(:,4),'c-',jd1097,rhwv10(:,5),'g-',[jd1097(1) jd1097(end)],[100 100],'k--')
legend([num2str(round(zave(1)*10)/10) ' m'],[num2str(round(zave(2)*10)/10) ' m'],[num2str(round(zave(3)*10)/10) ' m'],[num2str(round(zave(4)*10)/10) ' m'],[num2str(round(zave(5)*10)/10) ' m'])
xlabel('Year Day wrt Jan 1 1997');ylabel('RH wrt water (%)')
ylim([55 105])

fn=fn+1;
figure(fn)
subplot(3,1,1),plot(jd1097,ws10(:,1),'r-',jd1097,ws10(:,2),'b-',jd1097,ws10(:,3),'m-',jd1097,ws10(:,4),'c-',jd1097,ws10(:,5),'g-')
legend([num2str(round(zave(1)*10)/10) ' m'],[num2str(round(zave(2)*10)/10) ' m'],[num2str(round(zave(3)*10)/10) ' m'],[num2str(round(zave(4)*10)/10) ' m'],[num2str(round(zave(5)*10)/10) ' m'])
xlabel('Year Day wrt Jan 1 1997');ylabel('Wind Speed (m/s)')
ylim([0 30])
subplot(3,1,2),plot(jd1097,wd10(:,1),'r-',jd1097,wd10(:,2),'b-',jd1097,wd10(:,3),'m-',jd1097,wd10(:,4),'c-',jd1097,wd10(:,5),'g-')
legend([num2str(round(zave(1)*10)/10) ' m'],[num2str(round(zave(2)*10)/10) ' m'],[num2str(round(zave(3)*10)/10) ' m'],[num2str(round(zave(4)*10)/10) ' m'],[num2str(round(zave(5)*10)/10) ' m'])
xlabel('Year Day wrt Jan 1 1997');ylabel('Wind Direction (deg)')
ylim([0 360])
subplot(3,1,3),plot(jd1097,ww10(:,1),'r-',jd1097,ww10(:,2),'b-',jd1097,ww10(:,3),'m-',jd1097,ww10(:,4),'c-',jd1097,ww10(:,5),'g-',[jd1097(1) jd1097(end)],[0 0],'k--')
legend([num2str(round(zave(1)*10)/10) ' m'],[num2str(round(zave(2)*10)/10) ' m'],[num2str(round(zave(3)*10)/10) ' m'],[num2str(round(zave(4)*10)/10) ' m'],[num2str(round(zave(5)*10)/10) ' m'])
xlabel('Year Day wrt Jan 1 1997');ylabel('Vertical velocity (m/s)')
ylim([-2 2])

fn=fn+1;
figure(fn)
subplot(3,1,1),plot(jd1097,lwd10,'r-',jd1097,lwu10,'b-')
legend(['LW_d'],['LW_u'])
xlabel('Year Day wrt Jan 1 1997');ylabel('W/m^2')
ylim([110 340])
subplot(3,1,2),plot(jd1097,swd10,'r-',jd1097,swu10,'b-')
legend(['SW_d'],['SW_u'])
xlabel('Year Day wrt Jan 1 1997');ylabel('W/m^2')
ylim([0 820])
subplot(3,1,3),plot(jd1097,org10,'r-',jd1097,ncr10,'b-')
legend(['optical raingauge'],['ncar tipping bucket'])
xlabel('Year Day wrt Jan 1 1997');ylabel('mm/h')
ylim([0 10])

