% Dan Wolfe's updates for VOCALS 2008

% filecat('E:\CRUISES_RAW_data\GasExIII\flux\Processed\flux_5hf*.txt','E:\CRUISES_RAW_data\GasExIII\flux\flux_5hf.txt')
disp('epicz_read_rhb_2008');

clear all
close all
read_parameters_VOCALS; %Reading parameter files
tit = [cruise ' ' year];
x=load(['H:\RHB\flux\Processed\cat\flux_5hf_VOCALS2008.txt']);%flux_5hf_GasexIII_08.txt
%read file with 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) (z=15 m)
qse=x(:,8);%sea surface specific humidity,from snake (g/kg)
qa=x(:,9);%air specific humidity, PSD (g/kg) (z=15 m)
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);%ship speed, SCS gps (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, C (z=15 m)
qa_im=x(:,26);% IMET air specific humidity, g/kg (z=15 m)
s_ship=x(:,27);% IMET true wind speed, m/s (z=15 m)
dir_ship=x(:,28);% IMET true wind direction, deg (z=15 m)
psp_ship=x(:,29);%IMET solar flux, w/m^2
pir_ship=x(:,30);%IMET IR flux, w/m^2
barpress=x(:,31);%PSD atmospheric pressure, mb
RH=x(:,32);%PSD RHvais, %

zu_etl=18.5;%height of mean wind sensor, 17.7 m
zt_etl=15.;%height of mean air temperature sensor, 15 m
zq_etl=15.;%height of mean air humidity sensor, 15 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;
set(gcf,'defaultaxesfonts',14)
set(gca,'fontname','time','fontsize',14);
map_epicz_VOCALS_08;
%   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=286;je=298;% need to cHANGE???? VOCALS
%js=63;je=97;% need to cHANGE????
iv=ones(length(qa),1)*NaN;
latd=lat;lond=lon;
%stdt=62; %stratus07
stdt=294; %VOCALS
endt=floor(jdy(length(jdy)))+1; %STRATUS started Oct 14, 2006
%stdt=306;endt=floor(jdy(length(jdy)))+1;
jdys=jdy;
epic01_Rsclear_5x;%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);

timelab=['Year Day ' year];
sprin=[way_images_flux cruise year(3:4)];
figure;
plot(jdy,tsnk,'b.',jdy,ta,'g');
set(gca,'fontsize',14);
axis([stdt endt 15 35]);
xlabel(timelab);
ylabel('Temperature (C)');
title(tit)
legend('SST','Air temp');
print_buffer = [sprin '_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 20]);
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 -20 20]);
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 -20 20]);
xlabel(timelab);
ylabel('U_{E} Spd (ms^{-1})');
print_buffer = [sprin '_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 20]);
ylabel('True 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('True Wind Dir (deg)');
xlabel(timelab);
print_buffer = [sprin '_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=[cruise year];


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

%stdt=255;endt=275;

% stdt=294;endt=floor(max(jdy))+1;%216;  % GOMECC07 Jul 9 - Aug 4
jdq=jdy;%-85/15/24;%(convert to local time)

epicz_weller_diurnal_08;
%%%%%%%%%%%%%%%%%%%%%%%%
%stdt=276;endt=278;
%Stratus
%stlon=-85;enlon=-70;
% AMMA07

stlon=-85;enlon=-70;
ii=find(jdy>294& 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('Speed (ms^{-1})');

subplot(3,1,2);
plot(lon(ii),tsnk(ii),'k.',lon(ii),ta(ii),'r.');
set(gca,'fontsize',10);
axis([stlon enlon 15 35]);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 300]);xlabel('Long(Deg)');ylabel('Heat Flux (Wm^{-2})');
legend('Hsen','Hlat','location','Northwest')
print_buffer = [sprin '_U_SSTT_Net_lon.jpg'];
print('-djpeg90 ', print_buffer);

stlat=-20;enlat=4;
ii=find(jdy>285& jdy<299);
figure;
subplot(3,1,1);plot(lat(ii),U(ii),'k.');
title(tit)
set(gca,'fontsize',12);axis([stlat enlat 0 20]);ylabel('Speed (ms^{-1})');

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

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 Flux (Wm^{-2})');
set(gca,'fontsize',8);legend('Hsen','Hlat','location','Best')
print_buffer = [sprin '_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(-60:.25:-35,-55:.25:-48);%GASEX III site
    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;

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

   % yass=load('E:\CRUISES_RAW_data\AMMA06\AMMA_2006\RHB\AssmanPsychrometer\ASSMANS_06_nohd_leg1.txt');
     yass=load([way_proc_data_flux 'ST8_assmans_nohds.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=[];assjk=length(jdass);psp_a=[];
    for i=1:length(jdass)
        ias=find(jdy>jdass(i));
        if ~isempty(ias)
            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);
            psp_a(i)=psp(ias(1)-1);
        else
            jk=i;
        end;
    end;
        %Psychometer correction to ETl and ship instrument height
        dltq=qst/.4.*(log(15/7)-psit_26(15./LL)+psit_26(7./LL));
        dltt=tst/.4.*(log(15/7)-psit_26(15./LL)+psit_26(7./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',[6 12],[6 12],'--');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',[13 20],[13 20],'--');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)

    
    [nanmean([ta_et' ta_sp' ta_ass+dltt']) std([ ta_et'-ta_sp' ta_et'-(ta_ass+dltt') ])] %Tair (C)
    [nanmean([qa_et' qa_sp' qa_ass+dltq']) std([ qa_et'-qa_sp' qa_et'-(qa_ass+dltq') ])] %Qair (g/kg)
     ii=find(psp_a<100);
    [nanmean([ta_et(ii)' ta_sp(ii)' ta_ass(ii)+dltt(ii)']) std([ ta_et(ii)'-ta_sp(ii)' ta_et(ii)'-(ta_ass(ii)+dltt(ii)') ])] %Tair (C)
end  % loop around Assman

%  stdt=263;  
%  endt=floor(jdy(length(jdy)))+1;

figure;plot(jdy,tsnk,'k.',jdy,tsg);
set(gca,'fontsize',14);
xlabel(timelab);
ylabel('Temperature (C)');
legend('PSD Sea Snake','SCS 5m TSG')
axis([stdt endt 15 35]);
title(tit)
print_buffer = [sprin '_SST_PSDSCS.jpg'];
print('-djpeg90 ', print_buffer);

%figure;plot(jdy,ta,'b',jdy,ta_im,jdass,ta_ass+dltt','o');
figure;plot(jdy,ta,'b');  % DEW 5/18/07
set(gca,'fontsize',14);
xlabel(timelab);
ylabel('Air Temperature (C)');
axis([stdt endt 15 35]);
title(tit)
print_buffer = [sprin '_Tair.jpg'];
print('-djpeg90 ', print_buffer);

%figure;plot(jdy,qa,'b',jdy,qa_im,jdass,qa_ass+dltq','o');
figure;plot(jdy,qa,'b');  % DEW 5/18.07
set(gca,'fontsize',14);
xlabel(timelab);
ylabel('Specific Humidty (g/kg)');
axis([stdt endt 5 20]);
title(tit)
print_buffer = [sprin '_SpecRH.jpg'];
print('-djpeg90 ', print_buffer);

% % ll=load('coast.mat');
% % figure;plot(ll.long,ll.lat);
% % axis equal
% % axis([-95 -70 -25 10]);
% % hold on;
% % scatter(lon,lat,5,jdy);
% % dj=floor(jdy(length(jdy))-jdy(1));
% % jdx=jdy(1);
% % for i=1:dj;
% %     ii=find(floor(jdy)==jdx);
% %     text(lon(ii(1)),lat(ii(1)),num2str(jdx));
% %     jdx=jdx+1;
% % end;


figure;
set(gcf,'defaultaxesfonts',12)
set(gca,'fontname','time','fontsize',12);
subplot(6,1,1);plot(jdy,U,'b',[68 68],[0 20],'k--',[73 73],[0 20],'k--',[81 81],[0 20],'k--',[96 96],[0 20],'k--');xlim([stdt endt]);ylabel('WS');title('5min Surface Meteorological Data 10/12/08 - 10/24/08');
subplot(6,1,2);plot(jdy,dir,'b',[68 68],[0 400],'k--',[73 73],[0 400],'k--',[81 81],[0 400],'k--',[96 96],[0 400],'k--');xlim([stdt endt]);ylabel('WD');
subplot(6,1,3);plot(jdy,ta,'b',jdy,tsnk,'r',[68 68],[15 35],'k--',[73 73],[15 35],'k--',[81 81],[15 35],'k--',[96 96],[15 35],'k--');axis([stdt endt 15 35]);ylabel('Ta, Ts (red)');
subplot(6,1,4);plot(jdy,psp,'b',jdy,rl,'r',[68 68],[-5 1200],'k--',[73 73],[-5 1200],'k--',[81 81],[-5 1200],'k--',[96 96],[-5 1200],'k--');axis([stdt endt -5 1200]);ylabel('SW , LW (red)');
subplot(6,1,5);plot(jdy,RH,'b',[68 68],[40 100],'k--',[73 73],[40 100],'k--',[81 81],[40 100],'k--',[96 96],[40 100],'k--');axis([stdt endt 40 100]);ylabel('RH');
subplot(6,1,6);plot(jdy,barpress,'b',[68 68],[900 1050],'k--',[73 73],[900 1050],'k--',[81 81],[900 1050],'k--',[96 96],[900 1050],'k--');axis([stdt endt 1000 1030]);ylabel('BP');
xlabel('2008 Yearday');

figure;
set(gcf,'defaultaxesfonts',12)
set(gca,'fontname','time','fontsize',12);
subplot(6,1,1);plot(jdy,taub,'b',[68 68],[-.05 1],'k--',[73 73],[-.05 1],'k--',[81 81],[-.05 1],'k--',[96 96],[-.05 1],'k--');ylabel('Stress');title('5min Flux 10/12/08 - 10/24/08');axis([stdt endt -.05 1])
subplot(6,1,2);plot(jdy,hnet,'b',[68 68],[-500 800],'k--',[73 73],[-500 800],'k--',[81 81],[-500 800],'k--',[96 96],[-500 800],'k--');ylabel('Hnet');axis([stdt endt -500 800])
subplot(6,1,3);plot(jdy,0.95*rs,'b',[68 68],[-100 800],'k--',[73 73],[-100 800],'k--',[81 81],[-100 800],'k--',[96 96],[-100 800],'k--');ylabel('Rs');axis([stdt endt -10 1200])
subplot(6,1,4);plot(jdy,rlnet,'b',[68 68],[-150 100],'k--',[73 73],[-150 100],'k--',[81 81],[-150 100],'k--',[96 96],[-150 100],'k--');ylabel('Rnl');axis([stdt endt -150 100])
subplot(6,1,5);plot(jdy,-hsb,'b',[68 68],[-200 100],'k--',[73 73],[-200 100],'k--',[81 81],[-200 100],'k--',[96 96],[-200 100],'k--');ylabel('Hsb');axis([stdt endt -200 100])
subplot(6,1,6);plot(jdy,-hlb,'b',[68 68],[-250 100],'k--',[73 73],[-250 100],'k--',[81 81],[-250 100],'k--',[96 96],[-250 100],'k--');ylabel('Hlb');axis([stdt endt -250 100])
xlabel('2008 Yearday');

[mean(U(~isnan(U))) std(U(~isnan(U))) max(U(~isnan(U)))  min(U(~isnan(U)))]
[mean(ta(~isnan(ta))) std(ta(~isnan(ta))) max(ta(~isnan(ta))) min(ta(~isnan(ta)))]
[mean(tsnk(~isnan(tsnk))) std(tsnk(~isnan(tsnk))) max(tsnk(~isnan(tsnk))) min(tsnk(~isnan(tsnk)))]
tt=find(RH>110);RH(tt)=RH(tt-1);
[mean(RH(~isnan(RH))) std(RH(~isnan(RH))) max(RH(~isnan(RH))) min(RH(~isnan(RH)))]
tt=find(barpress>1025);barpress(tt)=barpress(tt-1);
tt=find(barpress<970);barpress(tt)=barpress(tt-1);
[mean(barpress(~isnan(barpress))) std(barpress(~isnan(barpress))) max(barpress(~isnan(barpress))) min(barpress(~isnan(barpress)))]

[mean(taub(~isnan(taub))) std(taub(~isnan(taub))) max(taub(~isnan(taub))) min(taub(~isnan(taub)))]
[mean(hnet(~isnan(hnet))) std(hnet(~isnan(hnet))) max(hnet(~isnan(hnet))) min(hnet(~isnan(hnet)))]
rlud=0.95*rs;
[mean(rlud(~isnan(rlud))) std(rlud(~isnan(rlud))) max(rlud(~isnan(rlud))) min(rlud(~isnan(rlud)))]
tt=find(rlnet>50);rlnet(tt)=rlnet(tt-1);
[mean(rlnet(~isnan(rlnet))) std(rlnet(~isnan(rlnet))) max(rlnet(~isnan(rlnet))) min(rlnet(~isnan(rlnet)))]
[mean(-hsb(~isnan(-hsb))) std(-hsb(~isnan(-hsb))) max(-hsb(~isnan(-hsb))) min(-hsb(~isnan(-hsb)))]
[mean(-hlb(~isnan(-hlb))) std(-hlb(~isnan(-hlb))) max(-hlb(~isnan(-hlb))) min(-hlb(~isnan(-hlb)))]

%a=85;b=86;
a=294;b=305;
figure;
set(gcf,'defaultaxesfonts',14)
set(gca,'fontname','time','fontsize',14);
subplot(4,1,1);plot(jdy,org1,'b');ylabel('Rain(mm/hr)');xlim([a b]);ylim([0 2]);
subplot(4,1,2);plot(jdy,ta,'b',jdy,tsnk,'g.');ylabel('T (C)');xlim([a b]);ylim([15 35]);
subplot(4,1,3);plot(jdy,U,'b');ylabel('U (m/s)');xlim([a b]);
subplot(4,1,4);plot(jdy,RH,'b');ylabel('RH (%)');xlim([a b]);ylim([0 100]);
xlabel('2008 Yearday');

figure;
set(gcf,'defaultaxesfonts',14)
set(gca,'fontname','time','fontsize',14);
subplot(4,1,1);plot(jdy,0.95*rs,'b');ylabel('Rs (W/m^2)');axis([a b 0 1200]);
subplot(4,1,2);plot(jdy,rl,'b');ylabel('Rl (W/m^2)');xlim([a b]);
subplot(4,1,3);plot(jdy,-hsb,'b',jdy,-hlb,'g.',jdy,rlnet,'r.',jdy,-hrain,'c.',[a b],[0 0],'m--');ylabel('H (W/m^2)');axis([a b -110 100]);
subplot(4,1,4);plot(jdy,hnet,'b',[a b],[0 0],'m--');ylabel('Hnet (W/m^2)');xlim([a b]);
xlabel('2008 Yearday');


figure;
set(gcf,'defaultaxesfonts',14)
set(gca,'fontname','time','fontsize',14);
plot(jdy,U,'-',jdy,s_ship,'-');xlim([a b])
xlabel('Julian Day (2008)');
legend('PSD','IMET')
ylabel('Wind speed (m/s)')
title('True wind speed')
