read_co2
% VOCALS 2008 :: 2008-10-16 :: Simon de Szoeke
% Does VOCALS have CO2 bulk or covariance flux measurement?
% From GASEXIII: Opens Channel of the UConn Licor 7500
% CO2 and H2O Absorptions are saved for UConn
% Inputs data from ASCII data files

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

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

z=[];
jd_lic=[];jd1=[];
jax1=1;lic1=[];
Licor_CO2=[];
Licor_H2O=[];
Licor_Press=[];
Licor_Temp=[];

for jam=0:23,					    %cycle thru 24 hourly stats files
    if jam<10,
        hr=['0' num2str(jam)];
    else
        hr=num2str(jam);
    end;                            %end if jam
   
    dfl=[way_raw_data_flux, 'day',jd,'\' fname '0', year(3:4) jd hr '_raw.txt']; 
    disp(['Reading Licor file  ',dfl])
    flist=fopen(dfl,'r');
    z=[];temp=[];   
    if flist>0,					    %if the file exists, 
        while feof(flist)==0,
            clear temp;
            if strcmp(fname,'lis')==1
                temp =textscan(flist,'%2f%2f%3f %f %f %f %f %f %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty','Overflow');
            elseif strcmp(fname,'lin')==1
                temp =textscan(flist,'%2f%2f%2f%3f %f %f %f %f %f %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty','Overflow');
            elseif strcmp(fname,'li1')==1 | strcmp(fname,'li2')==1 | strcmp(fname,'li3')==1
                temp =textscan(flist,'%2f%2f%2f%3f %f %f %f %f %f %f %f %*[^\n]','delimiter', ', ', 'headerlines', 1,  'emptyvalue', NaN,'treatAsEmpty','Overflow');
            end;

            try
                z=[z cell2mat(temp)'];
            catch
                [nr1,nl1]=size(temp);
                for ll=1:nl1
                     [l10,c10]=size(temp{1,nl1});
                     if length(temp{1,ll})~=l10
                         temp{1,ll}(l10+1)=[];
                     end;
                 end;
             if isempty(temp{1,1});
                 temp=[];
             end;
                z=[z cell2mat(temp)'];
            end;  %end of catch/try
        end; %end of while feof

        if strcmp(fname,'lis')==1
            ii=find(z(8,:)<50); z(7:8,ii)=z(7:8,ii-1);
            [nr1,nl1]=size(z);
        elseif strcmp(fname,'lin')==1
            ii=find(z(9,:)<50); z(8:9,ii)=z(8:9,ii-1);
            [nr1,nl1]=size(z);
        elseif strcmp(fname,'li1')==1 | strcmp(fname,'li2')==1 | strcmp(fname,'li3')==1
            ii=find(z(11,:)<50); z(10:11,ii)=z(10:11,ii-1);
            [nr1,nl1]=size(z);
        end;        
        
        if nl1<36000
            z=[z ones(11,36000-nl1)*NaN];
            nl1=36000;
        end;
                       
        lic1(:,jax1:jax1+nl1-1)=z;
       
        if strcmp(fname,'lis')==1
            jd1=ddd + (jam+((lic1(1,jax1:jax1+nl1-1)+((lic1(2,jax1:jax1+nl1-1)+lic1(3,jax1:jax1+nl1-1)/1000)/60))/60))/24;
        elseif strcmp(fname,'lin')==1 | strcmp(fname,'li1')==1 | strcmp(fname,'li2')==1 | strcmp(fname,'li3')==1
            jd1=ddd + (lic1(1,jax1:jax1+nl1-1)+((lic1(2,jax1:jax1+nl1-1)+((lic1(3,jax1:jax1+nl1-1)+lic1(4,jax1:jax1+nl1-1)/1000)/60))/60))/24;
        end;
        jd_lic = [jd_lic jd1];
        jax1=jax1+nl1;  
        
        if jam==0 | jam==12 | jam==18

            if strcmp(fname,'lis')==1
                ii=find(~isnan(z(6,:)));
                spechum_lic  =z(6,ii)*18/1000; %g/m3
            elseif strcmp(fname,'lin')==1
                ii=find(~isnan(z(7,:)));
                spechum_lic  =z(7,ii)*18/1000; %g/m3
            elseif strcmp(fname,'li1')==1 | strcmp(fname,'li2')==1 | strcmp(fname,'li3')==1
                ii=find(~isnan(z(9,:)));
                spechum_lic  =z(9,ii)*18/1000; %g/m3
            end;
              
            disp(['Calculating licor humidity spectra for ',hr])
            [S_Lic,F]=psd2(detrend(spechum_lic),length(spechum_lic),flicor);
            [S_Lics,Fsl]=specsmoo(S_Lic,flicor);
            
            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([fname, ' Licor Hygrometer Spectrum:   Day ',jd,'   Hour ',hr]);
            xlabel('Frequency in Hz');
            ylabel('Spectrum * Frequency');
            print('-djpeg90 ',[way_images_flux, fname, '_LicH2O_spectrum' jd hr '.jpg']);
            
            if jam==12,
               specfile=[way_spec_data_flux,fname,'_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;
    end;                            %end if flist
end;

% % 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   

if strcmp(fname,'lis')==1
    Licor_diag =lic1(4,:);
    Licor_CO2  =lic1(5,:);%mmol/m^3
    Licor_H2O  =lic1(6,:);%mmol/m^3
    Licor_Temp =lic1(7,:);%C
    Licor_Press=lic1(8,:);%kPa
elseif strcmp(fname,'lin')==1
    Licor_diag =lic1(5,:);
    Licor_CO2  =lic1(6,:);%mmol/m^3
    Licor_H2O  =lic1(7,:);%mmol/m^3
    Licor_Temp =lic1(8,:);%C
    Licor_Press=lic1(9,:);%kPa
elseif strcmp(fname,'li1')==1 | strcmp(fname,'li2')==1 | strcmp(fname,'li3')==1
    Licor_diag =lic1(5,:);
    Licor_CO2 =lic1(7,:);%mmol/m^3
    Licor_H2O  =lic1(9,:);%mmol/m^3
    Licor_Temp =lic1(10,:);%C
    Licor_Press=lic1(11,:);%kPa
end;
          
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;

ii=find(Licor_Press==0 | Licor_Press>200);jj=find(Licor_Press>0 & Licor_Press<200);
Licor_Press(ii)=median(Licor_Press(jj));
ii=find(Licor_Temp==0 | Licor_Temp>1000);jj=find(Licor_Temp>0 & Licor_Temp<1000);
Licor_Temp(ii)=median(Licor_Temp(jj));

% ii=find(Licor_CO2==0 | Licor_CO2>1000);jj=find(Licor_CO2>0 & Licor_CO2<1000);
% Licor_CO2(ii)=median(Licor_CO2(jj));
ii=find(Licor_CO2>100);Licor_CO2(ii)=NaN;ii=find(Licor_CO2<0);Licor_CO2(ii)=NaN;

% ii=find(Licor_H2O==0 | Licor_H2O>2000);jj=find(Licor_H2O>0 & Licor_H2O<2000);
% Licor_H2O(ii)=median(Licor_H2O(jj));
ii=find(Licor_H2O>500);Licor_H2O(ii)=NaN;ii=find(Licor_H2O<0);Licor_H2O(ii)=NaN;

      
      
Licor_CO2_a=Licor_CO2*8.314*(nanmean(Licor_Temp)+273.15)/nanmean(Licor_Press); %mmol/m3 to umol/mol
Licor_H2O_a=Licor_H2O*18/1000;%mmol/m3 to g/m^3
Rgas=287.1;
% rhoa_lic=nanmean(Licor_Press)*1000/(Rgas*(nanmean(Licor_Temp)+273.16)*(1+0.61*nanmean(Licor_H2O_a)/1250));%air density in licor
%1250 seems to be an approximation of the dry air density.
rhoa_lic=(nanmean(Licor_Press)*1000/(Rgas*(nanmean(Licor_Temp)+273.16)))-(0.61*nanmean(Licor_H2O_a)/1000);%air density in licor


Licor_H2O_a=Licor_H2O_a/rhoa_lic;%g/m^3 to g/kg

if plotit
    figure;plot(jd_lic,Licor_H2O_a);title([fname, ' Licor specific humidity. ', cruise, year]);xlabel(['JD', year]);ylabel('q (g/kg)');axis([ddd ddd+1 -1 30]);
        print('-djpeg90 ',[way_images_flux, fname, '_Licor_spechum_' num2str(ddd) '.jpg']);
    figure;plot(jd_lic,Licor_CO2_a);title([fname, ' Licor CO2. ', cruise, year]);xlabel(['JD', year]);ylabel('CO2 (umol/mol)');axis([ddd ddd+1 0 800])
        print('-djpeg90 ',[way_images_flux, fname, '_Licor_CO2_' num2str(ddd) '.jpg']);
    figure;plot(jd_lic,Licor_Temp);title([fname, ' Licor Temperature [degC]. ', cruise, year]);xlabel(['JD', year]);xlim([ddd ddd+1])
        print('-djpeg90 ',[way_images_flux, fname, '_Licor_temp_' num2str(ddd) '.jpg']);    
    figure;plot(jd_lic,Licor_Press);title([fname, ' Licor Pressure [kPa]. ', cruise, year]);xlabel(['JD', year]);xlim([ddd ddd+1])
        print('-djpeg90 ',[way_images_flux, fname, '_Licor_press_' num2str(ddd) '.jpg']);    
%     figure;plot(jd_lic,Licor_CO2_a);title('Licor CO2 [umol/mol]');
    figure;plot(jd_lic,Licor_H2O_a,jd_pc,qvais,'o');title([fname, ' Licor (blue), Vaisala (circles). ', cruise, year]);xlabel(['JD', year]);ylabel('q (g/kg)');axis([ddd ddd+1 0 50])
        print('-djpeg90 ',[way_images_flux, fname, '_Licor_spechum_comp_' num2str(ddd) '.jpg']);   
    figure;plot(jd_lic,agc,[ddd ddd+1],[55 55],'r',[ddd ddd+1],[65 65],'r'); title([fname, ' AGC value. Typical "clean window" values are 55 to 65%']);axis([ddd ddd+1 50 100]);
        print('-djpeg90 ',[way_images_flux, fname, '_Licor_AGC_' num2str(ddd) '.jpg']);   
end;
