% Reads in W-Band nc moment data and plot
% range gates are 25 m
% timing is every 0.25 sec
% Version_3 and above read kongsberg motion data and corrects vertical velocity
%The read_parameters.m file sets the paths for the radar data
% 96  W-band radar files
% 97  way_raw_data_wband=[big_dir '\RHB\radar\wband\Raw\'];
% 98  way_raw_images_wband=[big_dir '\RHB\radar\wband\Raw_images\'];
% You must go to lines 97-98 and set the path for your data archive
%close all
clear

%%%%%%%%%%%%%%%%%%%Modifies current directory and search path%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fluxroot='D:\Data\DYNAMO_2011\Revelle\Scientific_analysis\programs\';  % 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. 

%%%%%%%%%%%  Read parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     read_parameters_DYNAMO_2011;     %reads the parameters for the cruise.
%     currentMname='DYNAMO_2011';
currentMname=mfilename;   %Gets the name of currently running M-file to pick 
currentMname=currentMname(17:end);
%     eval([ 'read_parameters_', currentMname(11:end)])

%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, '.'));


%################## Defines some options and parameters ##################################
Nrange=220;
RangeRes=0.031;
ht =single((1:1:Nrange)*RangeRes); % range gates in km
zmxx=single(ht(Nrange));
dbz0=-43;
yyyy='2011';
year='2011';
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=0;  %to save individual hourly graphs

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=1;  %to add mailbox data to the Wband plots (note you need to run the mailbox codes first)


%%%%%%%%%%%%%%%%%%%%%%%%% Starts process %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% t = help('netcdf');  %checks for Matlab capabilities
% if isempty(t)
    Matlab_netcdf=0;
% else
    Matlab_netcdf=1;
% end

ddds = input('Input yearday  ','s');
ddd = str2num(ddds);
[m,d]=yd2md(str2num(yyyy), ddd);
Vdate=[str2num(yyyy) m d];
shrs = input('Input start hr  ','s');
shr = str2num(shrs);
ehrs = input('Input end hr  ','s');
ehr = str2num(ehrs);

dr=dir([way_raw_data_wband yyyy ddds  '*.nc']);
film=dr.name;
% jamn=1; %useless
cldtop=[];DVx=[];Refx=[];hhx=[];SNRx=[];SWx=[];jdx=[];
for jam = shr:ehr
    %close all
    clear hh y z dt jd
    hr=sprintf('%02i',jam);     
%     film=dr(jamn).name
    fn = [way_raw_data_wband yyyy ddds hr film(10:11) 'MMCRMom.nc'];%[way_raw_data_wband film];
    disp(['Reads 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
%         jamn=jamn+1;
        if Matlab_netcdf
            base_time=netcdf.getVar(FID,0); % bug for multiple files in 1 hour!   0
            time_offset=netcdf.getVar(FID,1)';          
        else
            ncload(fn);   
        end;
        tt = time_offset + double(base_time);
        % calculate time form time_offset and base_time

        %Would be better to use the following but not sure if crashes on
        %DAS...
%         z = tt/(60*60*24)+datenum([1970,1,1,0,0,0]);
%         dt = datestr(z, 31);    %2008-10-07 16:35:53
%         [yr, month, day, hour, minutes, secs] = datevec(dt);
%         secs = mod(time_offset,60)';
%         jd = (datenum(yr,month,day,hour,minutes,secs)-datenum(yr-1,12,31))';
                % calculate time form time_offset and base_time
        for i = 1:length(tt)
            %X = base_time;
            y(i) =tt(i)/(60*60*24);
            z(i) =y(i)+datenum([1970,1,1,0,0,0]);
            dt = datestr(z(i), 31);    %2008-10-07 16:35:53
            yr = str2num(dt(1:4));
            month = str2num(dt(6:7));
            day = str2num(dt(9:10));
            hour = str2num(dt(12:13));
            minutes = str2num(dt(15:16));
            secs = mod(time_offset(i),60);%str2num(dt(18:19));
            jd(i) = datenum(yr,month,day,hour,minutes,secs)-datenum(yr-1,12,31);
            %correct for
        end
        hh = (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;
        RangeRes=ht(2)-ht(1);dbz0=-43;if RangeRes>0.04;dbz0=-51;end;
        
        DV=single([]);
        if motionread%correct doppler for
               read_Wband_motion_lb
               try
                DV_mot=single(interp1(hhk,vkons(9,:),hh));
               catch
                DV_mot=0*DVa(1,:)   ;
                end;
               ii=find(isnan(DV_mot));DV_mot(ii)=0;
               for im=1:length(ht);
                    DV(im,:)=DVa(im,:)+DV_mot;
               end;
        else
            DV=DVa;
        end;
        DV=single(DV);
        clear Reflectivity SpectralWidth MeanDopplerVelocity SignalToNoiseRatio Power NoiseLevel SkyNoiseLevel DVa DV_mot vkons CircularDepolarizationRatio
        clear RangeCorrectedPower AvgNoiseLevel Azimuth
        for k=1:length(Ref);
            kk=find(Ref(:,k)-dbz0-20*log10(ht')>8 & ht'>.1 & ht'<zmxx & Ref(:,k)>-30);
            if ~isempty(kk)
                htt(k)=single(ht(kk(length(kk))));
                reft(k)=single(mean(10.^(Ref(kk,k)/10)));
                reft2(k)=single(mean(10.^((Ref(kk,k)-dbz0-20*log10(ht(kk)'))/10)));
                wvel(k)=single(mean(DV(kk,k)));
                nnt(k)=single(length(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=find(htt>=y(1) & htt<=y(5));
        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=single(sum(nnt));
        sigh=single(sqrt(nanmean(htt(kk).^2)-httb.^2));
        ['Hour= ' num2str((hh(1))) '  Mean cloud top ht= ' num2str(httb) '  std cloud top height = ' num2str(sigh)]
        ['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)];
        %DVx=[single(DVx) single(DV)];
        Refx=[single(Refx) single(Ref)];
%         SNRx=([single(SNRx) single(SNR)]);
        %SWx=[single(SWx) single(SW)];
        %jdx=[single(jdx) (jd)];
        ymaxx=7.0;if RangeRes>0.04;ymaxx=ht(Nrange);end;
        
        figure(1)
        clf
        subplot(3,1,1),imagesc(hh,ht, Ref(1:Nrange,:),[-60, 20]);  %
        axis xy;ylim([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',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none')
        h = colorbar;
        axes(h)
        ylabel('Ref dBZ')

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

        subplot(3,1,3),imagesc(hh,ht, SW(1:Nrange,:),[0, 5]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        xlabel('Hrs (UTC)')
        h = colorbar;
        axes(h)
        ylabel('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'); 
        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]);xlim([ddd ddd+1]); end
               
               
        figure(2)
        clf
        subplot(3,1,1),imagesc(hh,ht, SNR(1:Nrange,:),[-25, 20]);  %
        axis xy;ylim([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',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,hr),'FontWeight','Bold','Interpreter','none')
        h = colorbar;
        axes(h)
        ylabel('SNR dBZ')

        subplot(3,1,2),imagesc(hh,ht, DV(1:Nrange,:),[-3, 5]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        %axis([15 15.25 0 3])
        h = colorbar;
        axes(h)
        ylabel('DopVel (m/s)')
%         
%         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;
%         axes(h)
%         ylabel('DopVel (m/s)')

        subplot(3,1,3),imagesc(hh,ht, SW(1:Nrange,:),[0, 5]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        xlabel('Hrs (UTC)')
        h = colorbar;
        axes(h)
        ylabel('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'); 
        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]);xlim([ddd ddd+1]); end

        %fclose(FID);

    end % test for file

end % hr loop

%       
%clear Ref SW SNR DV
%%%%%%%%%%%%%%  Some diagnostics CWF added on the vocals cruise
% ceilo=1;
if ceilo
%*******************  Load 1 hr ceilo file
    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
%*******************  Load 1 hr ceilo file
    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 15 s ceilo file              
   

%     micro=1;
if micro   
        %CALNEX_MB_sci_day
        %fn2 = [way_proc_data_micro 'vocals_rhb_lwp10min_v1.nc'];
        mby=single([]);mbyy=single([]);mb=single([]);
%         fn2 = [way_proc_data_mb 'DYNAMO_2011_micro_s.txt'];
        %reads daily file instead of concatenated file
        fn2=[way_proc_data_mb, 'DYNAMO_2011_mwr_' yyyy(3:4) sprintf('%02i%02i',Vdate(2),Vdate(3))  '.txt'];
        mby=single(load(fn2));
        %mby= DYNAMO_2011_micro_s;
        ii=find(~isnan(mby(:,1)));mb=mby(ii,:);
%         fn2 = [way_proc_data_micro 'vocals_rhb_lwp_v1.nc'];
%        FID2 = netcdf.open(fn2,'NC_NOWRITE');
%         mby(:,1)=netcdf.getVar(FID2,0);
%         mby(:,2)=netcdf.getVar(FID2,1);  
%         mby(:,3)=netcdf.getVar(FID2,2);
%         mby(:,4)=netcdf.getVar(FID2,3);  
%         mby(:,5)=netcdf.getVar(FID2,4);
%         ii=find(mby(:,3))<1e30;mbyy=mby(ii,:);
%         jk=find(diff(mbyy(:,1))<-0.01);mb=mbyy(1:jk(1)-1,:);jk=find(diff(mb(:,1))<0.0);mb(jk,1)=mb(jk,1)+1e-4;
        
        [b, m, n] = unique(mb(:,1));
        jdm=mb(m,1)+0/2/3600/24;%correct time to middle of 1 min window
         ii=find(mb(:,5)>1e3);mb(ii,5)=NaN; ii=find(mb(:,5)<-50);mb(ii,5)=NaN;

%         mbx=mb;%interp1(jdm,mb(m,:),jdx);
        jdmx=mb(:,1);
        vap=single(mb(:,2));ii=find(vap<0.1);vap(ii)=NaN;%path-integrated water vapor, cm (8.3f)
        liq=single(mb(:,3)*1e3);%path-integrated liquid water, Grams/Meter2 
%         vap=mbx(:,4);%path-integrated water vapor, cm (8.3f)
%         liq=mbx(:,5);%path-integrated liquid water, Grams/Meter2 
%         ii=find(mb(:,5)>1e3);mb(ii,5)=NaN; ii=find(mb(:,5)<-50);mb(ii,5)=NaN;
%         TB1=mbx(:,13);%brightness temperature from Lowest radiometer channel, deg. K 
%         TB2=mbx(:,14);%brightness temperature from next radiometer channel, deg. K 
        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=1:2:length(hhx);
        figure
        %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]);
        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');
        axes(h)
        ylabel('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');
%         axes(h)
%         ylabel('DopVel (m/s)')
        
        subplot(2,1,2);plot((jdmx-ddd)*24,liq);ylim([-20 300]);axis([shr ehr+1 -20 300]);
        xlabel('Year Day (UTC)');ylabel('LWP (g/m^2)');
        %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]);xlim([ddd ddd+1]); end

      

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;
    ii=find(cldtop(:,8)<100);cldtop(ii,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-');
  %legend('Wbm','WbB','WbT' ,'location','best')

  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('Cm','C1','C2' ,'location','best')
  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],'--' );
  Ylabel('dBZ and 10*W (m/s)');
  xlabel(['Day-Of-Year (' year ')']); 
%   print('-djpeg90 ',[way_cloud_data_wband '\cloud_stats_' num2str(ddd), '.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_cloud_data_wband cruise '_cloud_stats_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]);xlim([ddd ddd+1]); end

  if ceilo
    ii=find(floor(jdh)==ddd);
    hgtsi=interp1(jdh(ii),cloudhgt(ii),jdt,'nearest');
    figure;plot(jdt,cldtop(:,2)-hgtsi/1000);axis([ddd ddd+1 0 .5]);
    Ylabel('Cloud Thickness (m)');
    xlabel(['Day-Of-Year (' year ')']); 
%     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_cloud_data_wband cruise '_cloud_thick_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]);xlim([ddd ddd+1]); end

  else
      hgtsi=NaN*jdt;
  end;
  cdgf=[cldtop hgtsi/1e3];
  f=[[way_cloud_data_wband '\Radar_Cldtop_' ddds '.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;
 detect=zeros(size(Ref));
 for i=1:length(Ref)
      kk=find(Ref(:,i)-dbz0-20*log10(ht')>8 & ht'>.1 & ht'<zmxx );
      detect(kk,i)=1;
 end;
 xRef=Ref.*detect;
 xDV=DV.*detect;
 xSW=SW.*detect;
 xZZ=10.^(Ref/10).*detect;

 dk=length(Ref)-1;kk=1;
 xRefm=mean(xRef(:,kk:kk+dk)');
 xZZm=mean(xZZ(:,kk:kk+dk)');
 xDVm=mean(xDV(:,kk:kk+dk)');
 xSWm=mean(xSW(:,kk:kk+dk)');
 detectm=mean(detect(:,kk:kk+dk)');
 mm=find(detectm<5/dk);xRevm(mm)=NaN;xDVm(mm)=NaN;xZZm(mm)=NaN;
 dZBbar=10*log10(xZZm)./detectm;mm=find(dZBbar<-60);dZBbar(mm)=NaN;
 
 ij=find(ht>2);
 wmm=mean(mean(DV(ij,kk:kk+dk)'));wdm=mean(mean(SW(ij,kk:kk+dk)'));
 figure;
 subplot(1,2,1);plot(xRefm./detectm,ht,mean(Ref(:,kk:kk+dk)'),ht,dZBbar,ht,'--',max(-60,dbz0+20*log10(ht')),ht,[-17 -17],[0 3],'--',[0 0],[0 3]);Ylabel('Height (m)');xlabel('dBZ');%axis([-60 20 0 2]);
 subplot(1,2,2);plot(detectm,ht);xlabel('Return Fraction');%axis([0 1 0 2]);
    
 DVbar=xDVm./detectm-wmm;
 dBZ=xRefm./detectm;
 Widthbar=xSWm./detectm;
 wdbm=nanmean(Widthbar(ij));
 
 DVbar2=mean(DV(:,kk:kk+dk)')-wmm;
 dBZ2=mean(Ref(:,kk:kk+dk)');
 ZZ=mean(10.^((Ref(:,kk:kk+dk)')/10));
 dBZbar2=10*log10(ZZ);
 Widthbar2=mean(SW(:,kk:kk+dk)');
 
 figure;
 subplot(1,2,1);
 plot(xDVm./detectm-wmm,ht,mean(DV(:,kk:kk+dk)')-wmm,ht,'--',[0 0],[0 3]);Ylabel('Height (m)');xlabel('Mean Velocity (m/s)');%axis([-1 5 0 2]);
 subplot(1,2,2);plot(xSWm./detectm,ht,mean(SW(:,kk:kk+dk)'),ht,'--',[0 0],[0 3]);xlabel('Doppler Width (m/s)');%axis([0 1 0 2]);
 
 %From Frisch et al 1995
 A=1.2e-4;
 B=1e-5;
 
 %sigx=sqrt(Widthbar.^2-wdbm^2)./DVbar;
 sigx=sqrt(log(1+(Widthbar.^2-wdbm^2)./(DVbar+B/A).^2));
  ii=find(sigx>1);sigx(ii)=NaN;DVbar(ii)=NaN;
 yz=dBZ-132+26*sigx.^2;
 Ql=(A*DVbar).^(-3).*10.^(yz/10);ii=find(Ql>10);Ql(ii)=NaN;%g/m^3
 Fq=Ql.*((DVbar+B/A).*exp(-3*sigx.^2)-B/A)*1e-6*1000*3600;%mm/hr
 r0=1e3*(A*DVbar+B).*exp(-6.5*sigx.^2);%drizzle radius in mm
 
% sigx2=sqrt(Widthbar2.^2-wdm^2)./DVbar2;
 sigx2=sqrt(log(1+(Widthbar2.^2-wdm^2)./(DVbar2+B/A).^2));
 ii=find(sigx2.^2>1 | sigx2.^2<0);sigx2(ii)=NaN;DVbar2(ii)=NaN;dBZ2(ii)=NaN;
 yz2=dBZ2-132+26*sigx2.^2;yz2=min(yz2,-225);
 Ql2=(A*DVbar2).^(-3).*10.^(yz2/10);ii=find(Ql2>10);Ql2(ii)=NaN;%g/m^3
 Fq2=Ql2.*((DVbar2+B/A).*exp(-3*sigx2.^2)-B/A)*1e-6*1000*3600;%mm/hr
 r02=1e3*(A*DVbar2+B).*exp(-6.5*sigx2.^2);%drizzle radius in mm
% ii=find(r02<0);r02(ii)=NaN;
 figure;
 subplot(1,2,1);
 plot(Ql.*detectm,ht);Ylabel('Height (m)');xlabel('Liquid Water (g/m^3)');axis([0 1.5 0 zmxx]);
 subplot(1,2,2);plot(Fq,ht);xlabel('Rainrate (mm/h)');axis([0 .5 0 zmxx]);