disp('epicz_read_rhb_2007');

clear all
close all
tit = 'GOMECC 2007';
%x=load('E:\CRUISES_RAW_data\AMMA06\AMMA_2006\RHB\flux\Processed_new\flux_5hf_amma_06.txt');%read file with hour average data; set your local directory
x=load('S:\GOMECC_2007\RHB\flux\Processed\flux_5hf_AMMA07.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,PSD sonic (m/s)
dir=x(:,3);%true wind direction, PSD sonic (deg)
tsnk=x(:,4);%sea snake temperature,PSD, 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, PSD (C)
qse=x(:,8);%sea surface specific humidity,from snake (g/kg)
qa=x(:,9);%air specific humidity, PSD (g/kg)
psp=x(:,10)*1.01;%downward solar flux, PSD units (W/m^2)
rl=x(:,11);%downward IR flux,  PSD units (W/m^2)
org1=x(:,12);%rainrate, PSD 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, PSD (m/s)
reldir=x(:,16);%relative wind direction (from),clockwise rel ship's bow, PSD sonic (deg)
lat=x(:,17);%latitude, deg  (SCS pcode)
lon=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)
ta_im=x(:,25);%IMET air temp @21 m
qa_im=x(:,26);%IMET
s_ship=x(:,27);%imet
dir_ship=x(:,28);%IMET
psp_ship=x(:,29);%imet

zu_etl=18.5;%height of mean wind sensor, 17.7 m
zt_etl=15.;%height of mean air temperature sensor, 15.5 m
zq_etl=15.;%height of mean air humidity sensor, 15.5 m
%*****   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)

org=zeros(length(org1),1);
rs=psp;
%ii=find(psp_ship>-10);
%rsm=mean([psp psp_ship]')';
%rs(ii)=rsm(ii);
latd=lat;
lond=lon;
figure;
map_amma2007;  %DEW 5/19/2007 skip map for testing
%   lat=latdecode(latx);
%  lon=-latdecode(lonx);%negative for W long
yy=[];ygs=[];xvc=[];
tsea=tsnk;
%     ii=find(jdy<340.83);tsea(ii)=tsg(ii);qse(ii)=qsea(tsg(ii));
for i=1:length(x);
    xx=[U(i) 0 tsea(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_5(xx);
    yy(i,:)=y;
    y2=cor26_co2_oz_1(xx);
    ygs(i,:)=y2;
end;
zou=yy(:,4);
zoq=yy(:,6);
ust=yy(:,8);
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);
tsr=yy(:,9);

U10_et=U+ust/.4.*(log(10./zu_etl)-psit_26(10./L)+psit_26(zu_etl./L));
U10_im=s_ship+ust/.4.*(log(10./15)-psit_26(10./L)+psit_26(15./L));

k_co2=ygs(:,22);%raw xfer velocity, F=Sol k * DPCo2  k in m/s
V_co2=ygs(:,23);%Flux=V*DPco2, V in mole/m^2/sec/atm
k_oz=ygs(:,24);%raw xfer velocity, ozone
rsa=ygs(:,25);%atmospheric resistance
k_atmos=ust./rsa;
resoz=300;%removal resistance on the ocean side, s/m
vdepoz=(resoz+1./k_oz).^(-1);
%xvc=[jdy k_co2 V_co2 vdepoz];size(xvc);
%*******************  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=190;je=192;% need to cHANGE????
js=289;je=293;% need to cHANGE????
iv=ones(length(qa),1)*NaN;
latd=lat;lond=lon;
stdt=123; %AMMA07
endt=floor(jdy(length(jdy)))+1; %STRATUS started Oct 14, 2006
%stdt=306;endt=floor(jdy(length(jdy)))+1;
epic01_Rsclear_5;%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);  %commented out DEW 21 Oct 2006
zetx=min(zetx,50);   %commented out DEW 21 Oct 2006
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;
%ii=find(dq10n<-5);dq10n(ii)=qse(ii)-qp(ii);

jj=find(jdy>168.55 & jdy<174.7);ta(jj)=NaN;
figure;
plot(jdy,tsnk,'b.',jdy,ta,'g');
set(gca,'fontsize',14);
axis([stdt endt 20 30]);
xlabel('Year Day (2007)');
ylabel('Temperature (C)');
title(tit)
legend('SST','Air temp');
print_buffer = ['S:\GOMECC_2007\RHB\flux\Processed_Images\AMMA07_FLUX_SST_Tair.jpg'];
print('-djpeg90 ', print_buffer);

%print d:\docs\reports\figs\figxx_temps_all.jpg -djpeg
%stdt=315;endt=278;

figure;
subplot(3,1,1);
plot(jdy,U);
title(tit)
set(gca,'fontsize',14);
axis([stdt endt  0 12]);
ylabel('Wind Spd (ms^{-1})');

subplot(3,1,2);
plot(jdy,U.*cos(dir/180*pi),[stdt endt],[0 0]);
set(gca,'fontsize',14);
axis([stdt endt -12 12]);
ylabel('U_{N} Spd (ms^{-1})');

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 (2007)');
ylabel('U_{E} Spd (ms^{-1})');
print_buffer = ['S:\GOMECC_2007\RHB\flux\Processed_Images\AMMA07_FLUX_Sonic_U3.jpg'];
print('-djpeg90 ', print_buffer);

figure;
subplot(2,1,1);
plot(jdy,U);
title(tit)
set(gca,'fontsize',14);
axis([stdt endt  0 12]);
ylabel('Wind Spd (ms^{-1})');

subplot(2,1,2);
plot(jdy,dir,[stdt endt],[0 0]);
set(gca,'fontsize',14);
axis([stdt endt 0 360]);
ylabel('Wind Dir (T)');
print_buffer = ['S:\GOMECC_2007\RHB\flux\Processed_Images\AMMA07_FLUX_Sonic_SpdDir.jpg'];
print('-djpeg90 ', print_buffer);
stabeff=1./(1-sqrt(cdn10)/.4.*psiu_26(zetx*zu_etl/10)).^2;

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='AMMA2007';


%ophclr=ones(np,1);ophclr=ophclr*2600;
%metz_med_tau_plt2;%_30;
%metz_med_hl_plt4;
org=max(0,org1-.06);
jdq=jdy;%

if 0   % skip processing for testing DEW 5/19/2007
%epicz_weller_daily_07;
end

%stdt=255;endt=275;
if 1   % skip processing for testing DEW 5/19/2007
stdt=123;endt=150;  % AMMA07 May 3-30
jdq=jdy;%-85/15/24;%(convert to local time)

epicz_weller_diurnal_07;
end
%stdt=276;endt=278;
%Stratus
%stlon=-85;enlon=-70;
% AMMA07
stlon=-65;enlon=-20;

% ii=find(jdy>340 & jdy<353);
% ii=find(jdy>278 & jdy<293);
ii=find(jdy>123);
%ii=find(jdy>305);

figure;
subplot(3,1,1);
plot(lon(ii),U(ii),'k.');
title(tit)
set(gca,'fontsize',10);
axis([stlon enlon 0 20]);
ylabel('Spd (ms^{-1})');

subplot(3,1,2);
plot(lon(ii),tsnk(ii),'k.',lon(ii),ta(ii),'r.');
set(gca,'fontsize',10);
axis([stlon enlon 17 30]);ylabel('Temp (C)');
legend('SST','Ta','location','Northwest')

subplot(3,1,3);
plot(lon(ii),hsb(ii),'k.',lon(ii),hlb(ii),'r.');
set(gca,'fontsize',10);
axis([stlon enlon -20 250]);xlabel('Long(Deg)');ylabel('Heat Flx (Wm^{-2})');
legend('Hsen','Hlat','location','Northwest')
print_buffer = ['S:\GOMECC_2007\RHB\flux\Processed_Images\AMMA07_FLUX_U_SSTT_Net_lon.jpg'];
print('-djpeg90 ', print_buffer);

stlat=0;enlat=20;

figure;
subplot(3,1,1);plot(lat(ii),U(ii),'k.');
title(tit)
set(gca,'fontsize',12);axis([stlat enlat 0 15]);ylabel('Spd (ms^{-1})');

subplot(3,1,2);plot(lat(ii),tsnk(ii),'k.',lat(ii),ta(ii),'r.');
set(gca,'fontsize',12);axis([stlat enlat 17 30]);ylabel('Temp (C)');
legend('SST','Ta','location','Northwest')

subplot(3,1,3);plot(lat(ii),hsb(ii),'k.',lat(ii),hlb(ii),'r.');
set(gca,'fontsize',12);axis([stlat enlat -20 250]);xlabel('Lat(Deg)');ylabel('Heat Flx (Wm^{-2})');
legend('Hsen','Hlat','location','Northwest')
print_buffer = ['S:\GOMECC_2007\RHB\flux\Processed_Images\AMMA07_FLUX_U_SSTT_Net_lat.jpg'];
print('-djpeg90 ', print_buffer);

skipit=0;
if skipit
    %[xi,yi]=meshgrid(-71:.25:-66,41.5:.25:44.5);
    [xi,yi]=meshgrid(-115:.25:-85,-15:.25:15);%TAO
    ZI = griddata(lon,lat,U,xi,yi);
    figure;[c,h] = contour(xi,yi,ZI); clabel(c,h), colorbar;
    title('Wind Speed (m/s)');

    ZI = griddata(lon,lat,100*vdepoz,xi,yi);
    figure;[c,h] = contour(xi,yi,ZI); clabel(c,h), colorbar;
    title('Ozone Deposition Velocity (cm/s)');

    ZI = griddata(lon,lat,ta-tsnk,xi,yi);
    figure;[c,h] = contour(xi,yi,ZI); clabel(c,h), colorbar;
    title('Air-Sea T difference (C)');

    ZI = griddata(lon,lat,tsnk,xi,yi);
    figure;[c,h] = contour(xi,yi,ZI); clabel(c,h), colorbar;
    title('Sea T C)');


    ZI = griddata(lon,lat,stabeff,xi,yi);
    figure;[c,h] = contour(xi,yi,ZI); clabel(c,h), colorbar;
    title('Stability Effect on Fluxes at 10-m');


    %print d:\docs\reports\figs\figxx_long.jpg -djpeg
end;


if 0
    %ASSMAN intercomparison (Where does the Assman come from?)

    yass=load('E:\CRUISES_RAW_data\AMMA06\AMMA_2006\RHB\AssmanPsychrometer\ASSMANS_06_nohd_leg1.txt');
    hhmm=yass(:,2);   %hour minutes
    jdass=yass(:,1)+floor(hhmm/100)/24+(hhmm-100*floor(hhmm/100))/24/60;   %Julian Day Assman
    p_ass=yass(:,3);  %pressure in hPa or mbar
    ta_ass=yass(:,4);  % dry bulb temperature in degC
    rh_ass=yass(:,8);  %relative humidity in %
    qa_ass=qair3([p_ass ta_ass rh_ass]); %Assman Psychometer specific humidity in g/kg
    qa_et=[];qa_sp=[];ta_et=[];ta_sp=[];
    for i=1:length(jdass)
        ias=find(jdy>jdass(i));
        qa_et(i)=qa(ias(1)-1);
        qa_sp(i)=qa_im(ias(1)-1);
        ta_et(i)=ta(ias(1)-1);
        ta_sp(i)=ta_im(ias(1)-1);
        LL(i)=L(ias(1)-1);
        qst(i)=1e3*qsr(ias(1)-1);
        tst(i)=tsr(ias(1)-1);

        %Psychometer correction to ETl and ship instrument height
        dltq=qst/.4.*(log(15/9)-psit_26(15./LL)+psit_26(9./LL));
        dltt=tst/.4.*(log(15/9)-psit_26(15./LL)+psit_26(9./LL));

        ii=find(abs((ta_ass-ta_et')>1));ta_ass(ii)=ta_ass(ii)-1;

        [ppa,ssq]=polyfit(qa_et,qa_ass'+dltq,1);
        [popq,del2a]=polyval(ppa,qa_et,ssq);

        [pps,ssq]=polyfit(qa_et,qa_sp,1);
        [popq,del2s]=polyval(pps,qa_et,ssq);

        figure;plot(qa_et,qa_sp,'o',qa_et,qa_ass'+dltq,'x',[7 20],[7 20],'--');set(gca,'fontsize',12);xlabel('PSD q (g/kg)');ylabel('IMET and Assman q (g/kg)');
        title('Assman psychometer (x) and IMET ship (o) versus PSD');axis([13 20 13 20]);
        %[Ps,Ss]=polyfit(qa_et,qa_sp,1)
        %[Pa,Sa]=polyfit(qa_et,qa_ass',1)

        figure;plot(ta_et,ta_sp,'o',ta_et,ta_ass'+dltt,'x',[16 30],[16 30],'--');set(gca,'fontsize',12);xlabel('PSD Ta (C)');ylabel('IMET and Assman Ta (C)');
        title('Assman psychometer (x) and IMET ship (o) versus PSD');axis([20 30 20 30]);
        %[Ps,Ss]=polyfit(ta_et,ta_sp,1)
        %[Pa,Sa]=polyfit(ta_et,ta_ass',1)

    end % skip Assman section

    [mean([ta_et' ta_sp' ta_ass+dltt']) std([ ta_et'-ta_sp' ta_et'-ta_ass+dltt' ])] %Tair (C)
    [mean([qa_et' qa_sp' qa_ass+dltq']) std([ qa_et'-qa_sp' qa_et'-qa_ass+dltq' ])] %Qair (g/kg)
    stdt=123;  %AMMA07
    endt=floor(jdy(length(jdy)))+1;
    
end  % loop around Assman

figure;plot(jdy,tsnk,'.',jdy,tsg);
set(gca,'fontsize',14);
xlabel('Year Day (2007)');
ylabel('Temperature (C)');
legend('PSD Sea Snake','SCS 5m TSG')
axis([stdt endt 20 30]);
title(tit)
print_buffer = ['S:\GOMECC_2007\RHB\flux\Processed_Images\AMMA07_SST_PSDSCS.jpg'];
print('-djpeg90 ', print_buffer);

%figure;plot(jdy,ta,'.',jdy,ta_im,jdass,ta_ass+dltt','o');
figure;plot(jdy,ta,'.');  % DEW 5/18/07
set(gca,'fontsize',14);
xlabel('Year Day (2007)');
ylabel('Air Temperature (C)');
axis([stdt endt 15 30]);
title(tit)
print_buffer = ['S:\GOMECC_2007\RHB\flux\Processed_Images\AMMA07_FLUX_Tair.jpg'];
print('-djpeg90 ', print_buffer);

%figure;plot(jdy,qa,'.',jdy,qa_im,jdass,qa_ass+dltq','o');
figure;plot(jdy,qa,'.');  % DEW 5/18.07
set(gca,'fontsize',14);
xlabel('Year Day (2007)');
ylabel('Specific Humidty (g/kg)');
axis([stdt endt 10 25]);
title(tit)
print_buffer = ['S:\GOMECC_2007\RHB\flux\Processed_Images\AMMA07_FLUX_SpecRH.jpg'];
print('-djpeg90 ', print_buffer);
