% VOCALS 2008 :: 2008-10-16 :: Simon de Szoeke
% Read 1-minute data processed by the Campbell Scientific
% data logger.
disp(['read_pc_day: ', cruise, ' ' yyyy]);
fclose('all');
jd=sprintf('%03i',ddd);

sig_sb=5.67e-8;				        %Stefan Boltzmann constant
jax1=1;
jd_pc=[];st2=[];

for jam=0:23,					    %cycle thru 24 hourly stats files
    hr=sprintf('%02i',jam);   
    e=[way_raw_data_flux, 'day',jd,'\met0' year(3:4) jd hr '_raw.txt'];
    %disp(['Reading means file from hour ',hr]);
    flist=fopen(e,'r');
   
    if flist>0,					    %if the file exists, 
            clear stx;
                                     %1                        %10                            %20                           %30                           %40                           %50                           %60                           %70                 %77
              stx =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 %f %f %f %f %f %f %f %f %f %f %f %f %f %f %*[^\n]', 'delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty', 'Radiometers & Mean Met');
%                 for ii=1:79
%                     if length(stx{1,ii})~=60
%                         stx{1,ii}=[stx{1,ii};ones(60-length(stx{1,ii}),1)*NaN];
%                     end;
%                 end;    
                stx=cell2mat(stx)';
%                 stx(4:79,any(isnan(stx),1))=NaN;      %put NaN into any rows containing NaNs (to remove bad lines, gap...), except JD

                [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 = [jd_pc jd1];
    	        jax1=jax1+nl1;
    end;                            %end if flist
end;                                %end for jam
if length(st2)<1440,
    st2=[st2 ones(79,1440-length(st2))*NaN];
    jd_pc=[jd_pc ones(1,1440-length(jd_pc))*NaN];
end

org_carrier=st2(17,:);         %sti carrier voltage
aspir_on=st2(34,:);            %aspirator operational (+5V) or backflow (0V)

press=st2(36,:);    %pressure (mbar)
psp1=st2(30,:);
psp2=st2(31,:);
org=st2(32,:);
Tsea=st2(25,:); %SPdeS
Tvais=st2(26,:);
Rhvais=st2(27,:); 

% diagnose averages
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']);

disp(sprintf('VOCALS2008 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))

%Do some adjustment by eye-hand custom QCing... !!!!! Every day is different !!!!!
ii=find(Rhvais<0);Rhvais(ii)=NaN;

%do some bias corrections for PIR temperatures.
% !!!!! Every day is different !!!!!
% These are for GASEX08 !
% if ddd==63
%     Tc1=st2(21,:)+0.85;  %correct for bias
%     Td1=st2(22,:)+0.85;
%     Tc2=st2(23,:)+0.89;   
%     Td2=st2(24,:)+0.99;
% end;
% if ddd==67
%     Tc1=st2(21,:)-.5;  %correct for bias
%     Td1=st2(22,:)-.5;
%     Tc2=st2(23,:)-.3;   
%     Td2=st2(24,:)-.2;
% end;
% if ddd==68
%     Tc1=st2(21,:)-.6;  %correct for bias
%     Td1=st2(22,:)-.6;
%     Tc2=st2(23,:)-.42;   
%     Td2=st2(24,:)-.22;    
% end;
% if ddd==71 || 72
%     Tc1=st2(21,:)-.4;  %correct for bias
%     Td1=st2(22,:)-.4;
%     Tc2=st2(23,:)-.2;   
%     Td2=st2(24,:)-.1;
% end;
% if ddd==76
%     Tc1=st2(21,:)-.5;  %correct for bias
%     Td1=st2(22,:)-.5;
%     Tc2=st2(23,:)-.3;   
%     Td2=st2(24,:)-.2;
% end;
% if ddd==77
%     Tc1=st2(21,:)-.3;  %correct for bias
%     Td1=st2(22,:)-.3;
%     Tc2=st2(23,:)-.2;   
%     Td2=st2(24,:);
% end;

% 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,:);

% pir1=therm1+sig_sb*(Tc1+273.16).^4-4*sig_sb*((Td1+273.16).^4-(Tc1+273.16).^4);
% pir2=therm2+sig_sb*(Tc2+273.16).^4-4*sig_sb*((Td2+273.16).^4-(Tc2+273.16).^4);
% 273.16 is the TRIPLE POINT, NOT the conversion from Celsius to Kelvin!  SPdeS
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;

%recompute mean after adjustment
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']);
disp(sprintf('%12s%12.2f%12.2f%12.2f%12.2f','adjusted',avgpsp,avgTvais,avgQvais,avgsnake))

if plotit
    figure;plot(jd_pc,press,'r');title(['Vaisala Pressure. ', cruise yyyy]);xlabel(['JD' yyyy]);ylabel('Pressure (mb)');xlim([ddd,ddd+1]);%axis([ddd,ddd+1,nanmean(press)-5,nanmean(press)+5])
    if prtit; print('-djpeg90 ',[way_images_flux '\Pressure_' num2str(ddd) '.jpg']); 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(['PSD Temperatures. ', cruise yyyy]);legend('case 1','case 2','dome 1', 'dome 2','snake', 'Tair');xlabel(['JD' yyyy]);ylabel('Temperature (C)');xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\Temperatures_' num2str(ddd) '.jpg']); end;
    
    figure;plot(jd_pc,psp1,jd_pc,psp2);title(['PSD psp. ', cruise yyyy]);xlabel(['JD' yyyy]);ylabel('Solar FLux (W/m^2)');legend('PSP1','PSP2');xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\Solar_flux_' num2str(ddd) '.jpg']); end;
    
    figure;plot(jd_pc,pir1,jd_pc,pir2);title(['PSD pir. ', cruise yyyy]);xlabel(['JD' yyyy]);ylabel('IR Flux (W/m^2)');legend('PIR1','PIR2');xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\IR_flux_' num2str(ddd) '.jpg']); end;
    
    figure;plot(jd_pc,Rhvais);title(['PSD RH. ', cruise yyyy]);xlabel(['JD' yyyy]);ylabel('Relative Humidity (%)');xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\RH_' num2str(ddd) '.jpg']); end;
    
    figure;plot(jd_pc,org);title(['PSD org rain rate. ', cruise yyyy]);xlabel(['JD' yyyy]);ylabel('Rain Rate (mm/hr)');xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\rainrate_' num2str(ddd) '.jpg']); end;
    
    figure;plot(jd_pc,aspir_on);title(['PSD aspirator. ', cruise yyyy]);xlabel(['JD' yyyy]);ylabel('Backflow Indicator (Volt). If 0V, there is backflow');axis([ddd,ddd+1,0,6])
    if prtit; print('-djpeg90 ',[way_images_flux '\backflow_indicator_' num2str(ddd) '.jpg']); end;
    
    figure;plot(jd_pc,org_carrier);title(['PSD org carrier. ', cruise yyyy]);xlabel(['JD' yyyy]);ylabel('Rain Gauge Function (V)');axis([ddd ddd+1 0 10])
    if prtit; print('-djpeg90 ',[way_images_flux '\rain_gauge_function_' num2str(ddd) '.jpg']); end;
    
    figure;plot(jd_pc,Tsea);title(['PSD seasnake. ', cruise yyyy]);xlabel(['JD' yyyy]);ylabel('Temperature (degC)');xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\seasnake_' num2str(ddd) '.jpg']); end;
end

if saveit;
  np=length(st2);  
  prt_jas_means;
  prt_jas_rads;
  [r,c] = size(st2');
  
  st2(c+1,:) = jd_pc;
  met = st2';
  save H:\RHB\flux\Raw\met met
end;