disp('read_pc_day_rhb_07')
fclose('all');
% modified on 30-mar-01 by D.E.Lane for Kai 2001 cruise
% modified 10 nov 03 by cfairall for the weller 2003 cruise
% modified 09 jul 04 by hare
% modified April 06 by lbariteau for the AMMA 2006 cruise (changement of
% the variables in the Mean Campbell datalogger)

jd=num2str(ddd);
if ddd<100
   jd=['0' num2str(ddd)];
end;

if ddd<10
   jd=['00' num2str(ddd)];
end;

plotit=0;
prtit=1;

clear jad;
hr='00';							%input('gmt hr','s');
g=0.0098;						    %adiabatic lapse rate
eps_w=0.97;						    %emissivity of water
sig_sb=5.67e-8;				        %Stefan Boltzmann constant
blk=1;							    %10-min index 
jax=1;
backchk=0;
clear st2
total_num_fields=72;

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
   
   	e=[way_raw_data, 'day', jd,'\P2_' jd hr];
    disp(['Reading means file from hour ',num2str(jam)]);
    flist=fopen(e,'r');
   
    if flist>0,					    %if the file exists, 
              
        dmp=fgets(flist);		    %read header line
		
        while feof(flist)==0,
            if backchk==0,
            tzilch=jam;
            backchk=1;
            end;						%end if
            clear st1 stz;
            if isfinite(flist),
                stx=fscanf(flist,'%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g',[77,1]);
                %                          1                           2                          3                          4                         5
                [nr1,nl1]=size(stx);
                if nr1==77
                    st2(:,jax:jax+nl1-1)=stx;
      	            jax=jax+nl1;
	            else
   		            x=fgets(flist);clear x
                end;                %end nr1
            dmp4=fgets(flist);  %read time from blank line
            end;                    %end if st0
        end;                        %end while feof
    end;                            %end if flist
end;                                %end for jam

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

press=st2(34,:);
psp1=st2(28,:);
psp2=st2(29,:);
org=st2(30,:);
Tvais=st2(24,:);
ii=find(Tvais<-100);Tvais(ii)=-999;    %ii=logical(Tvais<-100);Tvais(ii)=-999;
Rhvais=st2(25,:);

avgpsp=mean(psp1)
avgrain=mean(org);
avgTvais=mean(Tvais)
avgRHvais=mean(Rhvais);avgQvais=avgRHvais/100*qsea(avgTvais)
pp=ones(length(Tvais),1)*1009;
qvais=qair_p([Tvais' Rhvais' pp]);

Tsea=st2(23,:)+.01;   %correct for lead resistance 20C

Tc1=st2(19,:);%-0.35;  %correct for bias, day 321
Td1=st2(20,:);%-0.3;
Tc2=st2(21,:);%-.1;    %correct for bias
%correct  for bad connection days 221, 222
ik=find(Tc2<-10);Tc2(ik)=Tc1(ik);
%***********
Td2=st2(22,:);%+.35;%oct 10 2005
k=find(Td2<-10);Td2(k)=Tc2(k)+Td1(k)-Tc1(k);       %NOTE PROBLEM WITH Tdome on unit 2

%Td2=Td1;        %NOTE PROBLEM WITH Tdome on unit 2
therm1=st2(26,:);
therm2=st2(27,:);

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);
pp=ones(length(Tvais),1)*1009;

for iii=2:length(org);
    if org(iii)==0,
        org(iii)=org(iii-1);
    end;
    if org_carrier(iii)==0,
        org_carrier(iii)=org_carrier(iii-1);
    end;
    if Tc1(iii)==0,
        Tc1(iii)=Tc1(iii-1);
    end;
    if Tc2(iii)==0,
        Tc2(iii)=Tc2(iii-1);
    end;
    if Td1(iii)==0,
        Td1(iii)=Td1(iii-1);
    end;
    if Td2(iii)==0,
        Td2(iii)=Td2(iii-1);
    end;
    if pir1(iii)==0,
        pir1(iii)=pir1(iii-1);
    end;
    if pir2(iii)==0,
        pir2(iii)=pir2(iii-1);
    end;
    if psp1(iii)==0,
        psp1(iii)=psp1(iii-1);
    end;
    if psp2(iii)==0,
        psp2(iii)=psp2(iii-1);
    end;
    if Tvais(iii)==0,
        Tvais(iii)=Tvais(iii-1);
    end;
    if Rhvais(iii)==0,
        Rhvais(iii)=Rhvais(iii-1);
    end;
    if Tsea(iii)==0,
        Tsea(iii)=Tsea(iii-1);
    end;
end;

pspm=psp1;%(psp1+psp2)/2;Note psp2 had dome condensation most of tao04

%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;
ii=find(isnan(pirm));
pirm(ii)=pir1(ii);
ii=find(isnan(pirm));pirm(ii)=-999;

avgsnake=mean(Tsea)

% 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));
jdn = 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))
jdn(1) = jdn(1)-1;
end    

%jad=str2num(num2str(jd))+((st2(1,:)-st2(1,1))/1000/60/60 +tzilch)/24;
jad = jdn;

% % DEW Oct 14, 2006
% if 0
% figure;plot(psp1,psp2,'.');title('ETL PSPs Solar FLux (W/m^2)');xlabel('PSP1');ylabel('PSP2');
% axis('square')
% axis([min(min(psp1),min(psp2)),max(max(psp1),max(psp2)),min(min(psp1),min(psp2)),max(max(psp1),max(psp2))])
% figure;plot(pir1,pir2,'.');title('ETL PIRs IR Flux (W/m^2)');xlabel('PIR1');ylabel('PIR2');
% axis('square')
% axis([min(min(pir1),min(pir2)),max(max(pir1),max(pir2)),min(min(pir1),min(pir2)),max(max(pir1),max(pir2))])
% corrcoef(psp1,psp2),corrcoef(pir1,pir2)
% mean(psp1),mean(psp2),mean(pir1),mean(pir2)
% 
% figure;plot(Tc1,Tc2,'.');title('ETL PSPs Case temp (C)');xlabel('PSP1');ylabel('PSP2');
% axis('square')
% axis([min(min(Tc1),min(Tc2)),max(max(Tc1),max(Tc2)),min(min(Tc1),min(Tc2)),max(max(Tc1),max(Tc2))])
% figure;plot(Td1,Td2,'.');title('ETL PSPs Dome temp (C)');xlabel('PIR1');ylabel('PIR2');
% axis('square')
% axis([min(min(Td1),min(Td2)),max(max(Td1),max(Td2)),min(min(Td1),min(Td2)),max(max(Td1),max(Td2))])
% mean(Tc1),mean(Tc2),mean(Td1),mean(Td2)
% corrcoef(Tc1,Tc2),corrcoef(Td1,Td2)
% end

figure;plot(jad,press,'r');title('Vaisala Pressure (case 1 (red),2(magenta); dome 1(black),);xlabel(['DOY', year]);ylabel('Pressure (mb)');
figure;plot(jad,Tc1,'r',jad,Tc2,'m',jad,Td1,'k',jad,Td2,'c',jad,Tsea,'b',jad,Tvais,'g');title('ETL Temperatures (case 1 (red),2(magenta); dome 1(black),2(cyan); snake(blue), Tair(green))');xlabel(['DOY', year]);ylabel('Temperature (C)');
figure;plot(jad,psp1,jad,psp2);title('ETL psp');xlabel(['DOY', year]);ylabel('Solar FLux (W/m^2)');legend('PSP1','PSP2')
figure;plot(jad,pir1,jad,pir2);title('ETL pir');xlabel(['DOY', year]);ylabel('IR Flux (W/m^2)');legend('PIR1','PIR2')
figure;plot(jad,Rhvais);title('ETL RH');xlabel(['DOY', year]);ylabel('Relative Humidity (%)');
figure;plot(jad,org);title('ETL org rain rate');xlabel(['DOY', year]);ylabel('Rain Rate (mm/hr)');
figure;plot(jad,aspir_on);title('ETL aspirator');xlabel(['DOY', year]);ylabel('Backflow Indicator (Volt). If 0V, there is backflow');%axis([max(jad)-1 max(jad) -0.2 5.2])
figure;plot(jad,org_carrier);title('ETL org carrier');xlabel(['DOY', year]);ylabel('Rain Gauge Function (V)');
figure;plot(jad,Tsea);title('PSD seasnake');xlabel(['DOY', year]);ylabel('Temperature (degC)');
save d:\data\ETLrad jad psp1 psp2 pir1 pir2  % save data for comparison with SCS using rad_stat_06.m

np=jax-1;


if prtit==1;
  prt_jas_means_rhb_06;
  prt_jas_rads_rhb_06;
end;

if plotit==1;
   Tll=-25;
   hll=70;
   uul=10;
   vul=10;
   Tul=-0;
   hul=100;
   sul=12;
   Rlul=30;
   Rlll=-100;
   Rtul=300;
   Rtll=125;
   Rsul=1100;
   Rsll=0;

   j_st=jad(1);
   j_en=jad(np);
   figure;plot(jad(1:np),Tvais(1:np),'b',jad(1:np),Tsea(1:np),'g',jad(1:np),Td1(1:np),'r',jad(1:np),Tc1(1:np),'m');xlabel(['Jday', year]);ylabel('Temperature (C)');title('Tair_v=blue; Tsea_s_n_a_k_e=green;  T_d_o_m_e_1=red;  T_c_a_s_e_1=magenta;'); 
   figure;h=plot(jad(1:np),Tvais(1:np),jad(1:np),Tsea(1:np));set(gca,'fontsize',14);set(h,'LineWidth',2);xlabel(['Jday', year]);ylabel('Temperature (C)');title('Tair_v=blue;     Tsea_s_n_a_k_e=green'); 
   figure;h=plot(jad(1:np),pir1(1:np));set(gca,'fontsize',14);set(h,'LineWidth',2);xlabel(['Jday', year]);ylabel('Pir (W/m^2)');title('Downward IR Flux'); 
   figure;h=plot(jad(1:np),psp1(1:np));set(gca,'fontsize',14);set(h,'LineWidth',2);xlabel(['Jday', year]);ylabel('Psp (W/m^2)');axis([j_st j_en Rsll Rsul]);title('Downward Solar Flux'); 

   hold on;
   figure;plot(jad(1:np),org(1:np));xlabel(['Jday', year]);ylabel('ORG Precipitation Rate (mm/hr)');axis([j_st j_en 0 50]);title('STI Optical Rain Gauge');
end; 

%clear jd jax jam jad j_st j_en hul hr hll g flist f eps_w e dmp blk backchk ans
%clear vul uul tzilch therm1 therm2 sul stz sty stx st2 st1 st0 sig_sb psp1 psp2 prtit prt_it plotit
%clear Tc1 Tc2 Td1 Td2 Tll Tsea Tul Tvais jazz m0 n0 nl1 np npz nr1 org pir1 pir2
%clear Rhvais Rlll Rlul Rsll Rsul Rtll Rtul