disp('fix_fastex')
%load C:\temp\finsnake.prn
%load c:\temp\ybulk.prn


%x=[4.5 0 29.2 27.6 24.8 18.1 203 413 0 600 1010 15 15 15];
n=3690;
y=zeros(n,21);
tims=finsnake(:,1);%time in big file
jday=tims;
u=finsnake(:,2);%relative wind speed
tss=finsnake(:,4);%sea surf temp
t=finsnake(:,5);%air temp
qss=finsnake(:,6);%sea surf spec hum
q=finsnake(:,7)*1.2;%air spec hum
rs=finsnake(:,8);%downward solar flux
rl=finsnake(:,9);%downward longwave flux
rain=finsnake(:,10);%rain rate (partial record)
zi=finsnake(:,11);%height of inversion (assumed)
p=finsnake(:,12);%atmos pressure
zu=finsnake(:,13);%height of anemmometer
zt=finsnake(:,14);%height of slow temp sensor
zq=finsnake(:,15);%height of mean humidity sensor
timb=ybulk(:,1);%time from small file
tsb=ybulk(:,4);%sea surf temp
tauc=-finsnake(:,44)*1.27;%covariance stress
tauc=sign(tauc).*min(abs(tauc),9.99);
hsc=finsnake(:,46)*1250;%covariance sensible heat
hlc=finsnake(:,47)*2500*1.27;%covariance latent heat (wq in g/kg m/s)
hlc=sign(hlc).*min(abs(hlc),999);
taui=finsnake(:,74).*finsnake(:,74)*1.27;
hsi=-finsnake(:,74).*finsnake(:,76)*1250;
hli=-finsnake(:,74).*finsnake(:,78)*1.27*2500*1000;
taui=sign(taui).*min(abs(taui),9.99);
hli=sign(hli).*min(abs(hli),999);
hsi=sign(hsi).*min(abs(hsi),999);
plat=finsnake(:,62);%rms motion corrections

shphead=mod(180/pi*finsnake(:,37)+360,360);%mean gyrocompass head (1 hr time res)
reldirx=mod(180/pi*finsnake(:,43)+360,360);%rel wind dir in file
uyr=u.*cos(reldirx*pi/180);%assumes reldirx computed atan2(x,y)
uxr=u.*sin(reldirx*pi/180)*1.05;%assumes 5% distortion on bow
reldir=mod(atan2(-uyr,-uxr)*180/pi+360,360);
unr=uxr.*cos(shphead*pi/180)+uyr.*sin(shphead*pi/180);%rotate to n
uer=uxr.*sin(shphead*pi/180)-uyr.*cos(shphead*pi/180);%rotate to e

un=finsnake(:,40);%n wind component, miscorrected for ship speed?
ue=finsnake(:,41);%e wind component, miscorrected for ship speed?

%correct cog and wind speed
coga=finsnake(:,39);%original cog in file
us=finsnake(:,38);%ship speed

%compute n/e components
use=us.*cos(coga);%reverse because coga calculated backwards in atan2
usn=us.*sin(coga);%reverse because coga calculated backwards in atan2
cog=mod(atan2(use,usn)*180/pi+360,360);%correct cog

%compute true wind speed/dir
utn=unr+usn;%n true wind, I think 
ute=uer+use;%e true wind
ut=sqrt(utn.^2+ute.^2);
wdir=mod(atan2(-ute,-utn)*180/pi+360,360);%true wind direction, meteo convention

%recompute flux time series, create edited dataset
j=1;
k=1;
kk=1;
for i=1:n
   if jday(i)<100
      jday(i)=jday(i)+365;  %blend jul day across year
   end;
   
   if tims(i)==timb(j)  %substitute ship ts for snake
      ts(i)=tsb(j);
      j=j+1;
      qs(i)=qsea(ts(i));
   else                  %keep snake
      ts(i)=tss(i);
      qs(i)=qss(i);
   end%end if tims
   if rl(i)>400;rl(i)=330;end;%despike rl data 
   x=[ut(i) 0 ts(i) t(i) qs(i) q(i) rs(i) rl(i) rain(i) zi(i) p(i) zu(i) zt(i) zq(i)];
   yy=coul_lkb(x);
   %     y=[hsb hlb tau zo zot zoq L usr tsr qsr dter dqer tkt RF hl_webb Cd Ch Ce Cdn_10 Chn_10 Cen_10];   
   %         1    2   3  4  5   6  7  8   9   10  11   12   13 14   15    16 17 18  19       20    21 
   
   y(i,:)=yy;
   hs(i)=yy(1);
   hl(i)=yy(2);
   tau(i)=yy(3);
   %set motion correction limits
   if ut(i)<10
      pltchk=1;
   else;
      pltchk=1+.5/10*(ut(i)-10);
   end;
   
   %build filtered data set
   if us(i)<2%shipspeed less than 2 m/s
      if reldir(i)<60 | reldir(i) >300 %rel wind direct +- 60 of bow
         if plat(i)<pltchk %ship motion corrections 
            if rain(i)<5   %rainrate 
               if abs(tauc(i))<9
               	if q(i)>1
            tim_f(k)=jday(i);
         	u_f(k)=ut(i);
         	tauc_f(k)=tauc(i);
         	taub_f(k)=tau(i);
         	hsc_f(k)=hsc(i);
         	hsb_f(k)=hs(i);
         	hlc_f(k)=hlc(i);
            hlb_f(k)=hl(i);
            taui_f(k)=taui(i);
            hsi_f(k)=hsi(i);
            hli_f(k)=hli(i);
            sigq_f(k)=finsnake(i,55);
            sigop_f(k)=finsnake(i,60);
         	mnop_f(k)=finsnake(i,61);
         	plat_f(k)=finsnake(i,62);
            if sigop_f(k)<25 & mnop_f(k)>2500
               tim_ff(kk)=tims(i);
               hlc_ff(kk)=hlc(i);
               hlb_ff(kk)=hl(i);
               hli_ff(kk)=hli(i);
               u_ff(kk)=ut(i);
               qs_ff(kk)=qs(i);
               q_ff(kk)=q(i);
               ts_ff(kk)=ts(i);
               t_ff(kk)=t(i);
               sigq_ff(kk)=finsnake(i,55);
               kk=kk+1;
            end;
            k=k+1;
            end;%end q
            end;%end tauc
            end;%end rain
      	end;   
      end;%end if reldir
   end;%end if us
   
end;%end for i
m=k-1;
figure;plot(taub_f(1:m),taub_f(1:m),taub_f(1:m),tauc_f(1:m),'.');axis([-.25 1.5 -.25 1.5]);xlabel('Bulk Stress (N/m^2)');ylabel('Cov Stess (N/m^2)')
figure;plot(tim_f(1:m),taub_f(1:m),tim_f(1:m),tauc_f(1:m),'.');axis([357 396 -.2 1.4]);xlabel('Julian Day (1996)');ylabel('Stess (N/m^2)')

   