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.5;t=Ta10_C;
sst=sst+1.6;
ssq=qair3([P_mb sst 100*ones(1,length(sst))' ]);
Q=q03/1e3;
    tdk=273.15;
    kon=0.4;
    Rgas = 287.1;
     Le   = (2.501)*1e6;
     cpa  = 1004.67;
     cpv  = cpa*(1+0.84*Q);
     rhoa = P_mb*100./(Rgas*(Ta10_C+tdk).*(1+0.61*Q));
     visa = 1.326e-5*(1+6.542e-3.*t+8.301e-6*t.^2-4.84e-9*t.^3);
   
 Th03=Ta03_C+.0098*3;
 Th10=Ta10_C+.0098*10;
 U=wspd10_m_s;
 z0=5e-3;
 z0t=1e-4;
  Dth=(sst-Th10);
  Dq=(ssq-q10);
  Sm0=2*Sm5-Sm10;
  rhs=Sm10/0.3;rhs=min(1,rhs);ii=find(jd<202);rhs(ii)=0.2;
  ssq2=qair3([P_mb Ts10 100*ones(1,length(sst))' ]);

  Dq2=(ssq2-q10);
  zet=0;
  for i=1:5
  usr=U*kon./(log(10./z0)-psiu_30(zet));
  tsr=-Dth*kon./(log(10./z0t)-psit_30(zet));
  qsr=-Dq*kon./(log(10./z0t)-psit_30(zet));
  qsr2=-Dq2*kon./(log(10./z0t)-psit_30(zet));
  zet=kon*10*9.83./(t+tdk).*(tsr)./usr.^2;
  ii=find(zet>20);zet(ii)=20;
  ii=find(zet<-20);zet(ii)=-20;
  end;
  
  taub=rhoa.*usr.*usr;  % wind stress
  hsb=-rhoa.*cpa.*usr.*tsr;     % sensible heat flux
  hlb=-rhoa.*Le.*usr.*qsr/1e3;      % latent heat flux
  hlb2=-rhoa.*Le.*usr.*qsr2/1e3;      % latent heat flux

 [X, YTa03, M, SD, SE, SDM, SEM, N] = binave(jdy,Ta03_C,1,-20, 40);
 [X, YTa10, M, SD, SE, SDM, SEM, N] = binave(jdy,Ta10_C,1,-20, 40);
 YTh03=YTa03+.0098*3;
 YTh10=YTa10+.0098*10;
 [X, Ysst, M, SD, SE, SDM, SEM, N] = binave(jdy,sst,1,-20, 40);
 [X, YTs5, M, SD, SE, SDM, SEM, N] = binave(jdy,Ts5,1,-20, 40);
 [X, YTs10, M, SD, SE, SDM, SEM, N] = binave(jdy,Ts10,1,-20, 40);
 [X, YMs5, M, SD, SE, SDM, SEM, N] = binave(jdy,Sm5,1,-20, 40);
 [X, YMs5, M, SD, SE, SDM, SEM, N] = binave(jdy,Sm10,1,-20, 40);
 
 [X, Ywspd03, M, SD, SE, SDM, SEM, N] = binave(jdy,wspd03_m_s,1,-20, 40);
 [X, Ywspd10, M, SD, SE, SDM, SEM, N] = binave(jdy,wspd10_m_s,1,-20, 40);

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


 ii=find(co2_ppm>380 & co2_ppm<430);[X, Yhl03, M, SD, SE, SDM, SEM, N] = binave(jdy(ii),hl03(ii),1,-20, 40);
 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);

 [X, Yrsu, M, SD, SE, SDM, SEM, N] = binave(jdy,rs_up,1,-20, 40);
 [X, Yrsd, M, SD, SE, SDM, SEM, N] = binave(jdy,rs_dn,1,-20, 40);
 [X, Yrlu, M, SD, SE, SDM, SEM, N] = binave(jdy,rl_up,1,-20, 40);
 [X, Yrld, M, SD, SE, SDM, SEM, N] = binave(jdy,rl_dn,1,-20, 40);
 
 DTs5=diff(Ts5);DTs5=[DTs5(1) DTs5']';
 G=rs_dn+rs_up+rl_dn+rl_up-hs03-hl03;
 [X, YDT, M, SD, SE, SDM, SEM, N] = binave(jdy,DTs5,1,-20, 40);

 
 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');
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) P_mb(kk) Sm0(kk) Sm5(kk) Sm10(kk) hlb(kk) hlb2(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)');

Rnm=avem(:,1)+avem(:,2)+avem(:,3)+avem(:,4);
Gm=avem(:,10);
Hlm=avem(:,7);
Hlpotm=avem(:,24);
Hlb2m=avem(:,25);
Hsm=avem(:,5);
Pm=avem(:,20);
SSTm=avem(:,11);
Sm0m=avem(:,21);
%%%%%   Penmann-Montieth evap equs
qs=qsat([SSTm Pm])/1e3;
qsm=qsat([SSTm-.1 Pm])/1e3;
es=qs.*Pm./(0.62+0.38*qs);
esm=qsm.*Pm./(0.62+0.38*qsm);
delpsy=10*(es-esm);%slope vapor press/T hPA/C
gamm=1006/0.622/2.5e6*Pm;%phychr const hPA/C
figure;plot(jdm,1.26*delpsy./(delpsy+gamm).*(Rnm-Gm),jdm,Hlm);xlabel('Day of Year (2016)');ylabel('Latent Heat flux (W/m^2)');legend('P-M saturated','Measured');






