%run conv_cor1 first
%   eyes(ibg,:)=[jt oph tilt jm phirel];
Le=2.5e6;
cp=1e3;
clear dtx Hix cenb;
ix=1;
for i=1:nx
   if eyes(i,1)<1
      if eyes(i,2)<18
         if eyes(i,3)<11
            if eyes(i,4)<5
               if eyes(i,5)<90 | eyes(i,5)>270
                  dtx(ix,:)=dt(i,:);
                  %stullx(ix,:)=stull(i,:);
                  if abs(dt(i,12)-Hi(i,2))>.5*dt(i,12)
                     Hix(ix,:)=dt(i,12);
                  else                  
                  	Hix(ix,:)=Hi(i,:);%Hi_j(i,:);
                  end;
                  cenb(ix)=dt(i,31);
                  ix=ix+1;
               end;
            end;
         end;
      end;
   end;
end;
zu=dtx(:,38);

dT=dtx(:,4)-dtx(:,5)+.0098*zu;
dqer=dtx(:,22);
dq=(dtx(:,6)-dtx(:,7))*1e-3-dqer;
q=dtx(:,7)*1e-3;
t=dtx(:,5)+273.16;

rho=1013*100./t/287./(1+0.61*q);
bo=cp/Le*dT./dq;
wq=(dtx(:,36)+Hix(:,2))/2./rho/Le;
wtb=dtx(:,12)./rho/cp;
wqb=dtx(:,12)./rho/Le;
wstr=(9.8./t.*(wtb+0.61.*t.*wqb)*600).^.333;
L=dtx(:,17);
zu=dtx(:,38);
u=dtx(:,2);
zo=dtx(:,14);
zoq=dtx(:,16);
us2=dtx(:,13)./rho;

ce=wq./sqrt(u.^2+1.25*wstr.^2)./dq;
cen=ce./(1+psiug(zu./L)./log(zu./zo))./(1+psitg(zu./L)./log(zu./zoq));
ce10=.16./(log(10./zu)+.16./cen./log(zu./zo))./log(10./zo);
figure;plot(u,ce10,'.',u,cenb);axis([ 0 12 0.0005 .0015]);
gf=(1+1.61*q).*(1+dT./t.*q./dq);
ceq=ce10.*gf;
figure;plot(u,ceq,'.',u,cenb);axis([ 0 12 0.0005 .0015]);xlabel('Wind Speed (m/s)');ylabel('C_{en}');
z0q=10./exp(.16./ceq./log(10./zo));

