% reads in microwve .sci files files and writes out % one datafile with all of the pertinent info in it. % program written by C. Fairall. Modified by D. E. Lane % on 5 March 2001. Added 10-min statistics. disp('Reading ceilo files') clear all; fclose('all'); s='d:\data\ceilo\ceiloraw\*.txt'; dr=dir(s); %lists all files in directory [n m]=size(dr); nx=n;% %# of files in data directory (includes . and ..) %****************** setup read data loop ********** ss='d:\data\ceilo\ceiloraw\'; jamx=1; for ibg = 3:nx %major read loop e=[ss dr(ibg).name]; %gets names of files in directory stux=load(e); %loads all data from files into stux %****** decode date ************ if isempty(stux) % if no data, end the program break; end; [nl ml]=size(stux); %nl=file length - m1 is number of cols ff=dr(ibg).name; %names of files month=str2num(ff(4:5)); %get month day=str2num(ff(1:2)); %get day year=str2num(ff(7:8)); %get year jd0=datenum(2000+year,month,day)-datenum(2000+year-1,12,31); % calculates julian date jd=jd0+((stux(:,5)+stux(:,6)/60+stux(:,7)/3600)/24); % calculate decimal jdate happ(jamx:jamx+nl-1,1:5)=[jd stux(:,1:4)]; % make data array jamx=jamx+nl end; % data line loop happ(:,3)=happ(:,3);%convert ft to meters % plot 15-sec data ik=find(happ(:,3)==0);happ(ik,2)=5; ik=find(happ(:,2)==4 & happ(:,3)<0);happ(ik,2)=5; ill=find(happ(:,2)<5 & happ(:,2)>0); figure(1);plot(happ(:,1),happ(:,2),'o'); xlabel('Julian Day'); ylabel('Number of Cloud Layers'); figure(2);plot(happ(ill,1),happ(ill,3),'.');xlabel('Julian Day');ylabel('Layer-1 Cloud Bottom (m)'); % save 15-sec data ceilo_30sec=happ; % this is the raw data with date save d:\data\ceilo\weller03_ceilo_30s_a.txt ceilo_30sec -ascii -tabs %****************** setup read data loop ********** ny=jamx-1; lam=0:24; % hours jdx=lam/24; % fraction of day jdd=happ(:,1); % date vector jd0=floor(jdd(1)); % first Julian Day of data series jxm=floor(jdd(ny)); % last Julian Day of data series jx=1; % counts number of hours of data kx=1; % counts number of 10-min blocks of data min=0:6; % 10-min parts of an hour jdmx=min/6/24; % 10-min fraction of hour while jd0<=jxm+1 % the first day is l.t.e. last day ii=find(jdd>jd0 & jdd<(jd0+1)); % find all data if isempty(ii) % if no data jd0=jd0+1; % increment the start date else jj=1; % if not empty start first hour calc. while jj<25 % not the end of the day pp=jdd(ii)-floor(jdd(ii)); % fractional part of day ij=find(pp>jdx(jj) & pp(pp(ij(1))+jdmx(kk)) & pp(ij)<(pp(ij(1))+jdmx(kk+1))); % select 10-min section if isempty(kj) fracm(kx)=NaN; hgtm(kx,1:2)=NaN; fracmin(kx)=NaN; fracmino(kx)=NaN; mcloudy(kx)=NaN; hgtmin(kx)=NaN; hgtsm1(kx)=NaN; hgtsm2(kx)=NaN; mindata(kx)=0; clearm(kx)=NaN; onecldm(kx)=NaN; mltcldm(kx)=NaN; obscurem(kx)=NaN; pobsm(kx)=NaN; fracclr(kx)=NaN; else [mindata(kx) c]=size(pp(kj)); % how much data in 10-min section mincloud=find(happ(ii(ij(kj)),2)>0 & happ(ii(ij(kj)),2)<4); % if cloud minobscure=find(happ(ii(ij(kj)),2)==4); % if vertical vis if isempty(mincloud); mcloud(kx)=0; else [mcloud(kx),dx]=size(mincloud); end; if isempty(minobscure) obscurem(kx)=0; else [obscurem(kx),dy]=size(minobscure); end; mcloudy(kx)=mcloud(kx)+obscurem(kx); % total possible cloud fracmin(kx)=mcloud(kx)/mindata(kx); % cloud fraction w/o obscure fracmino(kx)=mcloudy(kx)/mindata(kx); % cloud fraction 2/ obscure cleardata=find(happ(ii(ij(kj)),2)==0); % find clear samples [clearm(kx),dz]=size(cleardata); % number of clear samples fracclr(kx)=clearm(kx)/mindata(kx); % clear fraction monelay=find(happ(ii(ij(kj)),2)==1); % find times with one layer [onecldm(kx),ddx]=size(monelay); % number times one layer mmltcld=find(happ(ii(ij(kj)),2)==2 | happ(ii(ij(kj)),2)==3); [mltcldm(kx),ddy]=size(mmltcld); % number time multiple layers ispobsm=find(happ(ii(ij(kj)),2)==5); [pobsm(kx),ddz]=size(ispobsm); % number of partially obscured times if mcloudy(kx)==0 hgtmin(kx)=NaN; hgtsm1(kx)=NaN; hgtsm2(kx)=NaN; else if mcloudy(kx)>0 allmincover=[mincloud' minobscure']; %all possible cloud samples hgtmin(kx)=median(happ(ii(ij(kj(allmincover))),3)); % take median hgt chgtsm=happ(ii(ij(kj(allmincover))),3); %get hgts all poss. clouds yy=sort(chgtsm); %sort [nny mmy]=size(yy); jj1=max(floor(nny*.16),1); %calculates 1-sigma for height jj2=max(floor(nny*.84),1); hgtsm1(kx)=yy(jj1); hgtsm2(kx)=yy(jj2); end; end; % end if not cloudy end; jdmin(kx)=jd0+jdx(jj)+jdmx(kk); % integrate timestep kk=kk+1; kx=kx+1; end %while kk [alldata(jx) b]=size(pp(ij)); % all data in selected hour iscloud=find(happ(ii(ij),2)>0 & happ(ii(ij),2)<4);%if cloud isobscure=find(happ(ii(ij),2)==4);%if vertical vis if isempty(iscloud) cloud(jx)=0; else [cloud(jx),dx]=size(iscloud); end; if isempty(isobscure) obscure(jx)=0; else [obscure(jx),dy]=size(isobscure); end; cloudy(jx)=cloud(jx)+obscure(jx);%add cloud and 'totally obscured' cases ceilofrac(jx)=cloud(jx)/alldata(jx); % cloud fraction cloud only ceilofraco(jx)=cloudy(jx)/alldata(jx);% cloud fraction cloud+obscure isclr=find(happ(ii(ij),2)==0); % if clear [clr(jx),c]=size(isclr); % number of clear points clrfrac(jx)=clr(jx)/alldata(jx); % fraction of time clear isonecld=find(happ(ii(ij),2)==1); % one cloud layer [onecld(jx) d]=size(isonecld); % number of point one layer ismltcld=find(happ(ii(ij),2)==2 | happ(ii(ij),2)==3); [mltcld(jx) e]=size(ismltcld); % number of points multiple layer ispobs=find(happ(ii(ij),2)==5); % is part obscured [pobs(jx) f]=size(ispobs); % amt part obscured if cloudy(jx)==0 cloudhgt(jx)=NaN; hgts1(jx)=NaN; hgts2(jx)=NaN; else if cloudy(jx)>0 allcover=[iscloud' isobscure']; % all possible cloud cloudhgt(jx)=median(happ(ii(ij(allcover)),3)); %median hgt chgts=happ(ii(ij(allcover)),3); % all hgts poss. cloud y=sort(chgts); [ny my]=size(y); j1=max(floor(ny*.16),1); % calculate +/- 1-sigma hgts j2=max(floor(ny*.84),1); hgts1(jx)=y(j1); hgts2(jx)=y(j2); end; end; % end if not cloudy end;%end if isempty ij jdh(jx)=jd0+jdx(jj); jj=jj+1; jx=jx+1; end;%end while jj jd0=jd0+1 end;%end if isempty ii end;%end while jx % plot rest of data ii=find(cloudhgt==0);cloudhgt(ii)=NaN;ceilofrac(ii)=NaN;hgts1(ii)=NaN;hgts2(ii)=NaN; %jdh=jdh(1:length(ceilofrac)); figure(3);plot(jdh,ceilofrac,'o');hold on;plot(jdmin,fracmin,'r.') xlabel('Julian Day (2004)');ylabel('Hourly (o) and 10-min (.) Ceilo Cloud Fraction'); figure(4);plot(jdh,cloudhgt,'o',jdmin,hgtmin,'.'); xlabel('Julian Day (2004)');ylabel('Ceilo Cloudbase Hts');title('Hourly (o) and 10-min (.)'); figure(5);subplot(2,1,1);plot(jdh,cloudhgt,'.',jdh,hgts1,'o',jdh,hgts2,'x'); hold on;errorbar_dana(jdh,cloudhgt,hgts1,hgts2,'.'); ylabel('Hourly Hts (+/- sigma)'); title('Ceilo Cloud Base Heights : top (hourly) bottom (10-min)'); subplot(2,1,2);plot(jdmin,hgtmin,'.',jdmin,hgtsm1,'o',jdmin,hgtsm2,'x'); hold on;errorbar_dana(jdmin,hgtmin,hgtsm1,hgtsm2,'.');ylabel('10-min Hts (+/- sigma)'); xlabel('Julian Day (2004)'); % save the data in hourly and 10-min files ceilo_hourly=[jdh' alldata' clr' onecld' mltcld' obscure' pobs' clrfrac' ceilofrac' ceilofraco' ... cloudhgt' hgts1' hgts2']; save d:\data\ceilo\neaqs04_ceilo_h_a.txt ceilo_hourly -ascii -tabs ceilo_10min=[jdmin' mindata' clearm' onecldm' mltcldm' obscurem' pobsm' fracclr' fracmin' ... fracmino' hgtmin' hgtsm1' hgtsm2']; save d:\data\ceilo\neaqs04_ceilo_10min_a.txt ceilo_10min -ascii -tabs % clean up memory %clear %pack