% read_co2
% Xcruise 2009 :: 2009-09-07 :: LB

disp(['read_co2: ', cruise]);

fclose('all');
jd=sprintf('%03i',ddd);

% allocate z, but z can still grow
switch licor_prefix
    case {'lis', 'lic'}
        z=NaN+zeros(24*3600*flicor,8); % flicor=10Hz -->864000 samples
    case 'lin'
        z=NaN+zeros(24*3600*flicor,9);
    case {'li1', 'li2', 'li3'}
        z=NaN+zeros(24*3600*flicor,11);
end
lic1=[];
jax1=1; % initialise line index, incremented for each textscan, cumulative for whole day
jstarthour=NaN+zeros(24,1);
jendhour=NaN+zeros(24,1);
for jam=0:23				    %cycle thru 24 hourly files
    hr=sprintf('%02i',jam);
    dfl=[way_raw_data_flux, yyyy(3:4),jd,'\' licor_prefix '0', yyyy(3:4) jd hr '_raw.txt'];
    disp(['Reading Licor file  ',dfl])
    flist=fopen(dfl,'r');
    
    jstarthour(jam+1)=jax1; % save index of first line in hour
    if flist>0          %if the file exists
        %how long is the data message?
        % if ncol==6 then message is:
        % mmsslll status CO2[mmol/m3] H2O[mmol/m3] Temperature[C] pressure[kPa]
        % if ncol==8:
        % mmssll status CO2[umol/mol] CO2[mmol/m3] H2O[mmol/mol] H2O[mmol/m3] Temperature[C] pressure[kPa]
        
        %Next 4 lines changeds as per LBariteau 10/8/2009 DEW
%         fgetl(flist);
%         ncol=length(sscanf(fgetl(flist),'%f'));
%         frewind(flist);
ncol =8;
        while ~feof(flist)
            switch licor_prefix
                case {'lis', 'lic'} %VOCALS 2008 is lic, and has 2 data formats!...
                    if ncol==6
                        temparr=textscan(flist,'%2f%2f%3f %f %f %f %f %f %*[^\n]','delimiter', ', ', 'headerlines', 1, 'emptyvalue', NaN,'treatAsEmpty','Overflow');
                    elseif ncol==8 % ignore data positions 5 and 7 
                        temparr=textscan(flist,'%2f%2f%3f %f %*f %f %*f %f %f %f %*[^\n]','delimiter', ', ', 'headerlines', 1, 'emptyvalue', NaN,'treatAsEmpty','Overflow');
                    end
                case 'lin'
                    temparr=textscan(flist,'%2f%2f%2f%3f %f %f %f %f %f %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty','Overflow');
                case {'li1', 'li2', 'li3'}
                    temparr=textscan(flist,'%2f%2f%2f%3f %f %f %f %f %f %f %f %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty','Overflow');
            end

            nl1=length(temparr); % 8-11
            nr1=length(temparr{nl1}); % <~36000 (use length of LAST nl1 column)
            try
                %z(jstarthour:jendhour,1:nl1)=cell2mat(temparr);
                % more robust in case of incomplete lines:
                for col=1:nl1
                    z(jax1:jax1+nr1-1,col)=temparr{col}(1:nr1);
                end
                jax1=jax1+nr1; % increment line index after each successful textscan
            catch
                disp(sprintf('read_co2: problem reading lines %i-%i.', jax1, jax1+nr1-1))
                pack
            end
        end %end of while ~feof
        % save index of last line in hour after reading all lines in hourly file
        jendhour(jam+1)=jax1-1; % jax1 is already the next line!
        
        lic1=z; % save an unfilled copy
        %fill some bad values of z with previous values, for transforms
        switch licor_prefix % process all types of licor messages
            case {'lis', 'lic'}
                ii=find(z(:,8)<50);
                z(ii,7:8)=z(ii-1,7:8);
            case 'lin'
                ii=find(z(:,9)<50);
                z(ii,8:9)=z(ii-1,8:9);
            case {'li1', 'li2', 'li3'}
                ii=find(z(:,11)<50);
                z(ii,10:11)=z(ii-1,10:11);
        end

        if ismember(jam,[00 12 18]) && ~isempty(z(~isnan(z(:,ncol)),ncol)) % spectra for selected hours
            disp(['Calculating licor humidity spectra for ',hr])
            switch licor_prefix
                case {'lis', 'lic'}
                    zhour=z(jstarthour(jam+1):jendhour(jam+1),6)*18/1000; %g/m3
                case 'lin'
                    zhour=z(jstarthour(jam+1):jendhour(jam+1),7)*18/1000; %g/m3
                case {'li1' ,'li2' ,'li3'}
                    zhour=z(jstarthour(jam+1):jendhour(jam+1),9)*18/1000; %g/m3
            end;
            spechum_lic=zhour(isfinite(zhour));            
            [S_Lic,F]=psd2(detrend(spechum_lic),length(spechum_lic),flicor);
            [S_Lics,Fsl]=specsmoo(S_Lic,flicor);

            if plotit
                figure;loglog(Fsl,S_Lics.*Fsl,Fsl(length(Fsl)-31:length(Fsl)),((Fsl(length(Fsl)-31:length(Fsl))).^(-2/3))./(Fsl(length(Fsl)-18).^(-2/3)).*S_Lics(length(Fsl)-18).*Fsl(length(Fsl)-18));
                title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor hygrometer (unit: %s),  spectrum: hour %02i',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,licor_prefix,jam),'FontWeight','Bold','Interpreter','none')
                xlabel('Frequency (Hz)');
                ylabel('Spectrum * Frequency (g^2 m^-^6)');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_Licor_H2O_spectra_hour' hr '_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end
            end
            if saveit && jam==12 % only save for 12Z
                specfile=[way_spec_data_flux,cruise '_PSD_' licor_prefix '_q_sp_hour' hr '_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd),'.txt'];
                fspec=fopen(specfile,'w');
                specs=[Fsl;S_Lics];
                fprintf(fspec,'%11.5f %11.5f \r\n',specs);
                fclose(fspec);
                clear fspec specs specfile
            end
            clear S_Lic F S_Lics Fsl spechum_lic
        end; % spectra only for 00 12 18 hours
        fclose(flist); % close the files!
    else
        disp(sprintf('\b<--File not found.'))
    end           % end if flist
end % for jam=0:23

% % 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(lic1(1,:)*10000*100e-9/86400));
% hh = str2num(dates(:,13:14));
% mm = str2num(dates(:,16:17));
% ss = str2num(dates(:,19:20));
% jd_lic= 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_lic(1) = jd_lic(1)-1;
% end

% How did we do? Diagnose jstarthour,jendhour:
%[jstarthour jendhour jendhour-jstarthour [0; jendhour(1:end-1)-jstarthour(2:end)]]
%should be jendhour(i)+1=jstarthour(i+1), but not so!

% truncate data and give sensible names
if ~isempty(lic1) && ~isnan(max(jendhour))
    switch licor_prefix
        case {'lis', 'lic'}
            Licor_diag =lic1(1:max(jendhour),4);
            Licor_CO2  =lic1(1:max(jendhour),5);%mmol/m^3
            Licor_H2O  =lic1(1:max(jendhour),6);%mmol/m^3
            Licor_Temp =lic1(1:max(jendhour),7);%C
            Licor_Press=lic1(1:max(jendhour),8);%kPa
        case 'lin'
            Licor_diag =lic1(1:max(jendhour),5);
            Licor_CO2  =lic1(1:max(jendhour),6);%mmol/m^3
            Licor_H2O  =lic1(1:max(jendhour),7);%mmol/m^3
            Licor_Temp =lic1(1:max(jendhour),8);%C
            Licor_Press=lic1(1:max(jendhour),9);%kPa
        case {'li1','li2','li3'}
            Licor_diag =lic1(1:max(jendhour),5);
            Licor_CO2  =lic1(1:max(jendhour),7);%mmol/m^3
            Licor_H2O  =lic1(1:max(jendhour),9);%mmol/m^3
            Licor_Temp =lic1(1:max(jendhour),10);%C
            Licor_Press=lic1(1:max(jendhour),11);%kPa
    end;
    % create yearday time array
    switch licor_prefix
        case {'lis', 'lic'}
            jd1=ddd+(...
                (lic1(:,1)+(lic1(:,2)+lic1(:,3)/1000)/60)/60)/24;
            for jam=1:24
                if isfinite(jendhour(jam))
                    jd1(jstarthour(jam):jendhour(jam))=jd1(jstarthour(jam):jendhour(jam))+(jam-1)/24;
                end
            end
        case {'lin', 'li1', 'li2','li3'}
            jd1=ddd+(lic1(:,1)+...
                (lic1(:,2)+(lic1(:,3)+lic1(:,4)/1000)/60)/60)/24;
    end
    jd_lic=jd1(1:max(jendhour)); % truncate

    % interpret status diagnostic
    xx=Licor_diag;
    b7=floor(xx/128);
    y=b7*128;
    b6=floor((xx-y)/64);
    y=y+b6*64;
    b5=floor((xx-y)/32);
    y=y+b5*32;
    b4=floor((xx-y)/16);
    y=y+b4*16;
    agc=(xx-y)*6.25;

    % QC values
%     Licor_Press(Licor_Press<=0 | Licor_Press>200)=NaN;
%     Licor_Temp(Licor_Temp==0 | Licor_Temp>1000)=NaN;
%     Licor_CO2(Licor_CO2>100 | Licor_CO2<0)=NaN;
%     Licor_H2O(Licor_H2O>1500 | Licor_H2O<0)=NaN; % threshold raised from 500 in GASEX to 1500 in VOCALS

    % convert units from mol/m3 to mol/mol
    Rgas=287.1;
    Rgas_universal=8.314472; % Pa m3 /K /mol
    Licor_CO2_a=(Licor_CO2*1e-3).*Rgas_universal.*(Licor_Temp+273.15)./(Licor_Press*1e3); %mmol/m3 to mol/mol
    Licor_H2O_a=(Licor_H2O*1e-3).*Rgas_universal.*(Licor_Temp+273.15)./(Licor_Press*1e3); %mmol/m3 to mol/mol
    %micromol/mol=ppmv

    Mw=18.01528; % molar mass of H20 g/mol
    Md=28.964;
    epsilon=Mw/Md; % mass of water to mass of dry air 0.622
    rhoa_lic=(Licor_Press*1e3).*(1-Licor_H2O_a*(1-epsilon))./(Rgas.*(Licor_Temp+273.15));%air density in licor [kg/m3]
    Licor_H2O_m=(Licor_H2O*1e-3).*Mw./rhoa_lic; %mmol/m3 to g/kg
    %Licor_H2O_q=Licor_H2O_a*epsilon; % mol/mol to kg/kg (note different units to above)

    if plotit
        figure;plot(jd_lic,Licor_H2O_m);
        xlabel('Hour (UTC)');ylabel('Specific humidity (g/kg)');
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor H2O (unit: %s)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,licor_prefix),'FontWeight','Bold','Interpreter','none')
        axis([ddd ddd+1 median(Licor_H2O_m(~isnan(Licor_H2O_m)))-10 median(Licor_H2O_m(~isnan(Licor_H2O_m)))+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_Licor_H2O_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end

        figure;plot(jd_lic,Licor_CO2_a*1e6);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor CO2 (unit: %s)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,licor_prefix),'FontWeight','Bold','Interpreter','none')
        xlabel('Hour (UTC)');ylabel('CO2 (\mumol/mol)');
        axis([ddd ddd+1 median(Licor_CO2_a(~isnan(Licor_CO2_a))*1e6)-10 median(Licor_CO2_a(~isnan(Licor_CO2_a))*1e6)+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_Licor_CO2_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end

        figure;plot(jd_lic,Licor_Temp);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor box temperature (unit: %s)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,licor_prefix),'FontWeight','Bold','Interpreter','none')
        xlabel('Hour (UTC)');ylabel('Temperature (degC)');
        axis([ddd ddd+1 median(Licor_Temp(~isnan(Licor_Temp)))-10 median(Licor_Temp(~isnan(Licor_Temp)))+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_Licor_Temp_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end

        figure;plot(jd_lic,Licor_Press*10);
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor box pressure (unit: %s)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,licor_prefix),'FontWeight','Bold','Interpreter','none')
        xlabel('Hour (UTC)');ylabel('Pressure (hPa)')
        axis([ddd ddd+1 median(Licor_Press(~isnan(Licor_Press))*10)-10 median(Licor_Press(~isnan(Licor_Press))*10)+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_Licor_Pressure_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end

        figure;plot(jd_lic,agc,[ddd ddd+1],[55 55],'r',[ddd ddd+1],[65 65],'r'); 
        title(sprintf('%s (%04i-%02i-%02i, DOY%03i).  Licor AGC value (unit: %s)',cruise,Vdate(1),Vdate(2),Vdate(3),ddd,licor_prefix),'FontWeight','Bold','Interpreter','none')
        axis([ddd ddd+1 0 100]);
        xlabel('Hour (UTC)');ylabel('AGC (%)')
        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_Licor_AGC_' sprintf('%04i_%02i_%02i_%03i',Vdate(1),Vdate(2),Vdate(3),ddd) graphformat]); end
    end;
end

if  exist('Licor_H2O_m','var') && ~isempty(Licor_H2O_m)
    jd_pc=ddd+(0:1439)/60/24;
    Licor_H2O_mi=intavg(jd_lic,Licor_H2O_m,jd_pc); %g/kg
    Licor_CO2_ai=intavg(jd_lic,Licor_CO2_a*1e6,jd_pc); %umol/mol
    Licor_Tempi=intavg(jd_lic,Licor_Temp,jd_pc); %degC
    Licor_Pressi=intavg(jd_lic,Licor_Press*10,jd_pc); %hPa
    Licor_agci=intavg(jd_lic,agc,jd_pc); 
    
    if saveit;
        eval(['prt_jas_co2_' currentMname]);
    end;
end;



if clearit
    clear Licor_H2O_m jd_lic Licor_CO2_a Licor_Temp Licor_Press agc
    clear Licor_diag rhoa_lic Licor_H2O Licor_CO2 Licor_H2O_a
    clear b4 b5 b6 b7 bb bbn jd1 lic1 xx y z zhour
end
