% 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 all

%%%%%%%%%%%%%%%%%%%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, 'flux'])         %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 ##################################
motionread=1;%set to zero to skip motion correction
ht = (1:1:120)*.025; % range gates in km
zmxx=ht(120);
dbz0=-42;
yyyy='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




%%%%%%%%%%%%%%%%%%%%%%%%% Starts process %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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=0;
cldtop=[];DVx=[];Refx=[];hhx=[];SNRx=[];SWx=[];jdx=[];
for jam = shr:ehr
    %close all
    clear hh y z dt jd
    hr=sprintf('%02i',jam);
    
%     varname_all=[];
%     for ii=1:43
%     [varname,xtype,dimids,natts] = netcdf.inqVar(FID,ii);
%     varname_all=[varname_all, ' ', varname];
%     end
    
    jamn=jamn+1;
    film=dr(jamn).name
    fn = [way_raw_data_wband yyyy ddds hr film(10:11) 'MMCRMom.nc'];
    FID = netcdf.open(fn,'NC_NOWRITE');%fopen(fn,'r');
    if FID > 0  % Test to see if file exists
       % ncload(fn)
         base_time=netcdf.getVar(FID,0); % bug for multiple files in 1 hour!   0
        % read 1D vectors
        time_offset=netcdf.getVar(FID,1)';
        tt = time_offset + double(base_time);
        % 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
        z=netcdf.getVar(FID,27);% 27
        Ref=netcdf.getVar(FID,40);  %40
        SNR=netcdf.getVar(FID,36); %36
        SW=netcdf.getVar(FID,38);  % 38
        DVa=netcdf.getVar(FID,35); %35
        netcdf.close(FID);
%         Ref = Reflectivity';
%         SW = SpectralWidth';
%         DVa = MeanDopplerVelocity';
%         SNR = SignalToNoiseRatio';
        DV=[];
        if motionread%correct doppler for
               read_Wband_motion_lb
               DV_mot=interp1(hhk,vkons(9,:),hh);
               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;
        
        clear Reflectivity SpectralWidth MeanDopplerVelocity SignalToNoiseRatio
        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)=ht(kk(length(kk)));
                reft(k)=mean(10.^(Ref(kk,k)/10));
                reft2(k)=mean(10.^((Ref(kk,k)-dbz0-20*log10(ht(kk)'))/10));
                wvel(k)=mean(DV(kk,k));
                nnt(k)=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=nanmean(htt(kk));
        reftb=10*log10(nanmean(reft(kk)));
        reft2b=10*log10(nanmean(reft2(kk)));
        wvelb=nanmean(wvel(kk));
        nntb=sum(nnt);
        sigh=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=[hhx hh];
        DVx=[DVx DV];
        Refx=[Refx Ref];
        SNRx=[SNRx SNR];
        SWx=[SWx SW];
        jdx=[jdx jd];
        ymaxx=2.0;
        
        figure(1)
        clf
        subplot(3,1,1),imagesc(hh,ht, Ref(1:120,:),[-60, 20]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Height (km)')
        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:120,:),[-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:120,:),[0, 5]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        xlabel('Hrs (UTC)')
        h = colorbar;
        axes(h)
        ylabel('SpecWid (m/s)')
        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_' 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:120,:),[-25, 20]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        %axis([15 15.25 0 3])
        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:120,:),[-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:120,:),[0, 5]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        xlabel('Hrs (UTC)')
        h = colorbar;
        axes(h)
        ylabel('SpecWid (m/s)')
        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_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

%