% read_co2
% VOCALS 2008 :: 2008-10-16 :: Simon de Szoeke
% Does VOCALS have CO2 bulk or covariance flux measurement?
% pC02 of seawater is 350ppmv off Panama.
% We will use shipboard ocean CO2 of Bob Castle for our bulk analysis.

% Error on 20081016: all times but hour zero are filled with NaN.
% for testing:
%  plotit=true;saveit=false;prtit=false;

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

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

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, 'day',jd,'\' licor_prefix '0', year(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]
        fgetl(flist);
        ncol=length(sscanf(fgetl(flist),'%f'));
        frewind(flist);
        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]) % 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([licor_prefix, ' Licor Hygrometer Spectrum:   Day ',jd,'   Hour ',hr]);
                xlabel('Frequency in Hz');
                ylabel('Spectrum * Frequency');
                if prtit; print('-djpeg90 ',[way_images_flux, licor_prefix, '_LicH2O_spectrum' jd hr '.jpg']); end
            end
            if saveit && jam==12 % only save for 12Z
                specfile=[way_spec_data_flux,licor_prefix,'_q_sp',jd,hr,'.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)
    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); set(gca,'ylim',[0 30])
        title([licor_prefix, ' Licor specific humidity. ', cruise, year]);xlabel(['yearday ', year]);ylabel('q (g/kg)');axis([ddd ddd+1 -1 30]);
        if prtit; print('-djpeg90 ',[way_images_flux, licor_prefix, '_Licor_spechum_' num2str(ddd) '.jpg']); end

        figure;plot(jd_lic,Licor_CO2_a*1e6);title([licor_prefix, ' Licor CO2. ', cruise, year]);xlabel(['yearday ', year]);ylabel('CO2 (\mumol/mol)');axis([ddd ddd+1 0 800])
        if prtit; print('-djpeg90 ',[way_images_flux, licor_prefix, '_Licor_CO2_' num2str(ddd) '.jpg']); end

        figure;plot(jd_lic,Licor_Temp);title([licor_prefix, ' Licor Box Temperature [degC]. ', cruise, year]);xlabel(['yearday ', year]);xlim([ddd ddd+1]);ylim([20 40])
        if prtit; print('-djpeg90 ',[way_images_flux, licor_prefix, '_Licor_temp_' num2str(ddd) '.jpg']); end

        figure;plot(jd_lic,Licor_Press*10);title([licor_prefix, ' Licor Box Pressure [hPa]. ', cruise, year]);xlabel(['yearday ', year]);xlim([ddd ddd+1]);ylim([1000 1020])
        if prtit; print('-djpeg90 ',[way_images_flux, licor_prefix, '_Licor_press_' num2str(ddd) '.jpg']); end

        figure;plot(jd_lic,agc,[ddd ddd+1],[55 55],'r',[ddd ddd+1],[65 65],'r'); title([licor_prefix, ' AGC value. Typical "clean window" values are 55 to 65%']);axis([ddd ddd+1 50 100]);
        if prtit; print('-djpeg90 ',[way_images_flux, licor_prefix, '_Licor_AGC_' num2str(ddd) '.jpg']); end
    end;
end