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

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

% Declare Arrays
nr1=57;
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'];
    
    flist=fopen(e,'r');
   
    if flist>0,					    %if the file exists, 
        disp(['Reading means file from hour ',hr]);
        stx=[];
        while ~feof(flist) % end-of-file not reached         
            try
                                     %1                           %10                          %20                           %30                           %40                           %50                                                    
%               temp =textscan(flist,'%2f%2f%3f %2f/%2f/%4f %2f:%2f:%2f %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)'];
                %reads in as %s as potential NAN in the met files
                temp =textscan(flist,'%2f%2f%3f %2f/%2f/%4f %2f:%2f:%2f %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %*[^\n]', 'delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN);
                stx=[stx temp{1} temp{2} temp{3} temp{4} temp{5} temp{6} temp{7} temp{8} temp{9}];
                for ll=10:length(temp)
                        temp{ll}(ismember(temp{ll},'NAN'))=cellstr('NaN');
                        stx=[stx str2num(char(temp{ll}))];
                end;
%                 lud=find((ismember(temp{ll},''))); 
% [temp{1}(lud) temp{2}(lud) temp{3}(lud) ]
%sometimes a comma is
%                 missing thus offsetting the reading and empty value are
%                 found which crashes the code - LB - 10192013
%                 will edit manually the file for now...

            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;           
                stx2=[];
                stx2=[stx2 temp{1} temp{2} temp{3} temp{4} temp{5} temp{6} temp{7} temp{8} temp{9}];
                for ll=10:length(temp)
                        temp{ll}(ismember(temp{ll},'NAN'))=cellstr('NaN');
                        stx2=[stx2 str2num(char(temp{ll}))];
                end;
            end;   %end try/catch
        end; %end ~feof(flist)
        stx=stx';
        [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
        disp(['Reading means file from hour ',hr]);
         disp(sprintf('\b<--File not found.'))
        stx=ones(nr1,60)*NaN;
        [nr1,nl1]=size(stx);
        jd1=ddd + (jam+(0:59)/60)/24;
        jd_pc(:,jax1:jax1+nl1-1)=jd1;
        st2(:,jax1:jax1+nl1-1)=stx;
        jax1=jax1+nl1;
    end;                            %end if flist
end;                                %end for jam

org_carrier=st2(20,:);         %sti carrier voltage  ~4.9 VDC
aspir_on=st2(26,:);            %motor revolution (in rpm)
press=st2(25,:);    %pressure (mbar)
psp1=-st2(17,:);
psp2=-st2(19,:);
org=st2(22,:);
Tsea=st2(24,:); %snake (4W)

% Tvais=st2(26,:);
% Rhvais=st2(27,:); 

%removes offset on new ORG
% org=org-min(org);
% %also sets remaining drifts to zero
% org(abs(diff(org))<=0.031)=0;


% 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(5,:)+0;  % adjusted
Td1=st2(7,:)+0;
Tc2=st2(9,:)+0;
Td2=st2(11,:)+0;
% thermopile voltages
therm1=-st2(13,:);
therm2=-st2(15,:);

C2K=273.15;
pir1=therm1+sig_sb*(Tc1+C2K).^4-3*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;


%%%%%%%%%%%%%%%%%%%%%%%%  reads in T/RH
read_trh_WHOTS_2014;

if length(RH(isnan(RH)))~=86400;   %means no data files detected during the 24H period
    Rhvais=interp1(jd1(~isnan(jd1)),RH(~isnan(jd1)),jd_pc);          
    Tvais=interp1(jd1(~isnan(jd1)),T(~isnan(jd1)),jd_pc);         
else
    Rhvais=ones(1,1440)*NaN;          
    Tvais=ones(1,1440)*NaN;         
end;
  
  
% 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 ' / 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))


% jd_pc=jd_mean1;
if plotit 
    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]);
%      ylim([1000 1020])
    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,Td1,'k',jd_pc,Tsea,'b',jd_pc,Tvais,'g',jd_pc,Tc2,'rx',jd_pc,Td2,'kx');
    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','dome 1', 'Seasnake', 'Tair','case 2','dome 2','Location','SouthWest');
%     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])
    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 40 105])
    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]);ylim([0 80])
    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.5])
    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(jd_pc(isnan(jd_pc)))==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;