disp('read_stats_day_mod7_mult')
fclose('all');
clear all;
%jd=input('julian day, ddd','s');
jd1=115;
jd2=151;
screwup=0;
%screwup=1;   %if starting between 8/31  19:51:46  and 9/3 18:26:22
%prtit=input('append master data file? 0=n, 1=y');
prtit=1;
%yd=date;
%jif=datevec(yd);
%xx=[jif(1) jif(2) jif(3) 0 0 0];
%jd=num2str(julday(xx)-1);
g=0.0098;%adiabatic lapse rate
eps_i=0.99;%emissivity of snow
sig_sb=5.67e-8;%Stefan Boltzmann constant

%radiometer calibration coefficients
% the factor 860/894 was inserted on 990318 because the SW wires has been reversed
rsdc1=(894/880)*(860/894);rsdc2=(894/880)*(860/894);
rsuc1=(860/856)*(894/860);rsuc2=(860/856)*(894/860);
rldc1=369/365;rldc2=369/367.5;
rluc1=421/398;rluc2=421/416;

bd1=4.0;bd2=3.45;
bu1=4.5;bu2=3.7;

%radiometer thermistor corrections 
%tcdc=0.2;%summer
tcdc=-0.5;%winter

def_t=999;def_u=999;def_rad=999;def_rr=999;def_rh=999;def_dir=999;
def_z=[999 999 999 999 999];

mnp10=2;
smpp10=mnp10;
smpp10s=mnp10*600;
perphr=1;
 
stux=zeros(7,50);%basic stats file block size, 7 cols 50 lines
%some empirical offsets determined in halfassed manner by CWF
%Tcor=[+.05 -.05 -.05 +.05 +.05];
ucor=[+.04 -.1 -.05  +.00 +.3];
vcor=[-.2 .00 .00 +.00 -.3];
%some empirical offsets determined by taking the mean of the calibrations
% on JD105 and JD143 relative to level 1, and assuming that the absolute
% calibration of level 1 as given by CWF above is correct -- OP
%Tcor=[+.05 -.05 -.13 +.01 -.14];
Tcor=[+.00 -.00 -.0 +.00 -.0];
%T_ice_corr=-0.5;
%T_snw_corr=-0.5;
%T_GE_corr=+0.5;
T_ice_corr=0;
T_snw_corr=0;
T_GE_corr=0;

%****
% day loop
dytot=272+365-302;
for jdy=jd1:jd2
blk=1;%10-min index
%calculate day of experiment
if jdy < 290
  jdex=jdy+365-302;
else
  jdex=jdy-302;
end;
%calculate interpolated radiometer coefficients for this day
rsdc=(1-(jdex/dytot))*rsdc1+(jdex/dytot)*rsdc2;
rsuc=(1-(jdex/dytot))*rsuc1+(jdex/dytot)*rsuc2;
rldc=(1-(jdex/dytot))*rldc1+(jdex/dytot)*rldc2;
rluc=(1-(jdex/dytot))*rluc1+(jdex/dytot)*rluc2;
bd=(1-(jdex/dytot))*bd1+(jdex/dytot)*bd2;
bu=(1-(jdex/dytot))*bu1+(jdex/dytot)*bu2;

stux=zeros(7,50);%basic stats file block size, 7 cols 50 lines 
ihrfl=0;%reset flag for if data found on this day
%jd=num2str(jdy);
if jdy<10
      jd=['00'num2str(jdy)]
   elseif jdy<100,
      jd=['0'num2str(jdy)]
   else
      jd=num2str(jdy)
   end;
clear jaz
%determine magnetic orientation of tower
if jdy>370,
  or_mag=40;
  decl=25;
else
  or_mag=330;
  decl=30;
end
Tpl=zeros(1,6);
zlch=zeros(1,6);
hr='00';%input('gmt hr','s');
nhr=0;
for jam=0:23,%cycle thru 24 hourly stats files
   if jam<10,
      hr=['0'num2str(jam)];
   else
      hr=num2str(jam);
   end;
   
%Find the data 
%e=['C:\SHEBA\data\98' jd hr '.sta'];
%e=['F:\users\data\98' jd hr '.sta'];
%e=['E:\ascii\98' jd hr '.sta'];
%e=['C:\SHEBA\data\sta_zip_data\97' jd hr '.sta'];
%e=['D:\DasWin32\data\ascii\98' jd hr '.sta'];
e=['C:\Olas_Files\SHEBA\sta_zip_data\98' jd hr '.sta'];
%e=['G:\Data\98' jd hr '.sta'];
   flist=fopen(e,'r');
   if flist>0,   %if the file exists,
      nhr=nhr+1; 
      %find instrument heights
      for jj=1:25,
      	bum=fgets(flist);
      	if jj==9,
            bam=bum;
            if jdy>304|jdy<285;
               z(blk,1)=str2num(bum(7:10));
            else   
           	z(blk,1)=str2num(bum(7:7));
            end;%if jd
         end;%if jj
      end; %end jj
     z(blk,2)=3.38;%fixed heights levels 2-4
     z(blk,3)=5.33;
     z(blk,4)=9.12;
      if jdy>304 & jdy < 319;
         z(blk,5)=str2num(bum(7:11));
      elseif jdy>=319 | jdy==1
         z(blk,5)=str2num(bum(7:10));
      elseif jdy<285;
         z(blk,5)=str2num(bum(7:11));
      elseif jdy>=285 & jdy<=304
         z(blk,5)=str2num(bum(7:8));
      end;%if jd 
      blks=blk;
      fclose(flist);%close and re-open file to read full data

      flist=fopen(e,'r');
      
      loop=1;  %loop thru 10-min sections to extract and compute 1-hr aves
      while feof(flist)==0,%go until you run out of 10-min sections
         dmp=fgets(flist);%read header lines and skip blanks
         dmp1=fgets(flist);
   	   dmp2=fgets(flist);
         dmpx=fgets(flist);
   	   dat(blk,1:6)=datevec([dmpx(48:55) ' ' dmpx(57:64)]);
         if eval(dmpx(57:58))==12,
            dat(blk,4)=12;
         end;
         jad(blk)=str2num(jd)+dat(blk,4)/24+dat(blk,5)/24/60;
         jz(blk)=jad(blk);
   		dmp3=fgets(flist);
		dmp4=fgets(flist);
		dmp5=fgets(flist);
		dmp6=fgets(flist);
if dat(blk,2) == 8 & dat(blk,3) == 31 & dat(blk,4) == 19 & dat(blk,5) == 51 & dat(blk,6) == 46 ;
  screwup = 1;
end;
if dat(blk,2) == 9 & dat(blk,3) == 3 & dat(blk,4) == 18 & dat(blk,5) ==  26 & dat(blk,6) == 22 ;
  screwup = 0;
end; 
if screwup == 1;
  stux1=fscanf(flist,'%*s%*s%g%g%g%g%g%g%g',[7,16]);%read first part
  stux2=fscanf(flist,'%*s%*s%*s%g%g%g%g%g%g%g',[7,4]);%read messed up part
  stux3=fscanf(flist,'%*s%*s%g%g%g%g%g%g%g',[7,30]);%read last part
  stux= [stux1 stux2 stux3];
else
  stux=fscanf(flist,'%*s%*s%g%g%g%g%g%g%g',[7,50]);%read main datak
end
[aa,bb]=size(stux);%check size
      %zero all accumulators
      th_ac=zeros(1,5);u_ac=zeros(1,5);v_ac=zeros(1,5);w_ac=zeros(1,5);
      Rlu_ac=0;Rld_ac=0;Rsu_ac=0;Rsd_ac=0;
      Tcu_ac=0;Tcd_ac=0;Tdd_ac=0;Tdu_ac=0;
      Rhv_ac=zeros(1,5);Tice_ac=0;Tsnw_ac=0;Vbr_ac=0;org_ac=0;
      tn_ac=0;tdn_ac=0;ncr_ac=0;ts_ac=zeros(1,5);
      ntv_ac=zeros(1,5);nu_ac=zeros(1,5);nv_ac=zeros(1,5);nw_ac=zeros(1,5);
      nrh_ac=zeros(1,5);nts_ac=zeros(1,5); 
      ntan_ac=0; ntdan_ac=0; nrld_ac=0;nrlu_ac=0;nrsd_ac=0;nrsu_ac=0;ntcd_ac=0;ntcu_ac=0;
      ntdd_ac=0;ntdu_ac=0;ntice_ac=0;ntsn_ac=0;nor_ac=0;nncr_ac=0;nvbr_ac=0;
         if bb>49, %only compute on full data blocks
            for i=1:5; %five levels of instruments
% sonics
               u_cnt(blk,i)=stux(1,1+4*(i-1));
               v_cnt(blk,i)=stux(1,2+4*(i-1));
               w_cnt(blk,i)=stux(1,3+4*(i-1));
               ts_cnt(blk,i)=stux(1,4+4*(i-1));
                   uraw(blk,i)=stux(2,1+4*(i-1));
		   vraw(blk,i)=stux(2,2+4*(i-1));
		   wraw(blk,i)=stux(2,3+4*(i-1));
		   tsraw(blk,i)=stux(2,4+4*(i-1));

		   sgu(blk,i)=stux(3,1+4*(i-1));
		   sgv(blk,i)=stux(3,2+4*(i-1));
		   sgw(blk,i)=stux(3,3+4*(i-1));
		   sgT(blk,i)=stux(3,4+4*(i-1));
%Vaisalas
		   Tv_cnt(blk,i)=stux(1,26+1*(i-1));
		   Rh_cnt(blk,i)=stux(2,31+1*(i-1));
		   Tvs(blk,i)=stux(2,26+1*(i-1));
		   Rhv(blk,i)=stux(2,31+1*(i-1));
	
		   sgTv(blk,i)=stux(3,26+1*(i-1));
		   sgRh(blk,i)=stux(3,31+1*(i-1));
         end;%end for i
         Tan_cnt(blk)=stux(1,25);
	      Tan(blk)=stux(2,25)+T_GE_corr;%Andreas General Eastern T measure, empirical offset
         sgTan(blk)=stux(3,25);
         Tdan_cnt(blk)=stux(1,36);
         Tdan(blk)=stux(2,36);%Andreas General Eastern Td measure
         sgTdan(blk)=stux(3,36);

         Tsn_cnt(blk)=stux(1,23);
   	   Tsnow(blk)=stux(2,23)+T_snw_corr;%Tsnow, ice, and PIR empirical offsets .5 C
         sgTsn(blk)=stux(3,23);
         Tice_cnt(blk)=stux(1,24);
         Tice(blk)=stux(2,24)+T_ice_corr; %Tice at snow ice interface
         sgTice(blk)=stux(3,24);
         or_cnt(blk)=stux(1,39);
         org(blk)=stux(2,39);%STI optical rain gauge, mm/hr
         ncr_cnt(blk)=stux(1,40);
         ncr(blk)=stux(2,40);%NCAR weighing rain gauge, mm/min
         rlu_cnt(blk)=stux(1,45);
	   Rlu(blk)=stux(2,45);
         rld_cnt(blk)=stux(1,46);
         Rld(blk)=stux(2,46);
         Rlu(blk)=Rlu(blk)*rluc;
         Rld(blk)=Rld(blk)*rldc;
         vbr_cnt(blk)=stux(1,47);
	   Vbrn(blk)=stux(2,47);%CRREL Barnes Tir
         sgvbrn(blk)=stux(3,47);
         rsd_cnt(blk)=stux(1,48);
	   Rsd(blk)=stux(2,48);
         rsu_cnt(blk)=stux(1,49);
	   Rsu(blk)=stux(2,49);
         Rsd(blk)=Rsd(blk)*rsdc;
         Rsu(blk)=Rsu(blk)*rsuc;

         tcu_cnt(blk)=stux(1,41);
         tcd_cnt(blk)=stux(1,42);
         tdu_cnt(blk)=stux(1,43);
         tdd_cnt(blk)=stux(1,44);
	   Tcu(blk)=stux(2,41)+tcdc;
	   Tcd(blk)=stux(2,42)+tcdc;
	   Tdu(blk)=stux(2,43)+tcdc;
	   Tdd(blk)=stux(2,44)+tcdc;         
	   sgRlu(blk)=stux(3,45);
	   sgRld(blk)=stux(3,46);
	   sgRsu(blk)=stux(3,48);
	   sgRsd(blk)=stux(3,49);
         
	   sgTcu(blk)=stux(3,41);
	   sgTcd(blk)=stux(3,42);
	   sgTdu(blk)=stux(3,43);
           sgTdd(blk)=stux(3,44);
         	%accumulate for 1-hr aves, eliminate missing or "odd" 10-min periods
           for j=1:5
             if Tv_cnt(blk,j) >= smpp10 & Tvs(blk,j) < 20 & Tvs(blk,j)>=-60 & sgTv(blk,j)<4
                th_ac(j)=th_ac(j)+Tvs(blk,j);
                ntv_ac(j)=ntv_ac(j)+1;
             end;
             if u_cnt(blk,j) >= smpp10s & uraw(blk,j) < 30 & uraw(blk,j)>-30 & sgu(blk,j)<6
         	u_ac(j)=u_ac(j)+uraw(blk,j);
                nu_ac(j)=nu_ac(j)+1;
             end;
             if v_cnt(blk,j) >= smpp10s & vraw(blk,j) < 30 & vraw(blk,j)>-30 & sgv(blk,j)<6
         	v_ac(j)=v_ac(j)+vraw(blk,j);
                nv_ac(j)=nv_ac(j)+1;
             end;
             if w_cnt(blk,j) >= smpp10s & wraw(blk,j) < 5 & wraw(blk,j)>-5 & sgw(blk,j)<2
         	w_ac(j)=w_ac(j)+wraw(blk,j);
                nw_ac(j)=nw_ac(j)+1;
             end;
             if ts_cnt(blk,j) >= smpp10s & tsraw(blk,j) < 20 & tsraw(blk,j)>=-60 & sgT(blk,j)<4
         	ts_ac(j)=ts_ac(j)+tsraw(blk,j);
                nts_ac(j)=nts_ac(j)+1;
             end;
             if Rh_cnt(blk,j) >= smpp10 & Rhv(blk,j) < 110 & Rhv(blk,j)>0 & sgRh(blk,j)<10
         	Rhv_ac(j)=Rhv_ac(j)+Rhv(blk,j);
                nrh_ac(j)=nrh_ac(j)+1;
             end;
             end;%for j=1:5
             if Tice_cnt(blk) >= smpp10 & Tice(blk) > -60 & Tice(blk)<20 & sgTice(blk)<4
         	Tice_ac=Tice_ac+Tice(blk);
                ntice_ac=ntice_ac+1;
             end;
             if Tsn_cnt(blk) >= smpp10 & Tsnow(blk) > -60 & Tsnow(blk)<20 & sgTsn(blk)<4
                Tsnw_ac=Tsnw_ac+Tsnow(blk);
                ntsn_ac=ntsn_ac+1;
             end;
             if Tan_cnt(blk) >= smpp10 & Tan(blk) > -60 & Tan(blk)<20 & sgTan(blk)<4
         	    tn_ac=tn_ac+Tan(blk);
                ntan_ac=ntan_ac+1;
             end;
             if Tdan_cnt(blk) >= smpp10 & Tdan(blk) > -60 & Tdan(blk)<20 & sgTdan(blk)<4
         	    tdn_ac=tdn_ac+Tdan(blk);
                ntdan_ac=ntdan_ac+1;
             end;
             if or_cnt(blk) >= smpp10 & org(blk) < 150 & org(blk) >= 0
         	org_ac=org_ac+org(blk);
                nor_ac=nor_ac+1;
             end;
             if ncr_cnt(blk) >= smpp10 & ncr(blk) < 150 & ncr(blk) >= 0
         	ncr_ac=ncr_ac+ncr(blk);
                nncr_ac=nncr_ac+1;
             end;
             if rld_cnt(blk) >= smpp10 & Rld(blk) < 50 & Rld(blk) > -150 & sgRld(blk)<30
        	Rld_ac=Rld_ac+Rld(blk);
                nrld_ac=nrld_ac+1;
             end;
             if rlu_cnt(blk) >= smpp10 & Rlu(blk) < 50 & Rlu(blk) > -100 & sgRlu(blk)<30               
         	Rlu_ac=Rlu_ac+Rlu(blk);
                nrlu_ac=nrlu_ac+1;
             end;
             if rsd_cnt(blk) >= smpp10 & Rsd(blk) < 900 & Rsd(blk) > -10 & sgRsd(blk)<150
         	Rsd_ac=Rsd_ac+Rsd(blk);
                nrsd_ac=nrsd_ac+1;
             end;
             if rsu_cnt(blk) >= smpp10 & Rsu(blk) < 900 & Rsu(blk) > -10 & sgRsu(blk)<150
         	Rsu_ac=Rsu_ac+Rsu(blk);
                nrsu_ac=nrsu_ac+1;
             end;
             if tcd_cnt(blk) >= smpp10 & Tcd(blk) > -60 & Tcd(blk)<20 & sgTcd(blk)<4
         	Tcd_ac=Tcd_ac+Tcd(blk);
                ntcd_ac=ntcd_ac+1;
             end;
             if tcu_cnt(blk) >= smpp10 & Tcu(blk) > -60 & Tcu(blk)<20 & sgTcu(blk)<4
         	Tcu_ac=Tcu_ac+Tcu(blk);
                ntcu_ac=ntcu_ac+1;
             end;
             if tdd_cnt(blk) >= smpp10 & Tdd(blk) > -60 & Tdd(blk)<20 & sgTdd(blk)<4
         	Tdd_ac=Tdd_ac+Tdd(blk);
                ntdd_ac=ntdd_ac+1;
             end;
             if tdu_cnt(blk) >= smpp10 & Tdu(blk) > -60 &Tdu(blk)<10 & sgTdu(blk)<4 & abs(Tdu(blk)-Tcu(blk))<1.5
         	Tdu_ac=Tdu_ac+Tdu(blk);
                ntdu_ac=ntdu_ac+1;
             end;
             if vbr_cnt(blk) >= smpp10 & Vbrn(blk) < 1.5 & Vbrn(blk) > 0 & sgvbrn(blk)<0.1
         	Vbr_ac=Vbr_ac+Vbrn(blk);
                nvbr_ac=nvbr_ac+1;
             end;
   	     loop=loop+1;
      	     blk=blk+1;
             strk=fgets(flist);
         end;%end if bb

      %compute aves
      %Vaisalas
      i10=blk-1;
      for j=1:5
      if ntv_ac(j) >= perphr
         hlx=th_ac(j)/ntv_ac(j);
         Th(i10,j)=hlx+Tcor(j);%correction for calibrations
      else
        Th(i10,j)=def_t;
      end;
      if nrh_ac(j) >= perphr
        Rhvx(i10,j)=Rhv_ac(j)/nrh_ac(j);
      else
        Rhvx(i10,j)=def_rh;
      end;
%sonics
      if nu_ac(j)>= perphr
        ulx(i10,j)=u_ac(j)/nu_ac(j)+ucor(j);
      else
        ulx(i10,j)=def_u;
      end;
      if nv_ac(j)>= perphr
        vlx(i10,j)=v_ac(j)/nv_ac(j)+vcor(j);
      else
        vlx(i10,j)=def_u;
      end;
      if nw_ac(j)>= perphr
        wlx(i10,j)=w_ac(j)/nw_ac(j);
      else
        wlx(i10,j)=def_u;
      end;
      if nts_ac(j)>= perphr
        tsx(i10,j)=ts_ac(j)/nts_ac(j);
      else
        tsx(i10,j)=def_t;
      end;
      end;%for j=1:5
%other instruments
      if ntan_ac >= perphr
        Th_tn(i10)=tn_ac/ntan_ac;%Andreas' GE
      else
        Th_tn(i10)=def_t;
      end;
      if ntdan_ac >= perphr
        Tdh_tn(i10)=tdn_ac/ntdan_ac;%Andreas' GE
      else
        Tdh_tn(i10)=def_t;
      end; 
      if nrld_ac >= perphr
        Rldx(i10)=Rld_ac/nrld_ac;
      else
        Rldx(i10)=def_rad;
      end;
      if nrlu_ac >= perphr
        Rlux(i10)=Rlu_ac/nrlu_ac;
      else
        Rlux(i10)=def_rad;
      end;
      if nrsd_ac >= perphr
        Rsdx(i10)=Rsd_ac/nrsd_ac;
      else
        Rsdx(i10)=def_rad;
      end;
      if nrsu_ac >= perphr
        Rsux(i10)=Rsu_ac/nrsu_ac;
      else
        Rsux(i10)=def_rad;
      end;
      if ntcd_ac >= perphr
        Tcdx(i10)=Tcd_ac/ntcd_ac;
      else
        Tcdx(i10)=def_t;
      end;
      if ntcu_ac >= perphr
        Tcux(i10)=Tcu_ac/ntcu_ac;
      else
        Tcux(i10)=def_t;
      end;
      if ntdd_ac >= perphr
        Tddx(i10)=Tdd_ac/ntdd_ac;
      else
        Tddx(i10)=def_t;
      end;
      if ntdu_ac >= perphr
        Tdux(i10)=Tdu_ac/ntdu_ac;
      else
        Tdux(i10)=def_t;
      end;
      if ntice_ac >= perphr
        Ticex(i10)=Tice_ac/ntice_ac;
      else
        Ticex(i10)=def_t;
      end;
      if ntsn_ac >= perphr
        Tsnwx(i10)=Tsnw_ac/ntsn_ac;
      else
        Tsnwx(i10)=def_t;
      end;
      if nor_ac >= perphr
        orgx(i10)=org_ac/nor_ac;
      else
        orgx(i10)=def_rr;
      end;
      if nncr_ac >= perphr
        ncrx(i10)=ncr_ac/nncr_ac*60;%convert to mm/hr
      else
        ncrx(i10)=def_rr;
      end;
      if nvbr_ac >= perphr
        Vbrx(i10)=Vbr_ac/nvbr_ac;
      else
        Vbrx(i10)=def_t;
     end; 
  end;%end while feof
     for j=1:5
        z(blks+1:blk-1,j)=z(blks,j);
     end; %for j
     ihrfl=1;%at least one file of data found for this day
   else;
%no file found for this hour
%Vaisalas
for ib=blk:blk+5
   jz(ib)=str2num(jd)+(jam+(ib-blk)*10/60)/24;
   z(ib,:)=def_z;
      for j=1:5
        Th(ib,j)=def_t;
        Rhvx(ib,j)=def_rh;
%sonics
        ulx(ib,j)=def_u;
        vlx(ib,j)=def_u;
        wlx(ib,j)=def_u;
        tsx(ib,j)=def_t;
      end;%for j=1:5
%other instruments
        Th_tn(ib)=def_t;
        Tdh_tn(ib)=def_t;
        Rldx(ib)=def_rad;
        Rlux(ib)=def_rad;
        Rsdx(ib)=def_rad;
        Rsux(ib)=def_rad;
        Tcdx(ib)=def_t;
        Tcux(ib)=def_t;
        Tddx(ib)=def_t;
        Tdux(ib)=def_t;
        Ticex(ib)=def_t;
        Tsnwx(ib)=def_t;
        orgx(ib)=def_rr;
        ncrx(ib)=def_rr;
        Vbrx(ib)=def_t;
      end;%for ib
      blk=blk+6;
   end;%end if flist   
end;%end jam
%if ihrfl==1;
[npz,a]=size(z);
%npz=blk;
jz(1);

for j=1:npz
if Rldx(j)<def_rad & Tcdx(j)<def_t & Tddx(j)<def_t;
  Rtdx(j)=Rldx(j)+sig_sb*((Tcdx(j)+273.15)^4)+bd*sig_sb*((Tcdx(j)+273.15)^4-(Tddx(j)+273.15)^4);
else
  Rtdx(j)=def_rad;
end
if Rlux(j)<def_rad & Tcux(j)<def_t & Tdux(j)<def_t;
  Rtux(j)=Rlux(j)+sig_sb*((Tcux(j)+273.15)^4)+bu*sig_sb*((Tcux(j)+273.15)^4-(Tdux(j)+273.15)^4);
else
  Rtux(j)=def_rad;
end
if Rtux(j)<def_rad & Rtdx(j)<def_rad;
  Tsrx(j)=((Rtux(j)-(1-eps_i)*Rtdx(j))/sig_sb/eps_i)^0.25-273.15;
else
  Tsrx(j)=def_t;
end;
%esx=6.112*exp(Th.*17.502./(Th+241)).*Rhvx/100;
%Qx=0.622*esx./(1020-0.378.*esx)*1000;% water vapor mixing ratio
%rhox=102000./((Th+273.15).*(1+0.61*qx/1000)*287);% air density, kg/m^3
%Th=Th+g*z;%potential temperature
%Th_tn=Th_tn+g*.12;
%Barnes calculations
if Vbrx(j) < def_t
  Nrx=20.265*Vbrx(j)+10.453;
  Tbrx(j)=-96.888+5.5859*Nrx-0.11727*Nrx.*Nrx+1.2436e-3*Nrx.*Nrx.*Nrx;%Barnes surface T in Cent, uncorrected for emissivity/reflection
  Qdnx=min(Rtux(j),2*Rtdx(j)-Rtux(j));%empirical estimate of Barnes downward based on eppley downward
  Rbux=sig_sb*(Tbrx(j)+273.15).^4;
  Tsbx(j)=((Rbux-(1-eps_i)*Qdnx)/sig_sb/eps_i).^0.25-273.15;
else;
  Tbrx(j)=def_t;
  Tsbx(j)=def_t;
end;
end;%for j=1:npz

%calculate winds
for j=1:npz;
  for i=1:5 
   if ulx(j,i)<def_u & vlx(j,i)<def_u & wlx(j,i)<def_u
     um=ulx(j,i);
     vm=vlx(j,i);
     wm=wlx(j,i);
   %rotate to streamwise flow
     s=sqrt(um^2+vm^2);
     thet=atan2(vm,um);
     phi=atan2(wm,s);
   
     rot(1,1,i)=cos(phi)*cos(thet);
     rot(1,2,i)=cos(phi)*sin(thet);
     rot(1,3,i)=sin(phi);
     rot(2,1,i)=-sin(thet);
     rot(2,2,i)=cos(thet);
     rot(2,3,i)=0;
     rot(3,1,i)=-sin(phi).*cos(thet);
     rot(3,2,i)=-sin(phi).*sin(thet);
     rot(3,3,i)=cos(phi);
     x=[ulx(j,i) vlx(j,i) wlx(j,i)];
     rotx=rot(1:3,1:3,i);
     xr=rotx*x';
     ur(j,i)=xr(1);
     vr(j,i)=xr(2);
     wr(j,i)=xr(3);
     spdh(j,i)=ur(j,i);  %magnitude of 1-hr mean streamwise wind vector
     dirx(j,i)=atan2(-vm,um)*180/pi;  %wind direction relative to boom point direction (about N)
     if (dirx(j,i)<0) 
       dirx(j,i)=dirx(j,i)+360;
     end
     dirx(j,i)=dirx(j,i)+or_mag+decl;
     dirx(j,i)=mod(dirx(j,i)+360,360);
   else
     spdh(j,i)=def_u;
     dirx(j,i)=def_dir;
   end; %end if
   end; %end for i  
end;%end for j

j_st=floor(jad(1));
j_en=j_st+1;

if prtit==1;
% prt_stats;
  prt_stats_redo3_10min;

end;
 
clear jaz Tvs uraw vraw wraw tsraw Rhv Tice Tan Tdan org ncr Rld Rlu Rsd Rsu Tcd Tcu Tdd Tdu Vbrn;
clear Rldx Rlux Tcdx Tddx Tcux Tdux vlx ulx tsx wlx;
clear sgu sgv sgw sgT sgTv sgRh sgTan sgTdan sgTsn sgTice sgvbrn; 
clear sgRlu sgRld sgRsu sgRsd sgTcu sgTcd sgTdu sgTdd;
clear Tv_cnt u_cnt v_cnt w_cnt Rh_cnt Tice_cnt Tsn_cnt Tan_cnt Tdan_cnt or_cnt ncr_cnt rld_cnt;
clear rlu_cnt rsd_cnt rsu_cnt tcd_cnt tcu_cnt tdd_cnt tdu_cnt vbr_cnt ts_cnt;
clear jz z spdh dirx Ticex Tsrx Th_tn Tdh_tn Th Qx Rtdx Rtux Rsdx Rsux orgx ncrx Rhvx esx ur;
clear Vbrx Nrx Qdnx Tbrx Rbux Tsbx Tsnwx;
clear jazd;
end;% for jdy

