% Wband deployment on German ship - 2014
% Reads in W-Band nc moment data and plot
% range gates could be variable but currently set to 25m
% timing is every 0.25 sec
% Read kongsberg motion data and corrects vertical velocity
% read_parameters.m sets the paths for the radar data
% Use jamn for multiple files in 1 hour.
% Preallocate arrays to combine full-day time series.
% 01/08/2014 -LB

close all
clear all

try
    %##########################################################################
    %%%%%%%%%%%%%%%%%%%%%%%%% Cruise informations and parameters %%%%%%%%%%%%%%
    %##########################################################################

    Fluxroot='E:\AtlanticTransect_2014\PolarStern\Scientific_analysis\programs\Wband\';  % name of the path where the 'flux' folder + read_parameters file are installed.

    cd([Fluxroot])         %Changes current directory to fluxroot/flux
    addpath(genpath(Fluxroot))     %adds flux folder and all of its subdirectories to the top of the MATLAB search path.
    
    %Gets the last uploaded parameter file from Fluxroot directory
    dirfluxroot=dir([Fluxroot 'read_parameters_*.m']);
    MostrecentPARfile=dirfluxroot(cellfun(@(x) x==max([dirfluxroot.datenum]),{dirfluxroot(:).datenum})).name;
    disp(['Most recent parameters being used: ' MostrecentPARfile])
    eval(strtok(MostrecentPARfile, '.'));

    %%%%%%%%%%%%%%%%%%%Modifies current directory and search path%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Ceiloroot='D:\Data\DYNAMO_2011\Revelle\Scientific_analysis\programs\';  % name of the path where the 'flux' folder + read_parameters file are installed.
    cd([Fluxroot, 'Wband'])         %Changes current directory to fluxroot
%     addpath(genpath(Ceiloroot))     %adds flux folder and all of its subdirectories to the top of the MATLAB search path.
    
    graphformat='.png';  %select graphics format files
    graphdevice='-dpng'; %select graphic device
    prtit=true;  % for printing plots--ineffective unless plotit is also true going to a temp directory
    prtit_hr=1;  %to save individual hourly graphs
    modeflag=0;
    motionread=1;%set to zero to skip motion correction
    ceilo=0;  %to add ceilo data to the Wband plots (note you need to run the ceilo codes first)
    micro=0;  %to add mailbox data to the Wband plots (note you need to run the mailbox codes first)
    Matlab_netcdf=1;

    %##########################################################################
    %%%%  Starts reading and ourly plotting process %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %##########################################################################

    dv = clock;  % gets current date and time as date vector.
    Vdate=datevec(datenum(dv(1),dv(2),dv(3)-1));
    ddd=datenum(Vdate)-datenum(Vdate(1),1,0);  %same as md2yd(V(1),V(2),V(3))   
%     yyyy=sprintf('%04i',Vdate(1));
    ddds=sprintf('%03i',ddd);
    shr=0;shrs=sprintf('%02i',shr);
    ehr=23;ehrss=sprintf('%02i',ehr);

    cldtop=[];
    % preallocate x-extended variables
    nxmax=175000;
    Nrange=220; % may only use 190 of the ranges or less or more
    Refx=deal(single(zeros(Nrange,nxmax))+NaN);
    %[SNRx, SWx, DVx]=deal(Refx); % not used
    [hhx, rangemode]=deal(single(zeros(1,nxmax))+NaN);

    Refx=deal(single(zeros(Nrange,nxmax))+NaN);
    [Refh,SNRh, SWh, DVah,]=deal(single(zeros(Nrange,15000))+NaN);
    [hh]=deal(single(zeros(1,15000))+NaN);


    handl(1)=figure(1); % open 2 figures whether anything to plot
    handl(2)=figure(2);
    % set(handl(:),'WindowStyle','docked');

    nx=0; % count position in extended variables
    for jam=shr:ehr
        isfile=false;
        clf(handl(:));
        nh=0;
        clear y z dt jd
        hh=[]; Ref=[]; DVa=[]; SW=[]; SNR=[]; nrg=[];
        hr=sprintf('%02i',jam);
        dr=dir([way_raw_data_wband yyyy ddds hr '*MMCRMom.nc']);
        if isempty(dr)
            error([way_raw_data_wband yyyy ddds hr '*MMCRMom.nc' '<---file not found'])
        end
        for jamn=1:length(dr) % multiple files in one hour
            fn=[way_raw_data_wband dr(jamn).name];
            disp(['Reading file: ' fn])
            try
                if Matlab_netcdf
                    FID=netcdf.open(fn,'NC_NOWRITE');
                else
                    FID=fopen(fn,'r');
                end
            catch
                FID=-1;
            end

            if FID>0  % Test to see if file exists
                isfile=true;
                if Matlab_netcdf
                    base_time=netcdf.getVar(FID,0);
                    time_offset=netcdf.getVar(FID,1)';
                else
                    ncload(fn);
                end
                % calculate time from time_offset and base_time
                tt=time_offset+double(base_time);
                nta=length(tt); % time length of current file
                z=tt/(60*60*24)+datenum([1970,1,1,0,0,0]);
                jd=z-datenum(str2double(yyyy)-1,12,31);
                %         dt = datestr(z, 31);    %2008-10-07 16:35:53
                %         [yr, month, day, hour, minutes, secs] = datevec(dt);

                hh0=(jd-(fix(jd(1))))*24; % subtract off JD and convert to an hour

                if Matlab_netcdf
                    z=netcdf.getVar(FID,27);% 27
                    ht=z(:,1)/1e3;
                    Reflectivity=netcdf.getVar(FID,40)';  %40
                    SignalToNoiseRatio=netcdf.getVar(FID,36)'; %36
                    SpectralWidth=netcdf.getVar(FID,38)';  % 38
                    MeanDopplerVelocity=netcdf.getVar(FID,35)'; %35

                    Ref=single(Reflectivity)';
                    SW=single(SpectralWidth)';
                    DVa=single(MeanDopplerVelocity)';
                    SNR=single(SignalToNoiseRatio)';
                    nrg=netcdf.getVar(FID,18);
                        Nrange=nrg(1);
                    netcdf.close(FID);              
                else
                    z=single(Heights);
                    ht=z(1,:)/1e3;
                    Nrange=single(NumHeights(1));
                    Ref=single(Reflectivity)';
                    SW=single(SpectralWidth)';
                    DVa=single(MeanDopplerVelocity)';
                    SNR=single(SignalToNoiseRatio)';
                    fclose(FID);
                end
            end % test for file

            %concatenate multiples files into one hour    
            [m0,n0]=size(Ref);
            hh(nh+(1:nta))=hh0;
    %         hth(1:m0)=ht;
            Refh(1:m0,nh+(1:nta))=Ref;
            DVah(1:m0,nh+(1:nta))=DVa;
            SWh(1:m0,nh+(1:nta))=SW;
            SNRh(1:m0,nh+(1:nta))=SNR;
            nh=nh+nta;    

        end % one or multiple files in hour

        % truncate hourly arrays
        clear hh0 Ref DVa SW SNR 
        Ref=Refh(1:Nrange,1:nh);
        DVa=DVah(1:Nrange,1:nh);
        SW=SWh(1:Nrange,1:nh);   
        SNR=SNRh(1:Nrange,1:nh); 
    %     ht=hth(1:nh);
        clear Refh DVah SWh SNRh hth


        dbz0=-43;
        RangeRes=ht(2)-ht(1);
        if RangeRes>0.04;
            dbz0=-51;
        end

        % store 2 height variables for the 2 modes if necessary
        if length(ht)==220
            ht1=ht;
        elseif length(ht)==190
            ht0=ht;
        end


        if motionread %correct doppler for ship motion
        %                 if exist([way_raw_motion_wband yyyy ddds hr '00Kongsberg.txt'],'file');
           if ~isempty(dir([way_raw_motion_wband year ddds hr '*Kongsberg.txt'])); %no files at all
                read_Wband_motion_GS_2014
                pitchk=vkonsa(7,:);
                rollk=vkonsa(8,:);
                try
                    DV_mot=single(interp1(hhka,vkonsa(9,:),hh));            
                    Konmotion=1;
                catch  % interpolation failed
                    DV_mot=0*DVa(1,:);
                    Konmotion=0;
                end
                DV_mot(isnan(DV_mot))=0;
                DV=DVa+repmat(DV_mot,Nrange,1);
            else
                DV=DVa;% no motion file
                Konmotion=0;
                [hhka,pitchk, rollk]=deal(single(zeros(1,36000))+NaN);
            end
        else
            DV=DVa; %motion correction turned off
            Konmotion=0;
        end
        DV=single(DV);



    %--------------------------------------------------------------------------
    %%%%%    do hourly plots
    %--------------------------------------------------------------------------

        ymaxx=ht(Nrange);
        zmxx=ymaxx;

        figure(1); clf
        subplot(3,1,1),
        imagesc(hh,ht, Ref(1:Nrange,:),[-60, 20]);
        axis xy; axis([jam jam+1 0 ymaxx]);
        ylabel('Ht (km)')
        %         tit={[cruise ' W-Band'' YD ' ddds ' Hr ' hr]};
        %         title(tit,'Interpreter','none')
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i, Hr-%s).  W-Band (motionread=%i, Kongsberg=%i)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr,motionread,Konmotion),'FontWeight','Bold','Interpreter','none')
        h=colorbar;
        ylabel(h,'Ref dBZ')
        hold on

        subplot(3,1,2)
        imagesc(hh,ht, DV(1:Nrange,:),[-7, 7]);  %
        axis xy;axis([jam jam+1 0 ymaxx]);
        ylabel('Ht (km)')
        h=colorbar;
        ylabel(h,'DopVel (m/s)')
        hold on

        subplot(3,1,3),imagesc(hh,ht, SW(1:Nrange,:),[0, 5]);  %
        axis xy;axis([jam jam+1 0 ymaxx]);
        ylabel('Ht (km)')
        xlabel('Hrs (UTC)')
        h=colorbar;
        ylabel(h,'SpecWid (m/s)')
        % print('-djpeg90 ',[way_raw_images_wband '\Wband_ref_' num2str(ddd), '_' hr '.jpg']);
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
        hold on
        if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_ref_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]);end

        figure(2); clf
        subplot(3,1,1),imagesc(hh,ht, SNR(1:Nrange,:),[-25, 20]);  %
        axis xy; axis([jam jam+1 0 ymaxx]);
        ylabel('Ht (km)')
        %axis([15 15.25 0 3])
        %         tit={[cruise ' W-Band'' YD ' ddds ' Hr ' hr]};
        %         title(tit,'Interpreter','none')
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i, Hr-%s).  W-Band (motionread=%i, Kongsberg=%i)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr,motionread,Konmotion),'FontWeight','Bold','Interpreter','none')
        h=colorbar;
        ylabel(h,'SNR dBZ')
        hold on

        subplot(3,1,2),imagesc(hh,ht, DV(1:Nrange,:),[-3, 5]);  %
        axis xy;axis([jam jam+1 0 ymaxx]);
        ylabel('Ht (km)')
        %axis([15 15.25 0 3])
        h=colorbar;
        ylabel(h,'DopVel (m/s)')
        hold on
        %         subplot(3,1,3),imagesc(hh,ht, DVa(1:120,:),[-3, 5]);  %
        %         axis xy
        %         ylabel('Ht (km)')
        %         %axis([15 15.25 0 3])
        %         h=colorbar;
        %         ylabel(h,'DopVel (m/s)')

        subplot(3,1,3),imagesc(hh,ht, SW(1:Nrange,:),[0, 5]);  %
        axis xy;axis([jam jam+1 0 ymaxx]);
        ylabel('Ht (km)')
        xlabel('Hrs (UTC)')
        h=colorbar;
        ylabel(h,'SpecWid (m/s)')
        %print('-djpeg90 ',[way_raw_images_wband '\Wband_snr_' num2str(ddd), '_' hr '.jpg']);
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
        hold on
        if prtit_hr; print(graphdevice,[way_raw_images_wband cruise '_Wband_snr_' sprintf('%04i_%02i_%02i_%03i_%02i',Vdate(1),Vdate(2),Vdate(3),ddd,jam) graphformat]);                 end

        clear Reflectivity SpectralWidth MeanDopplerVelocity SignalToNoiseRatio Power NoiseLevel SkyNoiseLevel vkons CircularDepolarizationRatio
        clear RangeCorrectedPower AvgNoiseLevel Azimuth


    %--------------------------------------------------------------------------
         %insert data in extended daylong array(s)
    %--------------------------------------------------------------------------
        [m1,n1]=size(Ref);
        if m1>0 && n1>0
            hhx(:,nx+(1:n1))=single(hh);
            Refx(1:m1,nx+(1:n1))=single(Ref);
            rangemode(1,nx+(1:n1))=m1==220; % rangemode=1 for 220 gates, 0 for 190 gates
            nx=nx+nta;
        end
        modeflag=any(rangemode(1:nx)) && any(~rangemode(1:nx)); % true for dissimilar range modes in 1 day

    %--------------------------------------------------------------------------
    % computes some stats
    %--------------------------------------------------------------------------                
        if isfile
            % hourly height statistics
            [htt reft reft2 wvel nnt]=deal(single(zeros(1,size(Ref,2))));
            for k=1:size(Ref,2); % time dimension
                kk=Ref(:,k)-dbz0-20*log10(ht)>8 & ht>.1 & ht<zmxx & Ref(:,k)>-30; % logical for noise&height in bounds
                if sum(kk)
                    htt(k)=single(ht(find(kk,1,'last')));
                    reft(k)=single(nanmean(10.^(Ref(kk,k)/10)));
                    reft2(k)=single(nanmean(10.^((Ref(kk,k)-dbz0-20*log10(ht(kk)))/10)));
                    wvel(k)=single(nanmean(DV(kk,k)));
                    nnt(k)=single(sum(kk));
                else
                    htt(k)=NaN;
                    reft(k)=NaN;
                    reft2(k)=NaN;
                    wvel(k)=NaN;
                    nnt(k)=0;
                end
            end

            y=prctile(htt,[10 15 50 85 90]);
            kk=htt>=y(1) & htt<=y(5); % times when cloud top in bounds
            httb=single(nanmean(htt(kk)));
            reftb=single(10*log10(nanmean(reft(kk))));
            reft2b=single(10*log10(nanmean(reft2(kk))));
            wvelb=single(nanmean(wvel(kk)));
            nntb(1,1:Nrange)=single(sum(nnt));
            sigh=single(sqrt(nanmean(htt(kk).^2)-httb.^2));
            disp(['Hour= ' num2str((hh(1))) '  Mean cloud top ht= ' num2str(httb) '  std cloud top height=' num2str(sigh)])
            disp(['Mean dBZ= ' num2str(reftb) '   Mean dBZ-noisedBZ= ' num2str(reft2b) '   Mean Vel= ' num2str(wvelb)])
            cldtop=[cldtop' [hh(1) httb sigh (y(4)-y(3))/2 reftb reft2b wvelb nntb]']';
            %hhx=[single(hhx) single(hh)];
        end % if there is a file


    end %%%%%%%%%%%  --------->  END of  hr loop


    %##########################################################################
    %%%%%% truncate extended arrays and read ceilo/MWR if selected %%%%%%%%%%%%
    %##########################################################################
    Refx=Refx(:,1:nx);
    hhx=hhx(:,1:nx);
    rangemode=rangemode(:,1:nx);

    %--------------------------------------------------------------------------
    %*******************  Load 1 hr ceilo file
    %--------------------------------------------------------------------------
    if ceilo
        zz=single(load([way_proc_data_ceilo cruise '_ceilo_h_a.txt']));
    else
        zz=single(NaN*ones(1,length(Ref)));
    end
    jdh=(zz(:,1)+.5/24);%1    Day-Of-Year, time at the beginning of 1-hr ave
    allda=single(zz(:,2));%number of obs in 10-min
    nocld=single(zz(:,3));%no clouds
    onecld=single(zz(:,4));%one cloud
    mltcld=single(zz(:,5));%multiple clouds
    obscure=single(zz(:,6));%obscured
    pobs=single(zz(:,7));%partially obscured
    clrf=single(zz(:,8));%clear fraction
    ceilof=single(zz(:,9));%cloudy fraction
    ceilofobs=single(zz(:,10));%cloudy fraction, including obscured
    cloudhgt=single(zz(:,11));%median cloud base height
    hgts1=single(zz(:,12));%15% clouds lower height
    hgts2=single(zz(:,13));%85% clouds lower height
    %--------------------------------------------------------------------------
    %*******************  Load 10 min ceilo file
    %--------------------------------------------------------------------------
    if ceilo
        zz10=single(load([way_proc_data_ceilo cruise  '_ceilo_10min_a.txt']));
    else
        zz10=single(NaN*ones(1,length(Ref)));
    end
    jdh10=(zz10(:,1)+.16666/24);%1    Day-Of-Year, time at the beginning of 10-min ave
    ceilof10=single(zz10(:,9));%cloudy fraction
    cloudhgt10=single(zz10(:,11));%median cloud base height
    hgts110=single(zz10(:,12));%15% clouds lower height
    hgts210=single(zz10(:,13));%85% clouds lower height
    %--------------------------------------------------------------------------
    %*******************  Load MWR file
    %--------------------------------------------------------------------------
    if micro
        mby=single([]);mbyy=single([]);mb=single([]);
        fn2=[way_proc_data_mb cruise '_micro_s.txt'];
        mby=single(load(fn2));
        mb=mby(isfinite(mby(:,1)),:);

        [b, m, n]=unique(mb(:,1));
        jdm=mb(m,1)+0/2/3600/24;%correct time to middle of 1 min window
        mb(mb(:,5)>1e3 | mb(:,5)<-50,5)=NaN;

        jdmx=mb(:,1);
        vap=single(mb(:,2));vap(vap<0.1)=NaN;%path-integrated water vapor, cm (8.3f)
        liq=single(mb(:,3)*1e3);%path-integrated liquid water, Grams/Meter2
        clear mb mbx mby
    else
        jdmx=(ddd+hhx/24);
        liq=NaN*hhx;

    end


    %##########################################################################
    %%%%%%%%%%%%%%%%%% PLOT full time series %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %##########################################################################

    close('all');clear temp DVx;
    % ymaxx=7.0;
    % ii=~mod(0:length(hhx)-1,10); % every tenth
    ii=~mod(0:length(hhx)-1,1); % every points

    %--------------------------------------------------------------------------
    %%%%%%%%%%%  figure Wband_ref_LWP_summary 
    %--------------------------------------------------------------------------
    figure
    if modeflag  %fix graph in case of mode change during doy; plot modes separately on 1 axes
        ll0=~rangemode & ii;  %detects the 190 gates - low res
        ll1=rangemode & ii;  %detects the 220 gates
        subplot(2,1,1)
        imagesc(hhx(ll0),ht0, Refx(1:190,ll0),[-60, 20])
        hold on
        imagesc(hhx(ll1),ht1,Refx(1:220,ll1),[-60, 20]);  %
        axis xy;axis([floor(hhx(1)) ceil(hhx(end)) 0 ht0(end)]);
    else
        %subplot(2,1,1),imagesc(hhx,ht, SNRx(1:Nrange,:),[-25, 20]);  %
        subplot(2,1,1),imagesc(hhx(ii),ht, Refx(1:Nrange,ii),[-60, 20]);  %
        axis xy;ylim([0 ymaxx]);
    end
    ylabel('Ht (km)')
    hold on
    %         plot((zz10(:,1)-ddd)*24,cloudhgt10/1e3,'ok','MarkerSize',4,'MArkerFaceColor','k');
    plot((zz10(:,1)-ddd)*24,cloudhgt10/1e3,'*k','MarkerSize',5);
    hold off
    %axis([15 15.25 0 3])
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  W-Band',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    %         tit=[cruise ' W-Band  YD ' ddds ];
    %         title(tit,'Interpreter','none')
    h=colorbar;
    h2=get(gca,'position');
    ylabel(h,'Refectivity (dBZ)');

    %
    %         subplot(3,1,2),imagesc(hhx,ht, DVx(1:Nrange,:),[-3, 5]);  %
    %         axis xy;ylim([0 ymaxx]);
    %         ylabel('Ht (km)')
    %         %axis([15 15.25 0 3])
    %         h=colorbar;
    %       h2=get(gca,'position');
    %         ylabel(h,'DopVel (m/s)')

    subplot(2,1,2)
    [sjd,ojd]=sort(jdmx);
    plot((sjd-ddd)*24,liq(ojd));ylim([-20 300]);axis([floor(hhx(1)) ceil(hhx(end)) -20 300]);
    xlabel('Year Day (UTC)');ylabel('LWP (g/m^2)');grid
    %title('VOCALS 2008 LWP');
    h3=get(gca,'position'); h3(3)=h2(3); %sets same width for figure
    set(gca,'position',h3)
    %       print('-dpng ',[way_raw_images_wband '\Wband_summary' num2str(ddd), '.png']);
    annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
    if prtit; print(graphdevice,[way_raw_images_wband cruise '_Wband_ref_LWP_summary' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end

    %--------------------------------------------------------------------------
    %%%%%%%%%%%  figures cloud_stats and cloud thickness (if ceilo ON) by CF
    %--------------------------------------------------------------------------
    if length(cldtop)>22;
        %     Upper panel: cloud
        % top height from radar (circle with solid line; ±1 standard deviation, dashed line)
        % and cloud base height  from ceilometer (diamonds with solid line; 1 standard deviation,
        % dots).  Lower panel: mean  radar backscatter for the cloud/drizzle layer
        % (circle with solid line) and mean droplet fall speed (x’s with solid
        % line) multiplied by 10 (in m/s).
        %
        % The txt file has essentially the same info:
        % hour
        % Cloud top ht (km)
        % CLoud top ht std (km)
        % Cloud top ht 1/2 width distribution (km)
        % Mean dBZ
        % Mean dBZ-noise
        % Mean Dopper veloc (+ down)
        % Total number range gates with signif signal (0.3 s, 1hr)
        % Cloud thickness (km)
        %ddd=327;ddds=num2str(ddd);cldtop=load(['cldtop_' ddds '.txt']);

        jdt=ddd+cldtop(:,1)/24+.5/24;
        cldtop(cldtop(:,8)<100,2:7)=NaN;

        figure;
        subplot(2,1,1);plot(jdt,cldtop(:,2),'-o',jdt,cldtop(:,2)+cldtop(:,3),'--',jdt,cldtop(:,2)-cldtop(:,3),'--',jdt,cldtop(:,2)-cldtop(:,8)/12600*25/1e3,'x-');grid
        legend('Wbd Cloud top','Wbd Cloud top +1sigma','Wbd Cloud top -1sigma','','location','NorthWest')

        if ceilo
            hold on;subplot(2,1,1),plot(jdh,cloudhgt/1000,'-d',jdh,hgts1/1000,'.',jdh,hgts2/1000,'.');axis([ddd ddd+1 0 ymaxx]);
            legend('Wbd Cloud top','Wbd Cloud top +1sigma','Wbd Cloud top -1sigma','','Ceilo Cloud base','Ceilo Cloud base 15%','Ceilo Cloud base 85%','location','NorthWest')
        end
        ylabel('Cloud top/base height (km)');
        subplot(2,1,2);plot(jdt,cldtop(:,5),'-o',jdt,10*cldtop(:,7),'-x',[ddd ddd+1],[0 0],'-' ,[ddd ddd+1],[-17 -17],'--' );grid
        legend('Mean dBZ','Mean Dopper veloc (+ down)','location','NorthWest')
        ylabel('dBZ and 10*W (m/s)');
        xlabel(['Day-Of-Year (' yyyy ')']);
        %   print('-djpeg90 ',[way_cloud_data_wband '\cloud_stats_' num2str(ddd), '.jpg']);
        datetick('x','HH:MM','keepticks');
        annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
        if prtit; print(graphdevice,[way_proc_images_wband cruise '_cloud_stats_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end

        if ceilo
            ii=floor(jdh)==ddd;
            if sum(ii)>=2;
                hgtsi=interp1(jdh(ii),cloudhgt(ii),jdt,'nearest');
                if any(isfinite(cldtop(:,2)-hgtsi/1000))
                    figure;
                    plot(jdt,cldtop(:,2)-hgtsi/1000);
                    axis([ddd ddd+1 0 190+30*any(rangemode)]);
                    ylabel('Cloud Thickness (m)');
                    xlabel(['Day-Of-Year (' yyyy ')']); grid
                    datetick('x','HH:MM','keepticks');
                    %     print('-djpeg90 ',[way_cloud_data_wband '\cloud_thick_' num2str(ddd), '_' hr '.jpg']);
                    annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none');
                    if prtit; print(graphdevice,[way_proc_images_wband cruise '_cloud_thick_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]);end
                end
            end
        else
            hgtsi=NaN*jdt;
        end
        cdgf=[cldtop hgtsi/1e3];
        f=[way_proc_data_wband cruise '_Radar_Cldtop_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) '.txt'];
        flist1=fopen(f,'w');
        fprintf(flist1,'%6.1f %9.3f %9.3f %9.3f %9.3f %9.3f %9.3f %9.0f %9.3f  \r\n',cdgf');
        fclose(flist1);
    end

    %##########################################################################
    %%%%%% clear some variables %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %##########################################################################
     clear hh Ref DV SW SNR  vkonsa
     clear hhka pitchk rollk

%%%%%%%%%%%%%%%%%%%%%%%%%% catch error if any, and creates log report %%%%%%%%%%%%%%%%%%%%%%%%%    
catch Wbandstatus1    
    jd=datenum(dv(1),dv(2),dv(3))-datenum(ddv(1)-1,12,31);
    report = getReport(Wbandstatus1);
    fid=fopen([way_proc_data_ceilo 'Wband_logerror' sprintf('%03d',jd)  '.txt'],'w');
    fprintf(fid,'%s',report);
    fclose(fid);
end

exit