cd c:\data\cwf\Dropbox\division\Review_2015\NOAA_PSD_init\OregonProj\byron_data
load('da_PS01_all_PF_rad_soil_corrected');
sst=((-rl_up-0.02*rl_dn)/.98/5.67e-8).^0.25-273.15;
jd=dnum-datenum(2015,12,31);
jdy=24*(jd-floor(jd));
q03=qair3([P_mb Ta03_C rh03 ]);
q10=qair3([P_mb Ta10_C rh10 ]);
%bias correct Ta10 and sst
Ta10_C=Ta10_C-0.25;
sst=sst+1.6;
rs_up=0.9*rs_up;
%%%%%%%%%%%%%%%%%%%

Th03=Ta03_C+.0098*3;
Th10=Ta10_C+.0098*10;
DTs5=diff(Ts5);DTs5=[DTs5(1) DTs5']';
G=rs_dn+rs_up+rl_dn+rl_up-hs03-hl03;
hl03z=hl03; ii=find(co2_ppm<380 | co2_ppm>430);hl03z(ii)=0;
G=rs_dn+rs_up+rl_dn+rl_up-hs03-hl03z;
ik=find(abs(G)>400);G(ik)=0.25*hs03(ik);

%***********   diurnal averages for selected period
daysin=30;%length of interval
jj0=230;% start of interval
kk=find(jd>=jj0 & jd<jj0+daysin);



 [X, YTa03, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),Ta03_C(kk),1,-20, 40);
 [X, YTa10, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),Ta10_C(kk),1,-20, 40);
 [X, Ysst, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),sst(kk),1,-20, 40);
 [X, YTs5, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),Ts5(kk),1,-20, 40);

 [X, Ywspd03, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),wspd03_m_s(kk),1,-20, 40);
 [X, Ywspd10, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),wspd10_m_s(kk),1,-20, 40);
 
 [X, Yhs03, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),hs03(kk),1,-20, 40);
 [X, Yhs10, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),hs10(kk),1,-20, 40);
 [X, Ytau03, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),tauc03(kk),1,-20, 40);
 [X, Ytau10, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),tauc10(kk),1,-20, 40);


 ii=find(co2_ppm>380 & co2_ppm<430 & jd>=jj0 & jd<jj0+daysin);[X, Yhl03, M, SD, SE, SDM, SEM, N] = binave(jdy(ii),hl03(ii),1,-20, 40);

 [X, Yrsu, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),rs_up(kk),1,-20, 40);
 [X, Yrsd, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),rs_dn(kk),1,-20, 40);
 [X, Yrlu, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),rl_up(kk),1,-20, 40);
 [X, Yrld, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),rl_dn(kk),1,-20, 40);
 
 [X, YDT, M, SD, SE, SDM, SEM, N] = binave(jdy(kk),DTs5(kk),1,-20, 40);

 YTh03=YTa03+.0098*3;
 YTh10=YTa10+.0098*10;
 YG=Yrsd+Yrsu+Yrlu+Yrld-Yhs03-Yhl03;
 
figure;
subplot(2,1,1);plot(X,YTh03,X,YTh10,X,Ysst,'--',X,YTs5,'-o');grid;xlabel('Time of Day (hrs, GMT)');ylabel('Temperature (C)');legend('Ta03','Ta10','SST','Ts5');title(['JD ' num2str(jj0) ' to ' num2str(jj0+daysin)]); 
subplot(2,1,2);plot(X,Yhs03,X,Yhl03,X,Yrsd+Yrsu,'--',X,Yrld+Yrlu,'-o');grid;xlabel('Time of Day (hrs, GMT)');ylabel('Flux (W/m^2)');legend('Hs','Hl','Net Solar','Net IR');

h=0.05;%depth of Ts5 in m
cpg=750;%est specific heat of soil J/kg/C
rhog=3E3;%est density soil kg/m^3
delti=1800;%time interval for derivative
k=rhog*cpg*h/delti;
figure;plot(Ysst-YTs5,YG+k*YDT,'o',Ysst-YTs5,9*(Ysst-YTs5));grid;xlabel('SST-Ts5 (C)');ylabel('Ground Flux (W/m^2)');
daysin=15;
n=floor((500-155)/daysin);
jj0=155;avem=[];
for i= 1:n
    jjs=jj0+daysin*(i-1);
    jje=jj0+i*daysin;
    kk=find(jd>jjs & jd<jje &rs_dn>0 & rs_dn<30 & wspd03_m_s>2);
    kk=find(jd>jjs & jd<jje);
    avem(i,:)=nanmean([rs_dn(kk) rs_up(kk) rl_dn(kk) rl_up(kk) hs03(kk) hs10(kk) hl03z(kk) tauc03(kk) tauc10(kk) G(kk) sst(kk) Th03(kk) Th10(kk) q03(kk) q10(kk) wspd03_m_s(kk) wspd10_m_s(kk) jd(kk) Ts5(kk)] );
               %          1          2       3          4        5        6          7         8           9       10     11     12       13      14       15      16               17   
end;

jdm=avem(:,18);
figure;plot(1.25*1e3*2.5e-3*(avem(:,11)-avem(:,12)).*avem(:,17),avem(:,5),'o',[-40 120],[-40 120],'--');grid;xlabel('rhoa*cpa*CH*(sst-Th10+1.6)*U10 (W/m^2)');ylabel('hs03 (W/m^2)');
figure;plot(1.25*4e-3*avem(:,17).^2,avem(:,8),'o',1e-2*[0 35],1e-2*[-0 35],'--');grid;xlabel('rhoa*Cd*U10^2 (N/m^2)');ylabel('tauc03 (N/m^2)');




