% Xcruise 2009 :: 2009-09-07 :: LB
% Read 1-minute data processed by the Campbell Scientific
% data logger.
disp(['read_met_day: ', cruise]);
fclose('all');
jd=sprintf('%03i',ddd);

sig_sb=5.67e-8;				        %Stefan Boltzmann constant
jax1=1;

% Declare Arrays
nr1=65;
st2=ones(nr1,1440)*NaN;
jd_pc=ones(1,1440)*NaN;

for jam=0:23,					    %cycle thru 24 hourly stats files
    hr=sprintf('%02i',jam);   
    e=[way_raw_data_flux, yyyy(3:4),jd,'\met0' yyyy(3:4) jd hr '_raw.txt'];
    disp(['Reading means file  ',e]);
    flist=fopen(e,'r');
   
    if flist>0,					    %if the file exists, 
        stx=[];
        while ~feof(flist) % end-of-file not reached         
            try
                                     %1                        %10                            %20                           %30                           %40                           %50                           %60                         
              temp =textscan(flist,'%2f%2f%3f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %*[^\n]', 'delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN);
              stx=[stx cell2mat(temp)'];
            catch 
                [lnr1,cnr1]=size(temp{1,nr1});
                if lnr1==0 || cnr1==0; continue; end % read next batch of lines
                for ll=1:nr1
                    if length(temp{1,ll})~=lnr1
                        temp{1,ll}(lnr1+1)=[];
                    end;
                end;
                if isempty(temp{1,1});
                    temp{1,nr1}=temp{1,nr1}';
                end;
                stx=[stx cell2mat(temp)'];
            end;   %end try/catch
        end; %end ~feof(flist)
        
        [nr1,nl1]=size(stx);
        st2(:,jax1:jax1+nl1-1)=stx;
        jd1=ddd + (jam+(stx(1,:)+(stx(2,:)+stx(3,:)/1000)/60)/60)/24;
        jd_pc(:,jax1:jax1+nl1-1)=jd1;
        jax1=jax1+nl1;
    fclose(flist);
    else
        stx=ones(nr1,65)*NaN;
        [nr1,nl1]=size(stx);
        st2(:,jax1:jax1+nl1-1)=stx;
        jax1=jax1+nl1;
    end;                            %end if flist
end;                                %end for jam

org_carrier=st2(19,:);         %sti carrier voltage  ~4.9 VDC
aspir_on=st2(34,:);            %aspirator operational (+5V) or backflow (0V)
                               %or could be tge motor revolution (in Hz?)
press=st2(36,:);    %pressure (mbar)
psp1=st2(30,:);
psp2=st2(31,:);
org=st2(32,:);
Tsea=st2(25,:); %snake
Tvais=st2(26,:);
Rhvais=st2(27,:); 


% Flag bad values with comments
% DO NOT ADJUST DATA HERE, do SeaSnake T offset, etc. in timatch_5 and
% timatch_30.
Rhvais(Rhvais<0)=NaN;

% case and dome temperatures for PIR1 and PIR2
Tc1=st2(21,:);  % unadjusted
Td1=st2(22,:);
Tc2=st2(23,:);
Td2=st2(24,:);
% thermopile voltages
therm1=st2(28,:);
therm2=st2(29,:);

C2K=273.15;
pir1=therm1+sig_sb*(Tc1+C2K).^4-4*sig_sb*((Td1+C2K).^4-(Tc1+C2K).^4);
pir2=therm2+sig_sb*(Tc2+C2K).^4-4*sig_sb*((Td2+C2K).^4-(Tc2+C2K).^4);

pspm=(psp1+psp2)/2;

%check for bad values in PIR
ij=find(pir1<50 | pir1>500);pir1(ij)=NaN;
ik=find(pir2<50 | pir2>500);pir2(ik)=NaN;
pirm=(pir1+pir2)/2;

% % Time in flux data files is the number of 100-nano second intervals since
% % Jan 1, 1601 divided by 10000. It was chosen to maintain compatibility
% % with an HP-UNIX version of the code form long ago
% %  0             'dd-mmm-yyyy HH:MM:SS'  
% %01-Mar-2000 15:45:17
% dates = datestr(datenum(1601,1,1) + datenum(st2(1,:)*10000*100e-9/86400));
% hh = str2num(dates(:,13:14));
% mm = str2num(dates(:,16:17));
% ss = str2num(dates(:,19:20));
% jd_pc= ddd + (hh+((mm+(ss/60))/60))/24;
% 
% % test to see if first point is from the previous day
% if str2num(dates(1,1:2)) ~= str2num(dates(2,1:2))
%     jd_pc(1) = jd_pc(1)-1;
% end    

% jd_pc2=str2num(jd)+((st2(1,:)-st2(1,1))/1000/60/60)/24;

% DO NOT DO ADJUSTMENT HERE

%Compute mean
avgsnake=nanmean(Tsea);
avgpsp=nanmean(psp1);
avgrain=nanmean(org);
avgTvais=nanmean(Tvais);
avgRHvais=nanmean(Rhvais);
avgQvais=avgRHvais/100*qsea(avgTvais);
qvais=qair_p([Tvais' Rhvais' press']);
% Adjustment is done in timatch_5 and timatch_30 for NTAS 2009

disp(sprintf([cruise yyyy 'day %03i PSD Campbell average meteorology'], ddd))
disp('-----------------------------------------------------')
disp(sprintf('%12s','','PSP','Tvais','Qvais','Tsnake'))
disp(sprintf('%12s%12.2f%12.2f%12.2f%12.2f','raw',avgpsp,avgTvais,avgQvais,avgsnake))

if plotit &&  ~isempty(st2(~isnan(st2)))  %means all NaNs, no data to plot
    figure;plot(jd_pc,press,'r');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Vaisala pressure',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel('Atmospheric Pressure (mb)');axis([ddd ddd+1 median(press(~isnan(press)))-10  median(press(~isnan(press)))+10]);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    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_images_flux cruise '_PSD_Vaisala_Pressure_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
    
    figure;plot(jd_pc,Tc1,'r',jd_pc,Tc2,'m',jd_pc,Td1,'k',jd_pc,Td2,'c',jd_pc,Tsea,'b',jd_pc,Tvais,'g');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD temperatures',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    legend('case 1','case 2','dome 1', 'dome 2','Seasnake', 'Tair','Location','Best');
%     Tscale=median([Tvais(~isnan(Tvais)) Tsea(~isnan(Tsea)) Td2(~isnan(Td2)) Td1(~isnan(Td1)) Tc2(~isnan(Tc2)) Tc1(~isnan(Tc1))]);
%     Tscalestd=std([Tvais(~isnan(Tvais)) Tsea(~isnan(Tsea)) Td2(~isnan(Td2)) Td1(~isnan(Td1)) Tc2(~isnan(Tc2)) Tc1(~isnan(Tc1))]);
%     axis([ddd ddd+1 Tscale-10  Tscale+10])
    xlabel('Hour (UTC)');ylabel('Temperature (degC)');xlim([ddd ddd+1])
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    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_images_flux cruise '_PSD_Temperatures_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
    
    figure;plot(jd_pc,psp1,jd_pc,psp2);
    xlabel('Hour (UTC)');ylabel('Solar Flux (W/m^2)');legend('PSP1','PSP2','Location','Best');
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Precision Spectral Pyranometers',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    axis([ddd ddd+1 -5 1200]);
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    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_images_flux cruise '_PSD_Eppley_Solarflux_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
    
    figure;plot(jd_pc,pir1,jd_pc,pir2);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Precision Infrared Radiometers',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel('IR Flux (W/m^2)');legend('PIR1','PIR2','Location','Best');xlim([ddd ddd+1])
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    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_images_flux  cruise '_PSD_Eppley_IRflux_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
    
    figure;plot(jd_pc,Rhvais);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Relative Humidity',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel('Relative Humidity (%)');
    axis([ddd ddd+1 0 100 ])
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))   
    annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
    grid
    if prtit; print(graphdevice,[way_images_flux cruise '_PSD_Vaisala_RH_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
    
    figure;plot(jd_pc,org);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD Optical Rain Gauge',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel('Rain Rate (mm/hr)');xlim([ddd ddd+1])
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    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_images_flux cruise '_PSD_ORG_rainrate_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
    
    figure;plot(jd_pc,aspir_on);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD T/RH aspirator fan',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel('Aspirator Fan Speed (rpm)');xlim([ddd,ddd+1])
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    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_images_flux cruise '_PSD_backflow_indicator_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
    
    figure;plot(jd_pc,org_carrier);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD ORG signal carrier',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel('Rain Gauge Voltage (volt)');axis([ddd ddd+1 4.0 5.0])
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    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_images_flux cruise '_PSD_ORG_rain_gauge_function_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
    
    figure;plot(jd_pc,Tsea);
    title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  PSD seasnake',cruise,Vdate(1),Vdate(2),Vdate(3),ddd),'FontWeight','Bold','Interpreter','none')
    xlabel('Hour (UTC)');ylabel('Sea Surface Temperature (degC)');xlim([ddd ddd+1])
    set(gca(gcf),'XTick',ddd:2/24:ddd+1);datetick('x','HH:MM','keepticks');
    xax=get(gca(gcf),'XTickLabel');xax(end,:)='24:00';set(gca(gcf),'XTickLabel',xax(:,1:2))    
    grid
    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_images_flux cruise '_PSD_seasnake_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end;
end

if saveit;
  if length(isnan(jd_pc)==1)==1440;   %means no data files detected during the 24H period
        jd_pc=ddd+(0:1439)/60/24;       %creates a fake decimal DOY base so that timatch programs work properly
  end;
  np=length(st2);  
  eval([ 'prt_jas_means_', currentMname])
  eval([ 'prt_jas_rads_', currentMname])
%   [r,c] = size(st2');
end;