% 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
%##########################################################################
%%%%%%%%%%%%%%%%%%%%%%%%% Cruise informations and parameters %%%%%%%%%%%%%%
%##########################################################################

cruise='AtlanticTransect_2014';             % Cruise name  GC_2014
% year='2011';               % year (deprecated--make year into a double) 
% yyyy='2011';               % year (string)
ship='PolarStern';               % Research vessel  GS

read_parameters_AtlanticTransect_2014;     %reads the parameters for the cruise.

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=1;  %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 %%%%%%%%%%%%%%%%%%%%%%%%%%%
%##########################################################################

% use known day/hour; prompt for it if unknown
if ~exist('ddd','var')
    ddds=input('Input yearday  ','s');
    ddd=str2double(ddds);
else
    ddds=sprintf('%03i',ddd);
end
[m,d]=yd2md(str2double(yyyy), ddd);
Vdate=[str2double(yyyy) m d];
if ~exist('shr','var')
    shrs=input('Input start hr  ','s');
    shr=str2double(shrs);
else
    shr=0;
end
if ~exist('ehr','var')
    ehrs=input('Input end hr  ','s');
    ehr=str2double(ehrs);
else % use hours 0-23 by default
    ehr=23;
end

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_AtlanticTransect_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(1,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)')

% PLOTS FOR RADIOMETER OUT
%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
