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

sig_sb=5.67e-8;				        %Stefan Boltzmann constant

eval([ 'read_mean1_campbell_', currentMname]);% Reads Sea Snake, ORG and airflow
eval([ 'read_mean2_campbell_', currentMname]);% Reads PSP and PIRs (Eppley and K&Z bridge) and BP          
eval([ 'read_trh_', currentMname]);% Reads T/RH     
        
org_carrier=mean1(12,:);         %sti carrier voltage
aspir_on=mean1(15,:);            %rpm
org=mean1(14,:);   %mm/hr
org(org<0)=0;  %removes not rain detection
Tsea=mean1(17,:); %snake 4 wires used for stratus2011 (degC)

Tvais=trhdata(4,:);
Rhvais=trhdata(5,:); 

press=mean2(26,:);    %pressure (mbar=hPa)
psp1=mean2(18,:);  %Eppley PSP1
psp2=mean2(22,:);  %Eppley PSP2

% 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 PIR Eppley
Tc1=mean2(8,:);  %case temperature for PIR Eppley 1 (bridge)--  checked ok 7/8/2013
Td1=mean2(10,:);  %dome temperature for PIR Eppley 1 (bridge)
Tc2=mean2(12,:);  %case temperature for PIR Eppley 2 (bridge)
Td2=mean2(14,:);  %dome temperature for PIR Eppley 2 (bridge)

% thermopile voltages
therm1=mean2(16,:); %Eppley PIR #1 (bridge)
therm2=mean2(20,:); %Eppley PIR #2 (bridge)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Synchronization done here.

jd_pc=ddd+(0:1439)/60/24;
   
ii=find(~isnan(jd_mean2));
if~isempty(ii);
    press=intavg(jd_mean2(ii),press(ii),jd_pc)';
    psp1=intavg(jd_mean2(ii),psp1(ii),jd_pc)';
    psp2=intavg(jd_mean2(ii),psp2(ii),jd_pc)';
    Tc1=intavg(jd_mean2(ii),Tc1(ii),jd_pc)';
    Td1=intavg(jd_mean2(ii),Td1(ii),jd_pc)';
    Tc2=intavg(jd_mean2(ii),Tc2(ii),jd_pc)';
    Td2=intavg(jd_mean2(ii),Td2(ii),jd_pc)';
    therm2=intavg(jd_mean2(ii),therm2(ii),jd_pc)';
    therm1=intavg(jd_mean2(ii),therm1(ii),jd_pc)';
else
    press=NaN*ones(1,length(jd_pc));
    psp1=NaN*ones(1,length(jd_pc));
    psp2=NaN*ones(1,length(jd_pc));
    Tc1=NaN*ones(1,length(jd_pc));
    Td1=NaN*ones(1,length(jd_pc));
    Tc2=NaN*ones(1,length(jd_pc));
    Td2=NaN*ones(1,length(jd_pc));
    therm2=NaN*ones(1,length(jd_pc));
    therm1=NaN*ones(1,length(jd_pc));
end;    
    
ii=find(~isnan(jd_mean1));
if~isempty(ii);
    org_carrier=intavg(jd_mean1(ii),org_carrier(ii),jd_pc)';       
    aspir_on=intavg(jd_mean1(ii),aspir_on(ii),jd_pc)';             
    org=intavg(jd_mean1(ii),org(ii),jd_pc)';   
    Tsea=intavg(jd_mean1(ii),Tsea(ii),jd_pc)';   
%     Tvais=intavg(jd_mean1(ii),Tvais(ii),jd_pc)';   
%     Rhvais=intavg(jd_mean1(ii),Rhvais(ii),jd_pc)';   
else
    org_carrier=NaN*ones(1,length(jd_pc));    
    aspir_on=NaN*ones(1,length(jd_pc));        
    org=NaN*ones(1,length(jd_pc));  
    Tsea=NaN*ones(1,length(jd_pc));
%     Tvais=NaN*ones(1,length(jd_pc));   
%     Rhvais=NaN*ones(1,length(jd_pc));
end;

% ii=find(~isnan(jdtrh));
% if~isempty(ii);
%      Tvais=intavg(jdtrh(ii),Tvais(ii),jd_pc)';   
%      Rhvais=intavg(jdtrh(ii),Rhvais(ii),jd_pc)';   
% else
%      Tvais=NaN*ones(1,length(jd_pc));   
%      Rhvais=NaN*ones(1,length(jd_pc));
% end;
% 
% jd_pc=ddd+(0:1439)/60/24; 
Tvais(isnan(jdtrh))=[];Rhvais(isnan(jdtrh))=[];jdtrh(isnan(jdtrh))=[];
RHm=intavg(jdtrh,Rhvais,jd_pc);        
Tm=intavg(jdtrh,Tvais,jd_pc);         
clear Rhvais Tvais
Rhvais=RHm';Tvais=Tm';
clear Tm RHm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C2K=273.15;
pir1=therm1+sig_sb*(Tc1+C2K).^4-4.2*sig_sb*((Td1+C2K).^4-(Tc1+C2K).^4);%dome coeffs adjusted by cwf 7/8/2013
pir2=therm2+sig_sb*(Tc2+C2K).^4-3.8*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))

% jd_pc=jd_mean1;
if plotit 
    if ~isempty(press(~isnan(press))) && exist('press','var')
    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;
    end;
    
    figure;plot(jd_pc,Tc1,'r',jd_pc,Td1,'k',jd_pc,Tsea,'b',jd_pc,Tvais,'g',jd_pc,Tc2,'--',jd_pc,Td2,'--');
    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','North');
%     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]);ylim([20 35])
    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 Eppley (bridge)','PSP2 Eppley (bridge)','Location','North');
    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 Eppley (bridge)','PIR2 Eppley (bridge)','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;
    subplot(2,1,1);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
    subplot(2,1,2);plot(jd_pc,cumsum(org)/60);
        xlabel('Hour (UTC)');ylabel('Accumulated Rain (mm)');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');            
    orient tall
    if prtit; print(graphdevice,[way_images_flux cruise '_PSD_ORG_' 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(jd_pc);  
  eval([ 'prt_jas_means_', currentMname])
  eval([ 'prt_jas_rads_', currentMname])
%   [r,c] = size(st2');
end;