disp('epicz_red');

   x=load('d:\data\epic\flux30_sum_nohd.txt');%read file with hour average data; set your local directory 
   jdy=x(:,1);%julian day at beginning of time average
   U=x(:,2);%true wind,ETL sonic (m/s)
   dir=x(:,3);%true wind direction, ETL sonic (deg)
   tsnk=x(:,4);%sea snake temperature,ETL, 0.05 m depth (C)
   tsg=x(:,5);%tsg water temperature, 5 m depth, (C)
   sal=x(:,6);%tsg salinity (psu)
   ta=x(:,7);%air temperature, ETL (C)   
   qse=x(:,8);%sea surface specific humidity,from snake (g/kg)
   qa=x(:,9);%air specific humidity, ETL (g/kg)nic anemometer(W/m^2)
   rs=x(:,10);%downward solar flux, ETL units (W/m^2)
   rl=x(:,11);%downward IR flux,  ETL units (W/m^2)
   org1=x(:,12);%rainrate, ETL STI optical rain gauge, uncorrected (mm/hr)
   ushp=x(:,13);%doppler log, SCS (m/s)
   head=x(:,14);%ship heading, deg clockwise rel north, SCS laser ring gyro (deg)
   urel=x(:,15);%relative wind speed, ETL (m/s)
   reldir=x(:,16);%relative wind direction (from),clockwise rel ship's bow, ETL sonic (deg)
   latx=x(:,17);%latitude, deg  (SCS pcode)
   lonx=x(:,18);%longitude, deg (SCS pcode)
   zts=x(:,19);%depth for bulk flux Ts reference, =0.05 when snake is used
   sig_u=x(:,20);%std dev of ship speed, m/s (>.2 indicates maneuver during half hour average)
   taub=x(:,21);%bulk wind stress along mean wind, (N/m^2)
   hsb=x(:,22);%bulk sensible heat flux, (W/m^2)
   hlb=x(:,23);%bulk latent heat flux, W/m^2 (includes Webb et al. correction)
   hrain=x(:,24);%rain heat flux, as per Gosnell et al 1995, JGR, 18437-18442, (W/m^2)
   org2=x(:,25);
   org3=x(:,26);
   org4=x(:,27);
   org1_c=x(:,28);
   org2_c=x(:,29);
   org=org3;
   
   zu_etl=18.5;%height of mean wind sensor, 17.7 m
   zt_etl=15.5;%height of mean air temperature sensor, 15.5 m
   zq_etl=15.5;%height of mean air humidity sensor, 15.5 m
   zq_ass=6.7;
   %*****   ships imet and scs data
   sog=x(:,13);%speed over ground, SCS gps, (m/s)
   cog=x(:,14);%%course over ground, SCS gps, (m/s)
   
   lat=latdecode(latx); 
   lon=-latdecode(lonx);%negative for W long
   clear yy;
for i=1:length(x);  
   xx=[U(i) 0 tsnk(i) ta(i) qse(i) qa(i) rs(i) rl(i) org(i) 600  1010 zu_etl zt_etl zq_etl 10 1] ;		%set data for basic flux alogithm
   y=cor26(xx);
   yy(i,:)=y;
end;
zou=yy(:,4);
zoq=yy(:,6);
cdn10=yy(:,19);
cen10=yy(:,21);
hlb=yy(:,2);%account for old algorithm
taub=yy(:,3);%account for old algorithm
hsb=yy(:,1);%account for old algorithm
dter=yy(:,11);
dqer=yy(:,12)*1e3;
L=yy(:,7);
qsr=yy(:,10);



   
   %*******************  set constants  ****************
tdk=273.16;
grav=9.83;%grv(lat)
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);

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); 
js=255;je=263;
stdt=255;endt=floor(jdy(length(jdy)))+1;
epic01_Rsclear_30;%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);
hnet=0.95*rs+rlnet-hsb-hlb-hrain;

cu=.00087;bu=.00018;cx=.00085;bx=.00015;
dts=tsnk-ta-dter;
dqs=qse-qa-dqer;
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*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
     %wbar=1.61*hlb/Le/rhoa+(1+1.61*Q)*hsb/rhoa/cpa/ta;
hl_webb=yy(:,15);%rhoair.*wbar.*qa*Le/1000;


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)-psiug(zetx))./gf;
q10n=qa-qsb/.4.*(log(zq_etl/10)-psitg(zetx.*zq_etl./zu_etl));
dq10n=qse-q10n-dqer;
%ii=find(dq10n<-5);dq10n(ii)=qse(ii)-qp(ii);

figure;plot(jdy,tsnk,'o',jdy,tsg,'-',jdy,ta);set(gca,'fontsize',14);axis([stdt endt 24 32]);xlabel('Year Day (2001)');ylabel('Temperature (C)');   
stdt=253;endt=278;
figure;subplot(3,1,1);plot(jdy,U);set(gca,'fontsize',14);axis([stdt endt  0 20]);xlabel('Year Day (2001)');ylabel('Wind Speed (m/s)');   
subplot(3,1,2);plot(jdy,U.*cos(dir/180*pi),[stdt endt],[0 0]);set(gca,'fontsize',14);axis([stdt endt -12 12]);xlabel('Year Day (2001)');ylabel('U_{N} (m/s)');   
subplot(3,1,3);plot(jdy,U.*sin(dir/180*pi),[stdt endt],[0 0]);set(gca,'fontsize',14);axis([stdt endt -12 12]);xlabel('Year Day (2001)');ylabel('U_{E} (m/s)');   
bradd_f=1
if bradd_f
    ee=load('d:\data\epic\assmen_nohd.txt');
    jd_ass=ee(:,1);
    pr_ass=ee(:,2);
    ta_ass=ee(:,5);
    e_ass=ee(:,7);
    q_ass=ee(:,8);
    rh_ass=ee(:,9);
    rhoa_ass=ee(:,10);
    rhov_ass=ee(:,11);
       
end;
dqass=qsr*1e3/.4.*(log(zq_etl/zq_ass)-psitg(zq_etl./L)+psitg(zq_ass./L));

%figure;plot(jdy,qa,jd_ass,q_ass,'o');set(gca,'fontsize',14);axis([stdt endt 16 20]);xlabel('Year Day (2001)');ylabel('Specific Humidity (g/kg)');      
%figure;plot(jdy,ta,jd_ass,ta_ass,'o');set(gca,'fontsize',14);axis([stdt endt 22 30]);xlabel('Year Day (2001)');ylabel('Air Temperature (C)');      
na=length(q_ass);
j=1;
ma=[];m=[];
for i=1:na
    ii=find(jdy>jd_ass(i));
    if ~isempty(ii)
        m(j)=ii(1);
        ma(j)=i;
        j=j+1;
    end;
end;

figure;plot(jdy,qa-dqass,jd_ass,q_ass,'o');set(gca,'fontsize',14);axis([stdt endt 12 20]);xlabel('Year Day (2001)');ylabel('Specific Humidity (g/kg)');      
figure;plot(jdy,ta,jd_ass,ta_ass,'o');set(gca,'fontsize',14);axis([stdt endt 16 30]);xlabel('Year Day (2001)');ylabel('Air Temperature (C)');      
figure;plot(q_ass(ma),qa(m)-dqass(m),'o',[12 21],[12 21],q_ass(ma),q_ass(ma)+.32);xlabel('Assman q (g/kg)');ylabel('Vaisala q (g/kg)');
figure;plot(ta_ass(ma),ta(m),'o',[16 32],[16 32],ta_ass(ma),ta_ass(ma)-.23);xlabel('Assman T (C)');ylabel('Vaisala T (C)');
jj=find(qa(m)-dqass(m)>q_ass(ma)+.32);size(jj)/size(ma');
kk=find(ta(m)>ta_ass(ma)-.21);size(kk)/size(ma');


v1=u10n;
  %ii=find(abs((U-U_scs)./U)>.5);v1(ii)=NaN;
hl_webar=hl_webb;
%hlb=hlb+hl_webb;

yqs=qse;
yq=qa;
jdollars='EPIC2001';
   
   
   %ophclr=ones(np,1);ophclr=ophclr*2600;
%metz_med_tau_plt2;%_30;
%metz_med_hl_plt4;
jdq=jdy;%
epicz_daily;

stdt=255;endt=275;
jdq=jdy-5.9/24;%(convert to local time)
epicz_diurnal;

stdt=276;endt=278;
stlat=0;enlat=5;
figure;
subplot(3,1,1);plot(lat,U,'o');set(gca,'fontsize',12);axis([stlat enlat 0 10]);xlabel('Latitude(Deg)');ylabel('Wind Speed (ms^{-1})');      
subplot(3,1,2);plot(lat,tsnk,'o',lat,ta,'x');set(gca,'fontsize',12);axis([stlat enlat 16 28]);xlabel('Latitude(Deg)');ylabel('Temperature (C)');      
subplot(3,1,3);plot(lat,hsb,'o',lat,hlb,'x');set(gca,'fontsize',12);axis([stlat enlat 0 200]);xlabel('Latitude(Deg)');ylabel('Heat Flux (Wm^{-2})');      