% 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
%read_parameters_1;%this program just calls read_parameters_VOCALSz
%cd F:\CALNEX_2010\Atlantis\Scientific_Analysis\programs\Wband\
% read_parameters;
% disp('flocmoi');

%################## Cruise informations  ##################################
cruise='CalNEX';             % Cruise name  
year='2010';               % year (deprecated--make year into a double) 
yyyy='2010';               % year (string)
ship='Atlantis';               % Research vessel
big_dir='C:\Data';

%################## W-band radar files  ###################################
way_raw_data_wband=[big_dir '\CRUISES\DYNAMO_2011\data\wband\'];
way_raw_images_wband=[big_dir '\CRUISES\DYNAMO_2011\data\wband\'];
way_raw_motion_wband=[big_dir '\CRUISES\DYNAMO_2011\data\wband'];
way_cloud_data_wband=[big_dir '\CRUISES\DYNAMO_2011\data\wband\'];

%################## Ceilometer files  #####################################
way_images_ceilo=[big_dir '\CRUISES\DYNAMO_2011\data\ceilometer\']; % daily saved impates
way_raw_data_ceilo=[big_dir '\CRUISES\DYNAMO_2011\data\ceilometer\'];  % Raw Data 
way_proc_data_ceilo=[big_dir '\CRUISES\DYNAMO_2011\data\ceilometer\']; % Processed Data

%################## Radiometer files  #####################################
way_proc_data_micro=[big_dir '\CRUISES\DYNAMO_2011\data\radiometer\'];  % Raw Data 

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


%%%%%%%%%%%%%%%%%%%%%%%%% Starts process %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ddds = input('Input yearday  ','s');
ddd = str2num(ddds);
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(4)
        clf
        subplot(3,1,1),imagesc(hh,ht, Ref(1:120,:),[-60, 20]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
        tit = {['CALNEX' year ' W-Band'' YD ' ddds ' Hr ' hr]};
        title(tit)
        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)')
        print('-djpeg90 ',[way_raw_images_wband '\Wband_ref_' num2str(ddd), '_' hr '.jpg']);

        figure(5)
        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])
        tit = {['CALNEX' year ' W-Band'' YD ' ddds ' Hr ' hr]};
        title(tit)
        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)')
        print('-djpeg90 ',[way_raw_images_wband '\Wband_snr_' num2str(ddd), '_' hr '.jpg']);

        %fclose(FID);

    end % test for file

end % hr loop

%       

%%%%%%%%%%%%%%  Some diagnostics CWF added on the vocals cruise
ceilo=1;
if ceilo
    %*******************  Load 1 hr ceilo file
        zz=load([way_proc_data_ceilo, 'CALNEXceilo_h_a.txt']);
        jdh=zz(:,1)+.5/24;%1    Julian Day, time at the beginning of 1-hr ave
        allda=zz(:,2);%number of obs in 10-min
        nocld=zz(:,3);%no clouds
        onecld=zz(:,4);%one cloud
        mltcld=zz(:,5);%multiple clouds
        obscure=zz(:,6);%obscured
        pobs=zz(:,7);%partially obscured
        clrf=zz(:,8);%clear fraction
        ceilof=zz(:,9);%cloudy fraction
        ceilofobs=zz(:,10);%cloudy fraction, including obscured
        cloudhgt=zz(:,11);%median cloud base height
        hgts1=zz(:,12);%15% clouds lower height
        hgts2=zz(:,13);%85% clouds lower height
        
     %*******************  Load 10 min ceilo file    
        zz10=load([way_proc_data_ceilo, 'CALNEXceilo_10min_a.txt']);
        jdh10=zz10(:,1)+.16666/24;%1    Julian Day, time at the beginning of 10-min ave
         ceilof10=zz10(:,9);%cloudy fraction
        cloudhgt10=zz10(:,11);%median cloud base height
        hgts110=zz10(:,12);%15% clouds lower height
        hgts210=zz10(:,13);%85% clouds lower height
        
      %*******************  Load 15 s ceilo file
        
        
    end;

    micro=1;
if micro
        %CALNEX_MB_sci_day
        %fn2 = [way_proc_data_micro 'vocals_rhb_lwp10min_v1.nc'];
        mby=[];mbyy=[];mb=[];
        fn2 = [way_proc_data_micro 'CALNEX_2010_micro_s.mat'];
%          test = load([way_proc_data_micro 'CALNEX_2010_micro_m.txt']);
        load(fn2);
        mby= CALNEX_2010_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
        mbx=interp1(jdm,mb(m,:),jdx);
        jdmx=mbx(:,1);
        vap=mbx(:,2);%path-integrated water vapor, cm (8.3f)
        liq=mbx(:,3);%path-integrated liquid water, Grams/Meter2 
        ii=find(mb(:,5)>1e3);mb(ii,5)=NaN; ii=find(mb(:,5)<-50);mb(ii,5)=NaN;
%         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 
%         
    end;

%******************   PLOT full time series
        ymaxx=2.0;
        figure
        subplot(3,1,1),imagesc(hhx,ht, SNRx(1:120,:),[-25, 20]);  %
        axis xy;ylim([0 ymaxx]);
        ylabel('Ht (km)')
         hold on
        plot((zz10(:,1)-ddd)*24,cloudhgt10/1e3,'-ok');
        hold off
        %axis([15 15.25 0 3])
        tit = {['VOCALS' year ' W-Band'' YD ' ddds ' Hr ' hr]};
        title(tit)
        h = colorbar;
        axes(h)
        ylabel('SNR dBZ');
       

        subplot(3,1,2),imagesc(hhx,ht, DVx(1:120,:),[-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(3,1,3);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)



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 1.0]);
  %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(['Julian Day (' year ')']); 
  print('-djpeg90 ',[way_cloud_data_wband '\cloud_stats_' num2str(ddd), '.jpg']);

  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(['Julian Day (' year ')']); 
    print('-djpeg90 ',[way_cloud_data_wband '\cloud_thick_' num2str(ddd), '_' hr '.jpg']);
  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]);
 
