% 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'); %cd E:\Data\DYNAMO_2011\Revelle\Scientific_analysis\programs\Flux %################## Cruise informations ################################## cruise='DYNAMO'; % Cruise name year='2011'; % year (deprecated--make year into a double) yyyy='2011'; % year (string) ship='Revelle'; % Research vessel big_dir='D:\Data'; %################## W-band radar files ################################### way_raw_data_wband=[big_dir '\DYNAMO_2011\' ship '\radar\wband\raw\']; way_raw_images_wband=[big_dir '\DYNAMO_2011\' ship '\radar\wband\raw_images\']; way_raw_motion_wband=[big_dir '\DYNAMO_2011\' ship '\radar\wband\raw']; way_cloud_data_wband=[big_dir '\DYNAMO_2011\' ship '\radar\wband\processed']; %################## Ceilometer files ##################################### way_images_ceilo=[big_dir '\DYNAMO_2011\' ship '\ceilometer\proceessed\']; % daily saved impates way_raw_data_ceilo=[big_dir '\DYNAMO_2011\' ship '\ceilometer\raw\']; % Raw Data way_proc_data_ceilo=[big_dir '\DYNAMO_2011\' ship '\ceilometer\processed\']; % Processed Data %################## Radiometer files ##################################### way_proc_data_micro=[big_dir '\DYNAMO_2011\' ship '\radiometer\mailbox\processed\']; % Raw Data %%%%%%%%%%%%%%%%%%%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 ################################## Nrange=220; RangeRes=0.031; motionread=0;%set to zero to skip motion correction ht =single((1:1:Nrange)*RangeRes); % range gates in km zmxx=single(ht(Nrange)); dbz0=-43; %%%%%%%%%%%%%%%%%%%%%%%%% 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=1; 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];W % end film=dr(jamn).name fn = [way_raw_data_wband yyyy ddds hr film(10:11) 'MMCRMom.nc'];%[way_raw_data_wband film]; FID = fopen(fn,'r'); if FID > 0 % Test to see if file exists jamn=jamn+1; 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=single(Heights);%netcdf.getVar(FID,27);% 27 ht=single(z(1,:)/1e3); RangeRes=ht(2)-ht(1);dbz0=-43;if RangeRes>0.04;dbz0=-51;end; % nrg=netcdf.getVar(FID,18); Nrange=single(NumHeights(1)); % Ref=netcdf.getVar(FID,40); %40 % SNR=netcdf.getVar(FID,36); %36 % SW=netcdf.getVar(FID,38); % 38 % DVa=netcdf.getVar(FID,35); %35 fclose(FID); Ref = single(Reflectivity)'; SW = single(SpectralWidth)'; DVa = single(MeanDopplerVelocity)'; SNR = single(SignalToNoiseRatio)'; 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; pause(5); 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'-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 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: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']); 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 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: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']); %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=load([way_proc_data_ceilo cruise '_ceilo_h_a.txt']); else zz=NaN*ones(1,length(Ref)); end; 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 if ceilo %******************* Load 1 hr ceilo file zz10=load([way_proc_data_ceilo cruise '_ceilo_10min_a.txt']); else zz10=NaN*ones(1,length(Ref)); end; 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 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_mb 'DYNAMO_2011_micro_s.txt']; % test = load([way_proc_data_micro 'CALNEX_2010_micro_m.txt']); mby=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=mbx(:,1); vap=mbx(:,2);ii=find(vap<0.1);vap(ii)=NaN;%path-integrated water vapor, cm (8.3f) liq=mbx(:,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;if RangeRes>0.04;ymaxx=ht(Nrange);end; ii=1:2:length(hhx); figure(3) %subplot(2,1,1),imagesc(hhx,ht, SNRx(1:Nrange,:),[-25, 20]); % subplot(2,1,1),imagesc(hhx(ii),ht, SNRx(1:Nrange,ii),[-25, 20]); % axis xy;ylim([0 ymaxx]); ylabel('Ht (km)') hold on plot((zz10(:,1)-ddd)*24,cloudhgt10/1e3,'ok','MarkerSize',3); hold off %axis([15 15.25 0 3]) tit = {[cruise year ' W-Band YD ' ddds ]}; %tit = {[cruise year ' W-Band'' YD ' ddds ' Hr ' hr]}; title(tit) h = colorbar; h2=get(gca,'position'); axes(h) ylabel('SNR 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('Dynamo 2011 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']); 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(4); 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(['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(5); 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']); 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'2); wmm=mean(mean(DV(ij,kk:kk+dk)'));wdm=mean(mean(SW(ij,kk:kk+dk)')); figure(8); 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(8); 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(8); 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]);