disp('read_lindsey_1')
fclose('all');
clear all;
plotit=input('plot plots? 0=n, 1=y');
%plotit=0;
%prtit=input('append master data file? 0=n, 1=y');
prtit=0;
%yd=date;
g=0.0098;%adiabatic lapse rate
eps_i=0.98;%emissivity of snow
sig_sb=5.67e-8;%Stefan Boltzmann constant
blk=1;%10-min index 
nlines=14298;
stux=zeros(16,nlines);%basic stats file block size, 7 cols 50 lines
%some empirical offsets determined in halfassed manner by CWF
e=['D:\matlabstf\cwf\sheba_data\lindsey_rad\station_2.txt'];
flist=fopen(e,'r');
stux=fscanf(flist,'%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g',[16,nlines]);%read main data block
[aa,bb]=size(stux);%check size
 nlines=bb;        
 %get time stuff
yy=stux(1,:);
mon=stux(2,:);
dd=stux(3,:);
hh=stux(4,:);
mm=stux(5,:);
timex(1,:)=yy;
timex(2,:)=mon;
timex(3,:)=dd;
timex(4,:)=hh;
timex(5,:)=mm;
timex(6,:)=zeros(1,nlines);
jad=365*(yy-97)+juldayv(timex)-1;
Rld=stux(6,:);
Rsd=stux(12,:);
Tcd=stux(10,:);
Tdd=stux(11,:);
Rlt=Rld+sig_sb*(Tcd+273.16).^4+4*sig_sb*((Tcd+273.16).^4-(Tdd+273.16).^4);
sgRld=stux(9,:);
sgRsd=stux(15,:);
chk=stux(16,:);
         
blk=1;  %loop thru 10-min sections to extract and compute 1-hr aves
loop=1;
%zero all accumulators
while loop<nlines
	Rld_ac=0;Rsd_ac=0;Rlt_ac=0;Tcd_ac=0;Tdd_ac=0;sgrl_ac=0;sgrs_ac=0;
   n_ac=0;
   hour=hh(loop);
   jrd(blk)=jad(loop);
   yyx(blk)=yy(loop);
   monx(blk)=mon(loop);
   ddx(blk)=dd(loop);
   mmx(blk)=mm(loop);
   clk=chk(loop);
   	%accumulate for 1-hr aves
   while hh(loop)==hour   	
 	  	Rld_ac=Rld_ac+Rld(loop);
  	  	Rlt_ac=Rlt_ac+Rlt(loop);
      Rsd_ac=Rsd_ac+Rsd(loop);
   	Tcd_ac=Tcd_ac+Tcd(loop);
   	Tdd_ac=Tdd_ac+Tdd(loop);
   	sgrl_ac=sgrl_ac+sgRld(loop)*sgRld(loop);
      sgrs_ac=sgrs_ac+sgRsd(loop)*sgRsd(loop);
      clk=max(clk,chk(loop));
      n_ac=n_ac+1;
      loop=loop+1;
      if loop>nlines,
         break
       end;  
   end;%end while hh
      %compute 1-hr aves
   Rldx(blk)=Rld_ac./n_ac;
   Rsdx(blk)=Rsd_ac./n_ac;
   Tcdx(blk)=Tcd_ac./n_ac;
   Tddx(blk)=Tdd_ac./n_ac;
   sgRlx(blk)=sqrt(sgrl_ac./n_ac);
   sgRsx(blk)=sqrt(sgrs_ac./n_ac);
   Rltx(blk)=Rlt_ac./n_ac;
   chkx(blk)=clk;
   blk=blk+1;
end;%end while blk   
[npz,a]=size(jrd);
%10-min time series
np=blk-1;


%1-hr time series
j_st=floor(jrd(1));
j_en=ceil(jrd(np));

if prtit==1;
  prt_stats_lin;
end;

if plotit==1;
	%  set some plot limits
	Tll=-40;
	Tul=-10;
	Rtul=290;
	Rtll=110;
	Rsul=200;	
	Rsll=-10;

%plot T profiles
	figure;
	plot(jrd(1:np),Tcdx(1:np),jrd(1:np),Tddx(1:np));
	xlabel('Jday (1997)');
	ylabel('T_{cs}, T_{dm} (C)');
	axis([j_st j_en Tll Tul]);
	Title('case=blu; dome=grn'); 
	
	figure;
	plot(jrd(1:np),Rltx(1:np),jrd(1:np),chkx(1:np)*100);
	xlabel('Jday (1997)');
	ylabel('R_{ld} (W/m^2)');
	axis([j_st j_en Rtll Rtul]);
	Title('Pyrgeomete=blu'); 
	 
end; 
