disp('microwave radiometer cloud liquid bias correction');
%This routine uses the ceilometer data to located clear periods and then
%uses the values of liq in those periods to estimate a bias correction
   %************************  daily averages  ***********	
 clear jdd jdqd liq0dx latd lond liq0
 jj0=floor(jdq(1));
[nn,nx]=size(jdq);

%first, find daily average clear liq
mx=floor(nn/24/6);jj=0;
ji=floor(jdq);ki=find(diff(ji)>0);mx=length(ki);
for i=1:mx
   ii=find(floor(jdq)==jj0);
   jdd(i)=mean(jdy(ii));
   jdqd(i)=mean(jdq(ii));
   jdqd_0(i)=mean(jdq(ii));
   ij=find(ceilof(ii)<.2 & isfinite(liq(ii)));
   if ~isempty(ij)
      liq0dx(i)=median(liq(ii(ij)));
      nl0(i)=length(ij);
  	else
        liq0dx(i)=NaN;
        nl0(i)=0;
   end;
   
   jj=jj+144;
      jj=ii(length(ii))+1;
   jj0=floor(min(jdq(jj),nn));   
end;

%if dkip, do some controls on liq0d
dkip=1;
if dkip
liq0dx(1)=.005;
for i=2:length(liq0dx)
   if isnan(liq0dx(i))
      liq0dx(i)=liq0dx(i-1);
   end;
   if liq0dx(i)>0.03
      liq0dx(i)=liq0dx(i-1);
   end;
end;
end;

%eliminate clear liq outliers from liq0d, 75 g/m^2
for i=1:mx
	ii=find(ceilof<.4 & isfinite(liq) & floor(jdq)==floor(jdqd(i)));
	kk=find(abs(liq(ii)-liq0dx(i))>200/1e4);
	if ~isempty(kk)
      liq(ii(kk))=liq0dx(i);
      nrj(i)=length(kk);
   end;
end;
%interpolate clear time series
ii=find(ceilof<.4 & isfinite(liq) );
liqx=liq(ii);jdx=jdq(ii);
%liq0=interp1(jdx,liqx,jdq);
liq0=interp1(jdqd,liq0dx,jdq);


ii=1:3817;%fall 99
jj=3818:7573;%spring 00
kk=7574:11288;%fall 00
ll=11299:15875;%spring 01
mm=15876:26450;%fall 01
nn=26451:30338;%spr 02
pp=30339:34082;
pq=34083:36926;
pr=36927:41210;

dexi=ii;figure;plot(jdq(dexi),liq(dexi)*1e4,jdq(dexi),liq0(dexi)*1e4,jdqd,liq0dx*1e4,'-o');axis([jdq(dexi(1)) jdq(dexi(length(dexi))) -100 500]);
dexi=pq;figure;plot(jdq(dexi),liq(dexi)*1e4,jdq(dexi),liq0(dexi)*1e4,jdqd,liq0dx*1e4,'-o');axis([jdq(dexi(1)) jdq(dexi(length(dexi))) -100 500]);
dexi=pr;figure;plot(jdq(dexi),liq(dexi)*1e4,jdq(dexi),liq0(dexi)*1e4,jdqd,liq0dx*1e4,'-o');axis([jdq(dexi(1)) jdq(dexi(length(dexi))) -100 500]);

