disp('metz_red3');

load('c:\data\cwf\Dropbox\DOC\CLIMATE\oceobs\GSOP2012\s1to10_best_met_1hr.mat');%read file with hour average data; set your local directory 
qbo=qair_p3([atmp bpr hrh/100]);
%          x=load('c:\data\cwf\data\neaqs\flux\neaqs04flux_hr_3b.txt');%read file with hour average data; set your local directory 
%       x=load('c:\data\cwf\dropbox\data\stratus10\stratus2010_flux_10_cf2.txt');%read file with hour average data; set your local directory 
   %x=load('c:\data\cwf\data\texaqs_06\flux\flux_30hf_Texaqs2006a.txt');
          yx=load('c:\data\cwf\dropbox\data\VOCALS_2008\STRATUS_2001_2010_climatology_hr_4.txt');%read file with hour average data; set your local directory 
jdday=yx(:,1)+datenum(yx(:,55),1,1)-datenum(1999,12,31)-1;
x=interp1(jdday,yx,yday);
jdy=x(:,1);%julian day at beginning of time average
ushp=x(:,2);%doppler log, SCS (m/s)
U=x(:,3);%true wind,PSD sonic (m/s)
dir=x(:,4);%true wind direction, PSD sonic (deg)
urel=x(:,5);%relative wind speed, PSD (m/s)
reldir=x(:,6);%relative wind direction (from),clockwise rel ship's bow, PSD sonic (deg)
head=x(:,7);%ship heading, deg clockwise rel north, SCS laser ring gyro

tsnk=x(:,8);%sea snake temperature, PSD, 0.05 m depth (C)
ta=x(:,9);%air temperature, PSD (C)
qse=x(:,10);%sea surface specific humidity, from snake (g/kg)
qa=x(:,11);%air specific humidity, PSD (g/kg)
hsc=x(:,12);%sensible heat flux, covariance, PSD sonic anemometer(W/m^2)
hsib=x(:,13);%sensible heat flux, ID, PSD sonic anemometer(W/m^2)
hsb=x(:,14);%bulk sensible heat flux, (W/m^2)
hlc=x(:,15)*1.04;%latent heat flux, covariance, (W/m^2)  4% for sensor separation
hlib=x(:,16);%latent heat flux, ID, (W/m^2)
hlb=x(:,17);%bulk latent heat flux, W/m^2 (includes Webb et al. correction)
taucx=x(:,18);%covariance streamwise stress, PSD sonic anemometer (N/m^2)
taucy=x(:,19);%covariance cross-stream stress, PSD sonic anemometer (N/m^2)
tauib=x(:,20);%ID streamwise stress, PSD sonic anemometer (N/m^2)
taub=x(:,21);%bulk wind stress along mean wind, (N/m^2)
rs=x(:,22);%downward solar flux, PSD units (W/m^2)
rl=x(:,23);%downward IR flux, PSD units (W/m^2)
org=x(:,24);%rainrate, PSD STI optical rain gauge, uncorrected (mm/hr)
J=x(:,25);%ship plume contamination index
tiltx=x(:,26);%flow tilt at PSD sonic anemometer, earth frame
Jm=x(:,27);%ship maneuver index
ct=x(:,28);%ct^2 (K^2/m^.667)
cq=x(:,29);%cq^2 ((g/kg)^2/m^.667)
cu=x(:,30);%cu^2 ((m/s)^2/m^.667)
cw=x(:,31);%cw^2 ((m/s)^2/m^.667)
hrain=x(:,32);%rain heat flux,Gosnell et al 1995, JGR, 18437-18442, (W/m^2)
hlwebb=x(:,33);%correction to measured latent heat flux, Webb et al.

lats=x(:,34);%latitude, deg (SCS pcode)
lons=x(:,35);%longitude, deg (SCS pcode)
zu_etl=x(:,36);%height of mean wind sensor, 17.7 m
zt_etl=x(:,37);%height of mean air temperature sensor, 15.5 m
zq_etl=x(:,38);%height of mean air humidity sensor, 15.5 m
%***** ships imet and scs data
sog=x(:,39);%speed over ground, SCS gps, (m/s)
U_scs=x(:,40); %true wind speed, imet propvane anemometer (m/s)
dir_scs=x(:,41);%true wind direction (from),clockwise rel north, imet,(deg)
cog=x(:,42);%%course over ground, SCS gps, (m/s)
tsg=x(:,43);%tsg water temperature, 5 m depth, (C)
ta_im=x(:,44);%imet air temperature (C)
qs_tsg=x(:,45);%imet bulk water specific humidity (g/kg)
qa_im=x(:,46);%imet air specific humidity, (g/kg)
rs_im=x(:,47);%imet solar flux, (W/m^2)
rl_im=x(:,48);%imet IR flux (W/m^2) – not connected for neaqs
wco2_lic=x(:,49);%LICOR CO2 flux, (micatm m/s)
q_lic=x(:,50);%Specific humidity from LICOR (g/kg)
sgq_lic=x(:,51);%Standard deviation of specific humidity from LICOR (g/kg)
co2_lic=x(:,52); %CO2 concentration from Licor (umol/mol)
sgc_lic=x(:,53);%Standard deviation of specific humidity from LICOR (g/kg)
year=x(:,55);
press=x(:,54); %Atmospheric pressure (mb)
%qse=qsee([tsnk press]);
xco2=co2_lic;
%  ii=find(lat>100);lat(ii)=lat(ii-1);
% ii=find(lon>-70);lon(ii)=lon(ii+1);
 
       [n m]=size(x);
%    yr=ones(n,1)*2010;
    
    
   
   
   q_oph=NaN;
   ophm=NaN;
   sigoph=NaN;
   hlib_lic=hlib;
   
   qp=qa;
   tp=ta;
   aoz=1e3;%waterside reaction rate, s^-1
    H=3;%henry constant = 1/alpha
    alph=1/H;%dimensionless solubility
    scw=500;%water side schmidt #
   clear y yy yy5 qo;
tdk=273.16;
grav=grv(lats);%9.82;
Rgas=287.1;
cpa=1004.67;
rhoair=1010*100./(Rgas*(ta+tdk).*(1+0.61*qa/1000));
visa=1.326e-5*(1+6.542e-3*ta+8.301e-6*ta.*ta-4.84e-9*ta.*ta.*ta);

%%%%%Recompute Fluxes
rh=relhum5([ta qa press]);
P=press;
zref_u=4.;zref_t=2.75;zref_q=2.75;zi=600;
A=[];A=coare35vn_cool_refa(U,zu_etl,ta,zt_etl,rh,zq_etl,P,tsnk,rs,rl,lats,zi,org,NaN,NaN,zref_u,zref_t,zref_q);
%A=[usr tau hsb hlb hbb hsbb hlwebb tsr qsr zou zot zoq Cd Ch Ce  L zet ug dter dqer tkt Urf Trf Qrf RHrf UrfN Rnl Le rhoa UN U10 U10N Cdn_10 Chn_10 Cen_10 RF Qs Evap];
%   1   2   3   4   5   6    7      8   9    10  11  12 13 14 15  16 17 18  19  20  21  22  23   24   25 26  27  28  29  30    31     32     33   34 35  36
usb=A(:,1);
taub=A(:,2);
hsb=A(:,3);
hlb=A(:,4);
usr=A(:,1);
tsr=A(:,8);
qsr=A(:,9);
zou=A(:,10);
zot=A(:,11);
zoq=A(:,12);
L=A(:,16);
Cd=A(:,13);
Ce=A(:,14);
ug=A(:,18);
Gs=sqrt(U.^2+ug.^2)./U;
U10n=usr./.4./sqrt(Gs).*log(10./zou);
DT10n=-tsr./.4.*log(10./zot);
Dq10n=-qsr./.4.*log(10./zoq);
uref=A(:,22);
Tref=A(:,23);
qref=A(:,24);



% jj=find(~isnan(U));    yy=NaN*ones(length(yday),22);
%     for i=1:length(jj);  
%    xx=[U(jj(i)) 0 tsnk(jj(i)) ta(jj(i)) qse(jj(i)) qa(jj(i)) rs(jj(i)) rl(jj(i)) org(jj(i)) 600  1010 zu_etl(jj(i)) zt_etl(jj(i)) zq_etl(jj(i)) lats(jj(i)) 1 0 1 1] ;		%set data for basic flux alogithm
%    y=cor30a(xx);
%    yy(jj(i),:)=y;
%    %xx5=[u us ts t qs q Rs Rl rain zi  P zu zt zq lat(ibg) jcool] ;
%    xx5=[U(i) 0 tsnk(i) ta(i) qse(i) qa(i) rs(i) rl(i) org(i) 600  1010 zu_etl(i) zt_etl(i) zq_etl(i) lats(i) 1] ;
%    %y=cor25(xx5);
%    %yy5(i,:)=y;
%     usr(jj(i))=y(8);hsb(jj(i))=y(1);hlb(jj(i))=y(2);
%     usrc(jj(i))=sqrt(abs((taucx(jj(i))+tauib(jj(i)))/2)./rhoair(jj(i)));
% 
%     if (reldir(jj(i))<100 | reldir(jj(i))>260) & ~isnan(usrc(jj(i))) & (abs(taucx(jj(i))-tauib(jj(i)))<taub(jj(i)))
%         usx=usrc(jj(i));
%     else
%         usx=usr(jj(i));
%     end;
%     usx=usr(jj(i));
% %     xxx=[U(i) tsnk(i) ta(i) qa(i) rl(i) 600 1010 zu_etl(i) usx hsb(i) hlb(i)  aoz alph scw];
%     yo=cor30_ks_oz_a(xxx);
%     qo(i,:)=yo;

% end;

% vt03=100*qo(:,5);
% ii=find(vt03>1e1);vt03(ii)=NaN;
% zou=yy(:,4);
% zoq=yy(:,6);
% cdn10=yy(:,19);
% cen10=yy(:,21);
% nnx=length(zou);%nnx=3508
% hlb(1:nnx,1)=yy(1:nnx,2);%account for old algorithm
% taub(1:nnx,1)=yy(1:nnx,3);%account for old algorithm
% hsb(1:nnx,1)=yy(1:nnx,1);%account for old algorithm
dtery=A(:,17);
%dqery=yy(:,12)*1e3;   

Lon=lons;Lat=lats;
Ts=tsnk;   
zts=.05*ones(length(jdy),1);
Hsb=hsb;Hlb=hlb;Taub=taub;
warn_calc2;
   %*******************  set constants  ****************

iv=3;
ii=find(isfinite(tsnk));Le=(2.501-.00237*median(tsnk(ii)))*1e6;
rlclr=(.52+.13/60*abs(lat)+(.082-.03/60*abs(lats)).*sqrt(qa)).*(5.67e-8*(ta+tdk).^4); 
stdt=10;endt=35;
figure;;
epic01_Rsclear_5;%epic99_da_Rsclear_hr;%get clear sky solar flux data
rsclr=Rscl;
rlnet=0.97*(rl-5.67e-8*(tsnk+tdk).^4);
rlnetc=0.97*(rlclr-5.67e-8*(tsnk+tdk).^4);
dtery=dt(:,21);
dqery=dt(:,22)*1e3;   

cu=.00087;bu=.00018;cx=.00085;bx=.00015;
dts=tsnk-ta-dtery;
dqs=qse-qa-dqery;
dtvs=dts+0.61e-3.*(ta+tdk).*dqs;



wtb=hsb./rhoair./cpa;
wqb=hlb./rhoair./Le*1e3;
wtvb=wtb+.61e-3.*wqb.*(ta+tdk);
wtb_son=wtb+.51e-3.*wqb.*(ta+tdk);
np=length(jdy);
for i=1:np
  	if wtvb(i)>0   
     	wstr(i)=  ( grav(i)*wtvb(i)*600./(ta(i)+tdk))^.333;
  	else
  		wstr(i)=.2/1.2;   
  	end;
end;
    
wgust=1.2*wstr';
gf=sqrt(1+wgust.^2./U.^2);
usb=sqrt(taub./rhoair.*gf);%
tsb=-wtb./usb;
tsbs=-wtb_son./usb;
qsb=-wqb./usb;
   
wbar=1.61*hlb/Le./(1+1.61*qa/1e3)./rhoair+hsb./rhoair/cpa./(ta+tdk);%formulation in hlb already includes webb
hl_webb=rhoair.*wbar.*qa*Le/1000;
ii=find(isnan(hlwebb));
hlwebb(ii)=hl_webb(ii);

dter=dt(:,21);
dqer=dt(:,22)*1e3;   
zetx=-.4*grav.*zu_etl.*wtvb./usb.^3./(ta+tdk);
Lb=zetx./zu_etl;
zetx=max(zetx,-50);
zetx=min(zetx,50);
u10n=U-usb/.4.*(log(zu_etl/10)-psiu_26(zetx))./gf;
q10n=qa-qsb/.4.*(log(zq_etl/10)-psit_26(zetx.*zq_etl./zu_etl));
dq10n=qse-q10n-dqer;
th10n=ta-tsb/.4.*(log(zq_etl/10)-psit_26(zetx.*zq_etl./zu_etl))+.0098*zq_etl;
dth10n=tsnk-th10n-dter;

ii=find(dq10n<-5);dq10n(ii)=qse(ii)-qp(ii);
 
   %**********************
ii=find(co2_lic<300);hlc(ii)=NaN;%check for dirty licor
hll=hlc;
ii=find(isnan(hlc) & isfinite(hlib));
hll(ii)=hlib(ii);
ii=find(isfinite(hlc) & isfinite(hlib));
hll(ii)=(hlib(ii)+hlc(ii))*.5;
ii=find(abs((hlc-hlib)./(hlib+hlc))>.4 | abs(abs(qsb).*phi_tt_nq(zetx)-sgq_lic)>.1);
hll(ii)=NaN;
yhsb=hsb;
yhsc=hsc;   
yhlc=hlc;
yhlb=hlb;
v1=u10n;
  %ii=find(abs((U-U_scs)./U)>.5);v1(ii)=NaN;
hl_webar=hlwebb';
yqs=qse;
yq=qa;
jdollars='stratus10';
   
   
% metz_med_tau_plt2;%_30;
% hlcn=hll*0;%no null
% metz_med_hl_plt7_str7;
% 
% figure;h=loglog(taubmm,abs(taucxmm),'o',taubmm,abs(tauibmm),'d',[1e-3 1],[1e-3 1]);axis([1e-3 1 1e-3 1]);set(gca,'fontsize',13);set(h,'LineWidth',2);xlabel('Bulk Stress (Nm^{-2})');ylabel('Turbulent Stress (Nm^{-2})');%title('covariance - cirlce; ID - diamond');
% 
% ii=find(isfinite(hll) & (reldir>270 | reldir<90) & isfinite(dq10n) & org<.3);
% figure;h=plot(u10n(ii).*gf(ii).*dq10n(ii)*1.13e-3*Le/1e3.*1.18,(hll(ii)+hlwebb(ii)),'.',[-0 300],[-0 300]);
% set(h,'LineWidth',2);axis([-0 300 -0 300]);set(gca,'fontsize',13);xlabel('C_{e10} S_{10n} \Deltaq_{10n} \rho_{a} Le  (Wm^{-2})');ylabel('Hl (Wm^{-2})');title([ 'stratus10']);
% 
%  

