% timatch_5
% VOCALS 2008 :: 2008-10-17 :: Simon de Szoeke
% Aggregates/averages data in 5-minute windows.
disp(['timatch_5: ', cruise, ' ', year]);

fclose('all');

clear ztsg zt zsnk zo y x w vv vrsp uu usr un5 un1 un ue5 ue1 ue ts_tsg ts_snk ts tn tm_s tm_ot tm_n tm_in tcol taub ta_im ta_etl ta t sogy sogx
clear sogm5 sog sody5 sody sodx5 sodx sm skip_it sig_sp shp_spd shp_hed sal_tsg s3 s1 rh_im rh_etl rf relsp reldir rdcon qs qcol qa_im qa_etl qa psp_im psp1 psp prt_it pir1 pir org odec5 ode npz
clear n minn m lrgsm lrgs5 lrgs lrgm5 lrgcm lrgc5 lrgc lrg lon lm lat k jl jk jj jg jf jdx jd jazz jazy5 jazy1 jazy jazx5 jazx1 jazx jaz5 jaz1 jaz jax jaq jadson jadscs jadmns jad5 
clear jad0 j imuy imux imum5 imu imdm5 imd ii i hsb hr hnet hm hlb g frac flist3 flist2 flist f e dt dir1 d cogm5 cog c by bx ay ax ans V U Sm Rns
clear Rnl Lm Dt

% inherit saveit from calling flux_eval program
% saveit=true;     %set true to print data to file
clearit=false;     %set true to clear variables at the end (clear same ones cleared above)

% jazx=zeros(6,1440);
% jaz=zeros(16,1440);
% jazy=zeros(12,1440);

jd=sprintf('%03i',ddd);
f=[way_proc_data_flux, 'proc_son' jd  '.txt'];
% flist=fopen(f,'r');
% jazx=fscanf(flist,'%g %g %g %g %g %g',[6,inf]);
jazx=load(f)';

e=[way_proc_data_flux, 'proc_scs' jd  '.txt'];
% flist2=fopen(e,'r');
% jaz=fscanf(flist2,'%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g',[16,inf]);
jaz=load(e)';

g=[way_proc_data_flux, 'proc_pc' jd  '.txt'];
% flist3=fopen(g,'r');
% jazy=fscanf(flist3,'%g %g %g %g %g %g %g %g %g %g %g %g',[12,inf]);
jazy=load(g)';

zsnk=ones(1,288);ztsg=ones(1,288);

[ax,bx]=size(jazx);			%jazx
[c,d]=size(jaz);				%jaz
[ay,by]=size(jazy);
tn=60*24;														%number of minutes in a day
n=1:tn;
jad0=str2double(jd)+(n-1)/24/60;
jad0(tn)=jad0(tn)+1e-5;		%	don't ask
jadson=jazx(1,1:bx);
jadscs=jaz(1,1:d);
jadmns=jazy(1,1:by);

for ii=2:length(jadscs);
   if jadscs(ii)<0,jadscs(ii)=jadscs(ii-1)+.000001;end;
end;

%%%%%%%%**********   sonic loop
clear jazx1;
jazxp=interp1(jadson,jazx',jad0);
jazx1=jazxp';

df=0; % skip the following if...
if df
    %fill in beginning empty spots
    jk=1;
    while jad0(jk)<jadson(1),
        jazx1(:,jk)=jazx(:,1);
        jk=jk+1;
    end;

    for i=jk:tn;
        t=find(jadson>jad0(i));
        if ~isempty(t),
            j=t(1);
        else
            j=i;
        end;
        dt=jadson(j)-jadson(j-1);
        Dt=jad0(i)-jadson(j-1);
        jazx1(:,i)=jazx(:,j-1)*(1-Dt/dt)+jazx(:,j)*Dt/dt;
        if j==bx;
            break;
        end;
    end;

    jk=tn;
    while jad0(jk)>=jadson(bx),
        jazx1(:,jk)=jazx(:,bx);
        jk=jk-1;
    end;
end; % if df, not done

%%%%%%%%*********  means loop
clear jazy1;
%iis=find(jazy(7,:)==99999);jazy(7,iip)=NaN;%change ta bad value

%fill in beginning empty spots 
jk=1;
while jad0(jk)<jadmns(1),
   jazy1(:,jk)=jazy(:,1);
   jk=jk+1;
end;

for i=jk:tn;
   t=find(jadmns>jad0(i));
   if ~isempty(t)
      j=t(1);
   else
      j=2;
   end;
   dt=jadmns(j)-jadmns(j-1);Dt=jad0(i)-jadmns(j-1);
   jazy1(:,i)=jazy(:,j-1)*(1-Dt/dt)+jazy(:,j)*Dt/dt;
   if j==by;break;end;
end;
jk=tn;
while jad0(jk)>=jadmns(by),
   jazy1(:,jk)=jazy(:,by);
   jk=jk-1;
end;
%jazy1=interp1(jadmns,jazy',jad0)';

% ship data loop

%first setup to interpolate fricking angles
cog=jaz(4,:);
sog=jaz(5,:);
lrg=jaz(11,:);
imu=jaz(12,:);
imd=jaz(13,:);

iic=sog<0 | sog>10;
cog(iic)=NaN; sog(iic)=NaN; lrg(iic)=NaN;
iid=imu<0 | imu>50;
imu(iid)=NaN; imd(iid)=NaN;

rdcon=pi/180;
jaz(4,:)=sog.*cos(cog*rdcon);			%sogx
jaz(5,:)=sog.*sin(cog*rdcon);			%sogy
jaz(12,:)=imu.*cos(imd*rdcon);		%imux
jaz(13,:)=imu.*sin(imd*rdcon);		%imuy

lrgc=cos(lrg*rdcon);
lrgs=sin(lrg*rdcon);

clear jaz1;
												%fill in beginning empty spots 
jk=1;
while jad0(jk)<jadscs(1),
   jaz1(:,jk)=jaz(:,1);
   lrgcm(jk)=lrgc(1);
   lrgsm(jk)=lrgs(1);
   jk=jk+1;
end;

for i=jk:tn;
   t=find(jadscs>jad0(i));
   if~isempty(t)
      j=t(1);
      dt=jadscs(j)-jadscs(j-1);Dt=jad0(i)-jadscs(j-1);
      jaz1(:,i)=jaz(:,j-1)*(1-Dt/dt)+jaz(:,j)*Dt/dt;
      lrgcm(i)=lrgc(j-1)*(1-Dt/dt)+lrgc(j)*Dt/dt;
      lrgsm(i)=lrgs(j-1)*(1-Dt/dt)+lrgs(j)*Dt/dt;
   else
      break;
   end;
end;

jk=tn;

while jk > 0 && jad0(jk)>=jadscs(d),
   jaz1(:,jk)=jaz(:,d);
   lrgcm(jk)=lrgc(d);
   lrgsm(jk)=lrgs(d);
   jk=jk-1;
end;

% do 5 minute averages
jax=1440;
jk=288;

ode=jaz1(15,:);%if no ODEC, use SOG or NaN

sodx=ode.*lrgcm;
sody=ode.*lrgsm;
uu=jazx1(2,:)/.95;				%.95 is empirical distortion factor
vv=jazx1(3,:)/1.15;				%1.15 is other empirical distortion factor

un=uu.*lrgcm+vv.*lrgsm;
ue=uu.*lrgsm-vv.*lrgcm;

clear jaz5 jazy5 jazx5 lrgc5 lrgs5 sodx5 sody5 un5 ue5 vrsp jad5;
iip=find(jazy(2,:)==-999);jazy(2,iip)=NaN;%change pir bad value

for j=1:jk
   jl=5*(j-1)+1;jj=min(jax,5*j);
   jazx5(:,j)=mean(jazx1(:,jl:jj),2)';
   jazy5(:,j)=mean(jazy1(:,jl:jj),2)';
   jaz5(:,j)=mean(jaz1(:,jl:jj),2)';
   lrgc5(j)=mean(lrgcm(jl:jj),2)';
   lrgs5(j)=mean(lrgsm(jl:jj),2)';
   sodx5(j)=mean(sodx(jl:jj),2)';
   sody5(j)=mean(sody(jl:jj),2)';
   un5(j)=mean(un(jl:jj),2)';
   ue5(j)=mean(ue(jl:jj),2)';
   vrsp(j)=var(jaz1(5,jl:jj))+var(jaz1(6,jl:jj));
   jad5(j)=jad0(jl);
end;
   
%*********************  then uncomponentize angles
U=jazx5(2,:);										
V=jazx5(3,:);										
W=jazx5(4,:);
reldir=mod(atan2(-V,U)/rdcon+180+360,360);		%signs convert jeff coord to conventional
relsp=sqrt(U.*U+V.*V);
sogx=jaz5(4,:);
sogy=jaz5(5,:);
imux=jaz5(12,:);%imet true wind
imuy=jaz5(13,:);%imet true wind

sogm5=sqrt(sogx.*sogx+sogy.*sogy);					%1-min sog
cogm5=mod(atan2(sogy,sogx)/rdcon+360,360);
imum5=sqrt(imux.*imux+imuy.*imuy);					%1-min imet speed, m/s
imdm5=mod(atan2(imuy,imux)/rdcon+360,360);
lrgm5=mod(atan2(lrgs5,lrgc5)/rdcon+360,360);			%heading system
odec5=sqrt(sodx5.*sodx5+sody5.*sody5);

%**********  do some true wind calcs
un1=un5+sogx;     
ue1=ue5+sogy;     
s1=sqrt(un1.*un1+ue1.*ue1);
s3=sqrt(un1.*un1+ue1.*ue1+W.*W);
s_shp=sqrt(imux.*imux+imuy.*imuy);

dir1=mod(atan2(-ue1,-un1)*180/pi+360,360);%- sign converts to to from
dir_shp=mod(atan2(imuy,imux)*180/pi+360,360);

% sea temps
ts_snk=jazy5(6,:); %correct tsnake
ts_tsg=jaz5(6,:);

sal_tsg=jaz5(7,:);
ts_tsg(ts_tsg<0)=NaN;
ts_snk(ts_snk<0)=NaN;
ts=ts_snk;							%water temperature if snake is in
% leave for auditability of snake:
%ts(isnan(ts_snk))=ts_tsg(isnan(ts_snk));  % use tsg if snake is invalid
zt=zsnk*.05;						%depth of water temp measurement
zt(isnan(ts_snk))=5.0;
qs=qsea(ts);

% air temps
ta_etl=jazy5(7,:)-0.05;%est correction for etl
ta_im=jaz5(9,:);
ta=ta_etl;
iita=find( ta_etl>50 | ta_etl<0);ta(iita)=NaN;ta_etl(iita)=NaN;

% pressure
pressm=jazy5(12,:);

% relative humidity
rh_etl=jazy5(8,:)+1;%est correct
rh_im=jaz5(8,:);
iih=find(rh_im<0);rh_im(iih)=NaN;
qa_etl=qair3([pressm' ta_etl' rh_etl'])';
qa_im=qair3([pressm' ta_im' rh_im'])';
qa=qa_etl;

% solar fluxes
psp1=jazy5(3,:);
psp_im=jaz5(10,:);
psp=psp1;
ii=find(psp<-50);psp(ii)=psp_im(ii);

% IR fluxes
pir=jazy5(2,:);
pir_im=jaz5(16,:);

org=jazy5(9,:);
lat=jaz5(2,:);
lon=jaz5(3,:);

shp_spd=odec5;
shp_hed=lrgm5;
sig_sp=sqrt(vrsp);

for k=1:jk
   %   u    0  ts  ta      qs     qa   rs      rl   rain     zi  p   zu    zt   zq
   x=[s1(k) 0 ts(k) ta(k) qs(k) qa(k) psp(k) pir(k) org(k) 600 p zus zts zqs lat(k) 1 0 1 1];
   y=cor30a(x);
   hsb(k)=y(1);
   hlb(k)=y(2);
   taub(k)=y(3);
   tcol(k)=y(11);
   qcol(k)=y(12);
   rf(k)=y(14);
   usr(k)=y(8);
   zo(k)=y(4);
end;

Rns=psp*.955;
Rnl=.97*(pir-5.67e-8*(ts_snk+273.15).^4);
hnet=Rns+Rnl-hsb-hlb-rf;
%iihm=find(~isnan(hnet));
hm=nanmean(hnet);sm=nanmean(hsb);lm=nanmean(hlb);Sm=nanmean(Rns);Lm=nanmean(Rnl);

if plotit
    figure;plot(jad5,hnet);xlabel(['Julian Day (', year, ')']);ylabel('Net Heat Flux (W/m^2)');title(['Heat Fluxes:  Net= ',num2str(hm),'; H_l_a_t= ',num2str(lm),'; H_s_e_n_s= ',num2str(sm),'; R_s_h_o_r_t= ',num2str(Sm),'; R_l_o_n_g= ',num2str(Lm)]);xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\Net_heat_flux_' num2str(ddd) '.jpg']); end
    figure;plot(jad5,s1,jad5,s_shp,jad5,relsp,'.');xlabel(['JD ', year]);title(['Wind Speed (m/s). ', cruise, year]);legend('PSD true', 'IMET true', 'PSD relative');xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\True_Wnd_speed_' num2str(ddd) '.jpg']); end
    figure;plot(jad5,dir1,jad5,dir_shp,jad5,reldir,'.');xlabel(['JD ', year]);title(['Wind dir (deg). ', cruise, year]);legend('PSD true', 'IMET true', 'PSD relative');xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\True_Wnd_dir_' num2str(ddd) '.jpg']); end
    figure;plot(jad5,hsb, jad5, hlb,jad5, Rnl,jad5,rf,'g.');xlabel(['JD ', year]);title('Heat Flux Components (W/m^2)');legend('H_s_e_n_s', 'H_l_a_t', 'Rnet_l_o_n_g', 'Rain','location','southwest','orientation','horizontal');xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\Heat_flux_components_' num2str(ddd) '.jpg']); end
    figure;plot(jad5,usr);xlabel(['JD ', year]);ylabel('U*(m/s)');title(['Friction Velocity COARE algorithm 3.0. ', cruise, year]);xlim([ddd ddd+1])
    if prtit; print('-djpeg90 ',[way_images_flux '\Friction_velocity_' num2str(ddd) '.jpg']); end
end

if saveit;
    np=jk;
    prt_GASEXIII_08_flux_5;
end;

if clearit
    clear ztsg zt zsnk zo y x w vv vrsp uu usr un5 un1 un ue5 ue1 ue ts_tsg ts_snk ts tn tm_s tm_ot tm_n tm_in tcol taub ta_im ta_etl ta t sogy sogx
    clear sogm5 sog sody5 sody sodx5 sodx sm skip_it sig_sp shp_spd shp_hed sal_tsg s3 s1 rh_im rh_etl rf relsp reldir rdcon qs qcol qa_im qa_etl qa psp_im psp1 psp prt_it pir1 pir org odec5 ode npz
    clear np n minn m lrgsm lrgs5 lrgs lrgm5 lrgcm lrgc5 lrgc lrg lon lm lat k jl jk jj jg jf jdx jd jazz jazy5 jazy1 jazy jazx5 jazx1 jazx jaz5 jaz1 jaz jax jaq jadson jadscs jadmns jad5 
    clear jad0 j imuy imux imum5 imu imdm5 imd ii i hsb hr hnet hm hlb g frac flist3 flist2 flist f e dt dir1 d cogm5 cog c by bx ay ax ans V U Sm Rns
    clear Rnl Lm Dt
end;