ij=find(jdh>jd(1));zb=cloudhgt(ij(1));
ib=floor(zb/25);%cloudbase index

ztop=httb*1e3;
it=floor(ztop/25);

refz=10.^((Ref-180)/10);%convert dBZ to Z reflectivity 

xzz=25*sum(sqrt(refz(ib:it,:)));%Integrate sqrt refl over cloud
q=liq'.*sqrt(refz(ipl,:))./xzz;

ipl=it-1;%select height near cloud top
figure;plot(jd,q);xlabel('Year Day (UTC)');ylabel('Cloud liquid content (g/m^3)');
title(['LWC at z= ' num2str(it*25)]);

nn=100*1e6;%# drops per m^3
reffect=(3*q/4/pi/1e6/nn).^0.333*exp(-0.35^2)*1e6;%effective radius in mirons
figure;plot(jd,reffect);xlabel('Year Day (UTC)');ylabel('r_{eff} (\mum )');
title(['Effective radius at z= ' num2str(it*25)]);
