%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Compute boundary layer top from vertical-beam profiler (60 m gates) % from the height of maximum range-corrected signal-to-noise ratio (SNR) % above 800 m. Output 10-minutely 16-84th percentile-averaged daylong % boundary layer height files. % % Simon de Szoeke 18 December 2007 % Adapted from profiler915_BLHT.m by Dan Wolfe (need to use Dan's version % for different radar modes, i.e. multiple elevation angles) % and from pick_mbl_snrcr.pro by Bill Otto. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; close all; read_parameters_Stratus %jd=input('Input Jd for raw file to plot (ex: rbc07291.raw, enter 291): ', 's'); for yday=340:355%:355; % -357? loop year day of files to read %yday=340; fn=sprintf('rbc04%3d.r*w',yday); % read .raw or .rew adir=dir([way_raw_data_profiler fn]); cbfn = way_proc_data_ceilo; datev=datevec(datenum(2004,0,0,0,0,0)+yday); cbfn1=sprintf('%03d_%02d_%02d_%02d.txt',yday,datev(3:-1:2),datev(1)-2000); lf1 = length(way_proc_data_ceilo); lf2 = length(cbfn1); cbfn(lf1+1:lf1+lf2)=cbfn1; % rbclear 0 0 13 91500 %5 75 45 7 123 0 0 26 %61 700 196 101 55 110 30 128 % -4.308 13.240 62.462 3.114 28.150 0 % -4.622 12.690 51.212 4.098 17.450 0 fid=fopen([way_raw_data_profiler adir.name],'r'); %frewind(fid); fgetl(fid); jd=[];hr=[];mn=[];sec=[];snr=[];width=[]; k=0; % start with no scans, increment k for desired (vertical) scans while ~feof(fid) % loop over scans k=k+1; %for pp = 1:20 %short test loop DEW temp1=fgetl(fid); data=sscanf(temp1,'%d'); %bb(k) = data(1); % beam index--arbitrary number el(k)=data(2); % elevation angle, degrees az(k)=data(3); % azimuth angle, degrees %yr(k)=data(4); % years since 1900 jd(k)=data(5); % append each scan hr(k)=data(6); mn(k)=data(7); sec(k)=data(8); %clear temp1 data; temp2=fgetl(fid); data=sscanf(temp2,'%d'); ipp(k)=data(1); % inter-pulse period (ms) pw(k)=data(2); % pulse width (ns) gate1(k)=data(3); % in m gate_space(k)=data(4); last_gate(k)=data(5); %clear temp2 data; for j=1:last_gate(k) % read a line for each gate temp=fgetl(fid); data=sscanf(temp,'%f'); % DEW 5/17/2007 changed from beam test to EL(elevation angle) and PW %if el==90 && pw==400 % Simon eliminated conditions if j==1 %k=k+1; ht(k,1:last_gate(k))=(gate1(k)+(0:1:last_gate(k)-1)*gate_space(k))/1000; % km % store gate heights for each scan end snr(k,j)=data(2); %pow(k,j)=data(3); %end %clear temp data; end %if el==90 && pw==400 % snr_cr=snr+20*log10(ht); %end end fclose(fid); % select only 400 ns pulses, not 700 ns ipulse=(pw(1:k)==400); iz=[true([1 52]) false([1 3])]; snr_cr=NaN+zeros(sum(ipulse),sum(iz)); snr_cr=snr(ipulse,iz)+20*log10(ht(ipulse,iz)); time=hr(ipulse)'*3600+mn(ipulse)'*60+sec(ipulse)'; % find BL top from max snr_cr M=spatialfilter(snr_cr,[1 2 3 2 1],[1 2 1]); % filtered SNR M(isnan(M))=snr_cr(isnan(M)); cb=load(cbfn); %Load Ceilometer cloud base txt file cb_time=cb(:,5)+((cb(:,6)+(cb(:,7)/60.))/60); [junk,ind]=unique(cb_time); iis=isfinite(cb(:,2)) & sparse(ind,1,ones(length(ind),1),length(cb),1); cbase=interp1(cb_time(iis)*3600,cb(iis,2)/1e3,time); cbase(isnan(cbase)|cbase>1.500)=0.500; % default level to look above=500m % look above 800m for inversion height, except for when % cloud base is low, look above ceilometer cloud base h=ht(ipulse,iz); ii=(1:size(h,2))'; ihmax=NaN+zeros(size(snr_cr,1),1); ihmaxtest=ihmax; bltop=ihmax; bltoptest=ihmax; for itime=1:size(snr_cr,1) % search for BL top echo above sea/ship clutter if cbase(itime)<1.000 % cloud base is lower than 1 km. mask=h(itime,:)>cbase(itime) & h(itime,:)<1.500; isrch=find(h(itime,:)>cbase(itime),1,'first'); else mask=h(itime,:)>0.800 & h(itime,:)<2.500; isrch=find((diff(M(itime,:)).*mask(1:end-1))>0,1,'first'); % look above min (filtered) SNR end minin(itime)=h(itime,isrch); mask1=ii<=isrch; % exclude below SNR min. to exclude clutter [junk,ihmaxtest(itime)]=max(M(itime,:)-1000*(mask1'|~mask)); %filtered top [snrmax(itime),ihmax(itime)]=max(snr_cr(itime,:)-1000*(mask1'|~mask)); bltoptest(itime)=h(itime,ihmaxtest(itime)); bltop(itime)=h(itime,ihmax(itime)); end bltop(abs(bltop-bltoptest)>.080)=NaN; % toss when filtered top disagrees by more than one gate % Form the mean every 10 minutes of the 16-84th percentile. minbase=(0:10:24*60)'; % marks beginning of time bins (last bin is empty of data) secbase=minbase*60; bltop10=NaN+zeros(length(minbase)-1,1); for tbin=1:length(minbase)-1 iin=time>=secbase(tbin) & time=4; % 4 needed to toss away outliers ind=0.5:1:nobs(tbin); irange=ind>=0.16*nobs(tbin) & ind<=0.84*nobs(tbin); topsort=sort(bltop(iin)); bltop10(tbin)=mean(topsort(irange)); else bltop10(tbin)=NaN; end end % save 10-min BL top data daily ftop=fopen([way_proc_data_profiler cruise year '_top_' sprintf('%3d',yday) '.txt'],'w'); fprintf(ftop,'%3d %5d %6.04f\n',[repmat(yday,[length(bltop10),1]) secbase(1:end-1,1) bltop10]'); fclose(ftop); % plot to see how BL height looks %figure() clf set(gcf,'defaultaxesfonts',14) set(gca,'fontname','time','fontsize',14); [c,hcf]=contourf((repmat(time,[1 52]))/3600,h,snr_cr,16); set(hcf,'edgecolor','none'); axis xy; hcb=colorbar; set(get(hcb,'ylabel'),'string','rcSNR','fontsize',14); xlim([0 24]) ylim([0 2.5]) %ylim([0 4.29]) set(gca,'xtick',0:4:24) ylabel('Height (km)') xlabel('Hour (UTC)') tit = {[cruise ' ' num2str(year) ' day ' num2str(jd(k))] 'Wind profiler rcSNR'}; title(tit) hold on plot(secbase(1:end-1)/3600,bltop10,'m.') % print_buffer=[way_images_profiler, cruise,'_',year,'_',fn(7:9)]; % print('-djpeg90', print_buffer) % print('-depsc2',print_buffer) plot(cb_time,cb(:,2)/1000,'r.') axis([0,24,0,2.5]) tit(2) = {'Wind profiler rcSNR, cloud top (magenta) and ceilometer cloud base (red)'}; title(tit) % print_buffer=[way_images_profiler,'wind_profiler_cloud_', cruise,year,'_',fn(7:9)]; % print('-djpeg', print_buffer) % print('-depsc2', print_buffer) end