disp('metz_red3');

%   x=load('c:\data\cwf\data\fluxdatas\metz_hr_3.txt');%read file with hour average data; set your local directory 
%          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\VOCALS_2008\STRATUS_2001_2010_climatology_hr_2.prn');%read file with hour average data; set your local directory 
   %x=load('c:\data\cwf\data\texaqs_06\flux\flux_30hf_Texaqs2006a.txt');
jdyx=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)
%qse=qsee([tsnk press]);
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)
x(:,15)=x(:,15)*1.00;%acount for licor-sonic separation
hlc=x(:,15);%latent heat flux, covariance, (W/m^2)
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.
lat=x(:,34);%latitude, deg (SCS pcode)
lon=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)
press=x(:,54); %Atmospheric pressure (mb)
sgC_lic=x(:,53);%Standard deviation of CO2 concentration from LICOR (microatm)
yearx=x(:,55);
% ii=find(press==1000);press(ii)=1015;
% ij=find(yearx==2006);tsnk(ij)=tsnk(ij)-0.6;
% qse=qsee([tsnk press]);%
jdy=jdyx+datenum(yearx,12,31)-datenum(2001,12,31);

%        ys=load('c:\data\cwf\data\gasex3\GasExIII 08flux_10_lis.txt');%read file with hour average data; set your local directory 
%        yn =load('c:\data\cwf\data\gasex3\GasExIII 08flux_10_lin.txt');%read file with hour average data; set your local directory 
%        y1x=load('c:\data\cwf\data\gasex3\GasExIII 08flux_10_li1.txt');%read file with hour average data; set your local directory 
%        y2x =load('c:\data\cwf\data\gasex3\GasExIII 08flux_10_li2.txt');%read file with hour average data; set your local directory 
%        y3x=load('c:\data\cwf\data\gasex3\GasExIII 08flux_10_li3.txt');%read file with hour average data; set your local directory 
%   


 

  wco2_lic_1=x(:,49);
  sgq_lic_1=x(:,51);
  hlc1=x(:,15);%latent heat flux, covariance, (W/m^2)
  hlib1=x(:,16);%latent heat flux, ID, (W/m^2)
  xco2_lic=x(:,52); %CO2 concentration from sample Licor C9 (umol/mol) 
  xco2a=xco2_lic;
  xco2=xco2a-10;
  %cq=y1(:,29);%cq^2 ((g/kg)^2/m^.667)

 % eqtemp=x(:,56); 
%   eqpress=x(:,57); 
%   xco2=x(:,58);  %aqueous Co2 (umol/mol)?  
%   xco2a=x(:,59); %air Co2 (umol/mol) 
%   fco2=x(:,60); %aqueous Co2 (umol/mol)? 
%   fco2a=x(:,61); %aqueous Co2 (umol/mol)? 
%   sss=x(:,62); %salinity 
%    wco2_lic_C6=x(:,63);

       [n m]=size(x);
    yr=ones(n,1)*2008;
    
    
   
   
   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(lat);%9.82;
Rgas=287.1;
cpa=1004.67;
rhoair=press.*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);
yyc=solco2(tsnk);solc=yyc(:,1);alc=yyc(:,2);schm=yyc(:,3);
    
    for i=1:length(x);  
   xx=[U(i) 0 tsnk(i) ta(i) qse(i) qa(i) rs(i) rl(i) org(i) 600  press(i) zu_etl(i) zt_etl(i) zq_etl(i) lat(i) 1 0 1 1] ;		%set data for basic flux alogithm
   y=cor30a(xx);
   yy(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  press(i) zu_etl(i) zt_etl(i) zq_etl(i) lat(i) 1] ;
       usr(i)=y(8);hsb(i)=y(1);hlb(i)=y(2);cd10n(i)=y(19);wbaral(i)=y(15);
       usx=usr(i);
     xxx=[xx(1) xx(3:4) xx(6) xx(8) xx(10:12) usx hsb(i) hlb(i) cd10n(i)];%form input for ozdep model
    %xxx=xxp(i,:);
    yo=cor30_ks_co2_ke(xxx);
    qoc(i,:)=yo;
    yo=cor30_ks_dms_ke(xxx);
    qoms(i,:)=yo;
    
    
end;


vtco2=qoc(:,5);%coare xfer velocity in m/s for CO2
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=yy(:,11);
dqery=yy(:,12)*1e3;   

Lon=lon;Lat=lat;
Ts=tsnk;   
zts=.05*ones(length(jdy),1);
Hsb=hsb;Hlb=hlb;Taub=taub;
warn_calc;
   %*******************  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(lat)).*sqrt(qa)).*(5.67e-8*(ta+tdk).^4); 
stdt=292;endt=339;
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=yy(:,11);
dqery=yy(:,12)*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;
   
wbarQ=1.61*hlb/Le./(1+1.61*qa/1e3)./rhoair;
wbarT=hsb./rhoair/cpa./(ta+tdk);%formulation in hlb already includes webb
wbar=wbarQ+wbarT;
     %wbar=1.61*hlb/Le/rhoa+(1+1.61*Q)*hsb/rhoa/cpa/ta;
hl_webb=rhoair.*wbar.*qa*Le/1000;
ii=find(isnan(hlwebb));
hlwebb(ii)=hl_webb(ii);

dter=yy(:,11);
dqer=yy(:,12)*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);

ct_dim=ct.*zu_etl.^.667./tsbs.^2;
figure; semilogy(zetx,ct_dim./psi_ft(zetx),'.',[-10 .5],[1 1]);axis([-2 .5 .1 100]);xlabel('z/L');ylabel('Dimensionless C_t^2/psi_ft');

cq_dim=cq.*zu_etl.^.667./qsb.^2;
figure; semilogy(zetx,cq_dim./psi_ft(zetx),'.',[-10 .5],[1 1]);axis([-2 .5 .1 100]);xlabel('z/L');ylabel('Dimensionless C_q^2/psi_ft');


%hlcs=hlc; hlibs=hlib; wco2_lic_s=wco2_lic;
 ij=find(abs(abs(qsb).*phi_tt_nq(zetx)-sgq_lic_1)>.15 | abs(hlc1-hlib1)>30);hlc1(ij)=NaN;hlib1(ij)=NaN;wco2_lic_1(ij)=NaN;
 
%ij=find(sgq_lic_s>25*abs(qsb).*phi_tt_nq(zetx));hlcs(ij)=NaN;hlibs(ij)=NaN;wco2_lic_s(ij)=NaN;
%ij=find(sgq_lic>25*abs(qsb).*phi_tt_nq(zetx));hlc1(ij)=NaN;hlib1(ij)=NaN;wco2_lic_1(ij)=NaN;


%kk=find(wco2_lic_1<-1.5 | sgC_lic_1>.3);wco2_lic_1(kk)=NaN;
wco2_lic=wco2_lic_1;
wco2_lic_n=zeros(size(wco2_lic_1));
   %**********************8
hll=hlc1;hlc=hlc1;hlib=hlib1;
ii=find(isnan(hlc1) & isfinite(hlib));
hll(ii)=hlib(ii);
ii=find(isfinite(hlc1) & isfinite(hlib));
hll(ii)=(hlib(ii)+hlc1(ii))*.5;
ii=find(abs(hlc1-hlib)>.3*abs(hlib+hlc1)+10);
hll(ii)=NaN;
hlcn=0.1*ones(size(hlc));%null licor latent heat
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='stratusClim';

qq=sgq_lic.^2;
jj=find(org<1 & reldir>270 | reldir<90);
figure;semilogy(zetx(jj),sqrt(qq(jj))./abs(qsb(jj)),'.',zetx(jj),phi_tt_nq(zetx(jj)),'.');axis([-10 2 0.1 5e1])
xlabel('z/L');ylabel('Dimensionless \sigma_q');

   %ophclr=ones(np,1);ophclr=ophclr*2600;
metz_med_tau_plt2;%_30;
metz_med_hl_plt7_str7;
sigcd=sigu.*yyu./taubmm;
sigce=sig./dqm/Le./dum*1e3./gfm;
cen10mt=yyz.*hlgm;
%save c:\data\cwf\doc\climate\wcrp\highlatitude\section\sogsex.txt dum cdn10mt nm sigcd cen10mt nn sigce  -tabs -ascii;


%nar_daily
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');
  
   
 %webster;  
   
u10o=U.*log(10/1e-4)./log(zu_etl/1e-4);
u10=max(2.5,u10o);
uinv=1./u10;
cdo=1e-3*(.862+.088*u10-.00089*u10.^2);
cd1=1e-3*(.1034-.00678*u10+.0001147*u10.^2);
cd=cdo+(tsnk-ta).*cd1;
   
tau_m=rhoair.*cd.*u10.^2;

clo=1e-3*(.994+.061*u10-.001*u10.^2);
cl1=1e-3*(-.02+.69*uinv-.817*uinv.^2);
cl=clo+(tsnk-ta).*cl1;
hl_m=rhoair*Le.*cl.*u10.*(qse-qa)/1000;
   
cs=0.96*cl;
hs_m=rhoair*cpa.*cs.*u10.*(tsnk-ta);
 

 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)-hlcn(ii)),'.',[-50 300],[-50 300]);
set(h,'LineWidth',2);axis([-50 250 -50 250]);set(gca,'fontsize',13);xlabel('C_{e10} S_{10n} \Deltaq_{10n} \rho_{a} Le  (Wm^{-2})');ylabel('Hl (Wm^{-2})');title([ '']);

%  ii=find(isfinite(hll) & (reldir>270 | reldir<90) & sog<2 & org<.3 );
% figure;h=plot(jdy(ii),wco2_lic(ii)+wbarT(ii).*xco2a(ii),'o',jdy,vtco2.*alc.*(xco2-xco2a),jdy(ii),wco2_lic_n(ii),'x');axis([60 100 -.6 .2]);
% set(h,'LineWidth',2);set(gca,'fontsize',13);xlabel('Julian Day (2008)');ylabel('<w*CO2> (micmole/mole*m/s)');title([ 'GasExIII']);
% 
%  jj=find(isfinite(hll) & (reldir>270 | reldir<90) & sog>2 & org<.3);
% figure;h=plot(U(ii),42*1e-6*3600*24*365*wco2_lic_n(ii),'o',U(jj),42*1e-6*3600*24*365*wco2_lic_n(jj),'x',[0 10],[0 0],':k');
% %figure;h=plot(U(ii),42*1e-6*3600*24*365*wco2_lic_C6(ii),'o',U(jj),42*1e-6*3600*24*365*wco2_lic_C6(jj),'x',U,2.03-.007*U,'--',[0 10],[0 0],':k');
% set(h,'LineWidth',2);set(gca,'fontsize',13);xlabel('Wind Speed (m/s)');ylabel('NULL <w*CO2> (mole/m^2/yr)');title([ 'GasExIII BIAS: o, sog<2; x, sog>2']);

ub=2:2:22;
for i=1:length(ub)-1; 
ii=find(isfinite(hll) & (reldir>300 | reldir<60) & sog<5 & org<.3 & u10n>ub(i) & u10n<ub(i)+2 & abs(dq10n)>1 & xco2>340 & abs(q_lic-qa)<5);
 xyl = prctile(hll(ii)+hlwebb(ii),[2.5 15 50 85 97.5]);xyb = prctile(hlb(ii),[2.5 15 50 85 97.5]);xyu = prctile(u10n(ii),[2.5 15 50 85 97.5]);
 xyn = prctile(hlcn(ii),[2.5 15 50 85 97.5]);
  xyc = prctile(wco2_lic(ii),[2.5 15 50 85 97.5]);
%  xyc1 = prctile(wco2_lic_1(ii),[2.5 15 50 85 97.5]);
%  xycn = prctile(wco2_lic_n(ii),[2.5 15 50 85 97.5]);
 xyco=prctile(vtco2(ii).*alc(ii).*(xco2(ii)-xco2a(ii)),[2.5 15 50 85 97.5]);
 ubm(i)=xyu(3);
 ylx(1:3,i)=xyl(2:4);
 ybx(1:3,i)=xyb(2:4);
 ybn(1:3,i)=xyn(2:4);
 ybc(1:3,i)=xyc(2:4);
%  ybc1(1:3,i)=xyc1(2:4);
%  ybcn(1:3,i)=xycn(2:4);
 ybco(1:3,i)=xyco(2:4);
end;

figure;plot(ubm,ylx(2,:)-ybn(2,:),'-x',ubm,ybx(2,:),'o');

ii=find(isfinite(hll) & (reldir>270 | reldir<90) & isfinite(dq10n) & org<.3  & sog<5 );
 [ss,pp]=polyfit(hlb(ii),(hll(ii)+hlwebb(ii)),1)
  mk=1:length(x);
 ii=find( (reldir>270 | reldir<90) & xco2a>340 & org<.3 & sog<2 );figure;plot(mk(ii),wco2_lic(ii)-wco2_lic_n(ii)+xco2a(ii).*wbarT(ii),'o',mk,vtco2.*alc.*(xco2-xco2a));
 
 Csb=-vtco2.*alc./usb.*(xco2-xco2a);
 %figure;semilogy(mk,x(:,51),'o',mk,abs(qsb).*phi_tt_nq(zetx),'k');xlabel('JDY');ylabel('sig_q');ylim([ 1e-2 1]);
 %figure;semilogy(jdy,sgC_lic,'o',jdy,abs(Csb).*phi_tt_nq(zetx),'k');xlabel('JDY');ylabel('sig_{CO2}');ylim([1e-5 1e1]);

 figure;semilogy(mk,x(:,51),'o',mk, abs(qsb).*phi_tt_nq(zetx),'k');xlabel('10-min units');ylabel('sig_q');ylim([ 1e-2 2]);
 
 Xm=[ 0  10  20 30 40 50 60 70 75 80 85 90 95 100 105 110 120 130 140 150 160 180 200 210 220];ii=find( (reldir>270 | reldir<90) & xco2a>340 & org<.3 & sog<5 );[N1,X1]=hist(hlb(ii)',Xm);[N2,X2]=hist(hll(ii)+hl_webar(ii)',Xm);figure;h=plot(X1(1:length(X1)-1),N1(1:length(N1)-1)./diff(X1),'-x',X2(1:length(X2)-1),N2(1:length(N2)-1)./diff(X2),'-or')
 set(h,'LineWidth',2);xlabel('Latent Heat Flux (W/m^2)');ylabel('Number density');legend('Bulk Flux','Turbulent Flux');
 