%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
% revised for 2001: March 2008
% Adapted from profiler915_BLHT.m by Dan Wolfe
% 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=292:294; %290-308 loop year day of files to read
    fn = 'rbc107xxx.r*w'; % read .raw or .rew
    fn(7:9) = sprintf('%3d',yday);
    adir=dir([way_raw_data_profiler fn]);

    cbfn = way_proc_data_ceilo;
    datev=datevec(datenum(2007,0,0,0,0,0)+yday);
    cbfn1=sprintf('%03d_%02d_%02d_%02d.txt',yday,datev(2:3),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');
    fgetl(fid);

    jd=[];hr=[];mn=[];sec=[];snr=[];width=[];
    k=1;

    while ~feof(fid)
        temp1=fgetl(fid);
        data=sscanf(temp1,'%d');
        bb = data(1);
        el = data(2);
        az = data(3);
        jd(k)=data(5);
        hr(k)=data(6);
        mn(k)=data(7);
        sec(k)=data(8);
        %clear temp1 data;

        temp2=fgetl(fid);
        data=sscanf(temp2,'%d');
        pw = data(2);

        gate1=data(3); % in m
        gate_space=data(4);
        last_gate=data(5);
        % there are only vertical beams
        ht=(gate1+(0:1:last_gate-1)*gate_space)/1000; % in km
        %clear temp2 data;
        if el == 90 && pw == 400 % elevation angle, a
            for j=1:last_gate
                temp=fgetl(fid);
                data=sscanf(temp,'%f');
                snr(k,j)=data(2); % moved before incrementing k
                width(k,j)=data(4);
                %disp(['k=',num2str(k)])
            end
            k=k+1;
        end
    end
    fclose(fid);

    time=hr'*3600+mn'*60+sec';
    %k is the number of observations

    % Range-corrected SNR
    snr_cr=NaN+snr;
    for j=1:last_gate
        snr_cr(:,j)=snr(:,j)+20*log10(ht(j));
        %pow_cr(:,j)=pow(:,j)+20*log10(height(j));
    end

    % 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;
    
    % look above 800m for inversion height, except for days<=294, when
    % cloud base is low, look above ceilometer cloud base
    ii=(1:length(ht))';
    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 yday<=294
            mask=ht>cbase(itime) & ht<1.500;
            isrch=find(ht>cbase(itime),1,'first');
        else
            mask=ht>0.800 & ht<2.500;
            isrch=find((diff(M(itime,:)).*mask(1:end-1))>0,1,'first'); % look above min (filtered) SNR
        end
        minin(itime)=ht(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)=ht(ihmaxtest(itime));
        bltop(itime)=ht(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<secbase(tbin+1) & isfinite(bltop);
        nobs(tbin)=sum(iin);
        if nobs(tbin)>=4; % 4 needed to toss away outliers
            ind=0.5:1:nobs;
            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);
    imagesc((time-time(1))/3600,ht,snr_cr',[-10 20])
    axis xy;
    hcb=colorbar; set(get(hcb,'ylabel'),'string','rcSNR','fontsize',14);
    xlim([0 23])
    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 ' fn(7:9)] '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