% epicz_read
% Xcruise 2009 :: 2009-09-07 :: LB
% Reads flux_5hf_Xcruise_2009 and does plots and other stuff

clear all
close all
fclose('all');

dv = clock;  % gets current date and time as date vector.
Vdate=datevec(datenum(dv(1),dv(2),dv(3)-1));
    
%%%%%%%%%%%  Read parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
currentMname='DYNAMO_2011';   %Gets the name of currently running M-file to pick 
% Fluxroot='D:\STRATUS_2011\RHB\Scientific_analysis\programs\';  % name of the path where the 'flux' folder + read_parameters file are installed.
Fluxroot='C:\Data\CRUISES\DYNAMO_2011\Revelle\Scientific_analysis\programs\on DAS\flux\';  % name of the path where the 'flux' folder + read_parameters file are installed.

cd([Fluxroot, 'flux'])         %Changes current directory to fluxroot/flux
addpath(genpath(Fluxroot))     %adds flux folder and all of its subdirectories to the top of the MATLAB search path. 

%Gets the last uploaded parameter file from Fluxroot directory
read_parameters_DYNAMO_2011;

%%%%%%%%%%%  Defines some options  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jwave_default=0; % off by default
hwave_default=2; % wave height (m)
twave_default=10; % wave period (s)
ass=false; % compare to Assman psychrometry
graphformat='.png';  %select graphics format files
graphdevice='-dpng'; %select graphic device
prtit=true;  % for printing plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(['epicz_read: ' cruise]);

% concatenate flux_Xhf files if they haven't been already
if ~exist([way_proc_data_flux cruise '_PSD_flux_5min_all.txt'],'file')
    filecat([way_proc_data_flux cruise '_PSD_flux_5min_*.txt'],[way_proc_data_flux cruise '_PSD_flux_5min_all.txt']);
end
if ~exist([way_proc_data_flux cruise '_PSD_flux_30min_all.txt'],'file')
    filecat([way_proc_data_flux cruise '_PSD_flux_30min_*.txt'],[way_proc_data_flux cruise '_PSD_flux_30min_all.txt']);
end

x=load([way_proc_data_flux cruise '_PSD_flux_5min_all.txt']);
%read file with average data; set your local directory
% jdy=x(:,1);%julian day at beginning of time average
jdy=x(:,1)+datenum(Vdate(1),1,0);%julian day to Matlab format
% jdy=jdy+datenum(Vdate(1),1,0);  %converts time to Serial Date Number
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)
zsst=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
press=x(:,31);%PSD atmospheric pressure, mb
RH=x(:,32);%PSD RHvais, %

zu_etl=zus;%height of mean wind sensor, 17.7 m on RHB
zt_etl=zts;%height of mean air temperature sensor, 15 m  on RHB
zq_etl=zqs;%height of mean air humidity sensor, 15 m  on RHB
%*****   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)

stdt=floor(jdy(1));
endt=ceil(jdy(end)); 
xtickspace=2;  %defines spacing for plot on x-axis

% %plots for Jennifer McAdoo, 09/04 to 09/10
% % Uses SCS wind direction as 30deg angle correction was only performed on 09/09/2010 
% stdt=247;
% endt=254; 
% dir=dir_ship;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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);
title(sprintf('%s.  R/V %s track',cruise,ship),'FontWeight','Bold','Interpreter','none')
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
% map
% eval(['map_', currentMname,'_ep(lon,lat,cruise,ship,[Lonmin Lonmax Latmin Latmax],Fluxroot,jdy)'])
eval(['map_', currentMname,'(lon,lat,cruise,ship,[Lonmin Lonmax Latmin Latmax],Fluxroot)'])
if prtit; print(graphdevice,[way_images_flux cruise '_ship_track_all' graphformat]); end;

       
%   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  press(i) zu_etl zt_etl zq_etl lat(i) 1 0 twave_default hwave_default]; %set data for basic flux algorithm
    y=cor30b(xx); % cor30a or b requires wave height for roughness model
    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);

iv=ones(length(qa),1)*NaN;
latd=lat;
lond=lon;
jdys=jdy;
%epic01_Rsclear_5x;%get clear sky solar flux data
eval([currentMname,'_Rsclear_1' ])

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);


figure;
plot(jdy,tsnk,'b.',jdy,ta,'g');
set(gca,'fontsize',12);
axis([stdt endt 26 32]);
% xlim([stdt endt]);
xlabel('Day (UTC)');
ylabel('Temperature (degC)');grid on;%grid minor
% title(sprintf('%s.  PSD air and SCS sea temperatures',cruise),'FontWeight','Bold','Interpreter','none')
title(sprintf('%s.  PSD air and sea temperatures',cruise),'FontWeight','Bold','Interpreter','none')
set(gca(gcf),'XTick',stdt:xtickspace:endt);datetick('x','mm/dd','keepticks','keeplimits');
legend('Sea','Air','Location','Best');
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
% if prtit; print(graphdevice,[way_images_flux cruise '_PSD-SCS_Tair_comparison_all' graphformat]); end;
if prtit; print(graphdevice,[way_images_flux cruise '_PSD-SST_Tair_comparison' graphformat]); end;

figure;
subplot(3,1,1);set(gca,'fontsize',12);
plot(jdy,U);
title(sprintf('%s.  PSD True wind speed',cruise),'FontWeight','Bold','Interpreter','none')
axis([stdt endt  0 20]);
grid on;%grid minor
set(gca(gcf),'XTick',stdt:xtickspace:endt);datetick('x','mm/dd','keepticks','keeplimits');
ylabel('Wind speed (ms^{-1})');

subplot(3,1,2);
plot(jdy,U.*cos(dir/180*pi),[stdt endt],[0 0]);
set(gca,'fontsize',12);
axis([stdt endt -20 20]);
grid on;%grid minor
set(gca(gcf),'XTick',stdt:xtickspace:endt);datetick('x','mm/dd','keepticks','keeplimits');
ylabel('U_{N} speed (ms^{-1})');

subplot(3,1,3);
plot(jdy,U.*sin(dir/180*pi),[stdt endt],[0 0]);
set(gca,'fontsize',12);
axis([stdt endt -20 20]);
grid on;%grid minor
set(gca(gcf),'XTick',stdt:xtickspace:endt);datetick('x','mm/dd','keepticks','keeplimits');
xlabel('Day (UTC)');
ylabel('U_{E} speed (ms^{-1})');
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
if prtit; print(graphdevice,[way_images_flux cruise '_PSD_Sonic_speeds' graphformat]); end;

figure;
subplot(2,1,1);set(gca,'fontsize',14);
plot(jdy,U);
title(sprintf('%s.  PSD True wind speed and direction',cruise),'FontWeight','Bold','Interpreter','none')
grid on;%grid minor
set(gca(gcf),'XTick',stdt:xtickspace:endt);datetick('x','mm/dd','keepticks','keeplimits');
axis([stdt endt  0 15]);
ylabel('Speed (ms^{-1})');

subplot(2,1,2);
plot(jdy,dir,'.',[stdt endt],[0 0]);
set(gca,'fontsize',14);
axis([stdt endt 0 360]);
grid on;%grid minor
set(gca(gcf),'XTick',stdt:xtickspace:endt);datetick('x','mm/dd','keepticks','keeplimits');
ylabel('Direction (deg)');
xlabel('Day (UTC)');
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
if prtit; print(graphdevice,[way_images_flux cruise '_PSD_Sonic_SpdDir' graphformat]); end;



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 yyyy];

org=max(0,org1); %-.06); % adjust the optical rain gauge--do in timatch

jdq=jdy;%
%epicz_daily;
eval([ 'epicz_daily_', currentMname])

jdq=jdy;%-85/15/24;%(convert to local time)
%epicz_diurnal;
eval([ 'epicz_diurnal_', currentMname])

stlon=Lonmin;enlon=Lonmax;
stlon=78;enlon=90;

ii=find(jdy>=stdt & jdy<endt);

figure;
subplot(3,1,1);
plot(lon(ii),U(ii),'k.');
title(sprintf('%s',cruise),'FontWeight','Bold','Interpreter','none')
set(gca,'fontsize',10);
axis([stlon enlon 0 20]);
grid on;%grid minor
ylabel('True wind 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 20 35]);ylabel('Temperature (degC)');
grid on;%grid minor
legend('SST','Ta','location','Best')

subplot(3,1,3);
plot(lon(ii),hsb(ii),'k.',lon(ii),hlb(ii),'r.');
set(gca,'fontsize',10);
axis([stlon enlon -20 600]);xlabel('Long (deg)');ylabel('Heat Flux (Wm^{-2})');
legend('Hsen','Hlat','location','Best')
grid on;%grid minor
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
if prtit; print(graphdevice,[way_images_flux cruise '_U_SSTT_Net_lon' graphformat]); end;


% stlat=Latmin;enlat=Latmax;
stlat=-3;enlat=10;

figure;
subplot(3,1,1);plot(lat(ii),U(ii),'k.');
title(sprintf('%s',cruise),'FontWeight','Bold','Interpreter','none')
set(gca,'fontsize',12);axis([stlat enlat 0 20]);ylabel('Speed (ms^{-1})');
grid on;%grid minor

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

subplot(3,1,3);plot(lat(ii),hsb(ii),'k.',lat(ii),hlb(ii),'r.');
set(gca,'fontsize',12);axis([stlat enlat -20 600]);xlabel('Lat (deg)');ylabel('Heat Flux (Wm^{-2})');
set(gca,'fontsize',8);legend('Hsen','Hlat','location','Best')
grid on;%grid minor
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
if prtit; print(graphdevice,[way_images_flux cruise '_U_SSTT_Net_lat' graphformat]); end;


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');
end;

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;
% set(gca,'fontsize',14);
% plot(jdy,tsnk,'k.',jdy,tsg,'.');
% xlabel('Day (UTC)');grid on;%grid minor
% ylabel('Temperature (degC)');
% legend('PSD Sea Snake','SCS 2m TSG','location','Best')
% axis([stdt endt 20 35]);
% set(gca(gcf),'XTick',stdt:xtickspace:endt);datetick('x','mm/dd','keepticks','keeplimits');
% title(sprintf('%s',cruise),'FontWeight','Bold','Interpreter','none')
% annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
% if prtit; print(graphdevice,[way_images_flux cruise '_SST_PSDSCS' graphformat]); end;
% 

%figure;plot(jdy,ta,'b',jdy,ta_im,jdass,ta_ass+dltt','o');
% figure;
% set(gca,'fontsize',14);
% plot(jdy,ta,'b');  % DEW 5/18/07
% xlabel('Day (UTC)');grid on;%grid minor
% ylabel('Air Temperature (degC)');
% axis([stdt endt 20 35]);
% set(gca(gcf),'XTick',stdt:xtickspace:endt);datetick('x','mm/dd','keepticks','keeplimits');
% title(sprintf('%s',cruise),'FontWeight','Bold','Interpreter','none')
% annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
% if prtit; print(graphdevice,[way_images_flux cruise '_Tair' graphformat]); end;


%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('Day (UTC)');grid on;%grid minor
ylabel('Specific Humidty (g/kg)');
axis([stdt endt 0 30]);
set(gca(gcf),'XTick',stdt:xtickspace:endt);datetick('x','mm/dd','keepticks','keeplimits');
title(sprintf('%s',cruise),'FontWeight','Bold','Interpreter','none')
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
if prtit; print(graphdevice,[way_images_flux cruise '_SpecRH' graphformat]); end;


% % 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;

% xs=1; %step for x-axis tick
xs=xtickspace; %step for x-axis tick

figure;
set(gcf,'defaultaxesfonts',12)
set(gca,'fontname','time','fontsize',12);
subplot(6,1,1);plot(jdy,U,'b','LineWidth',2);
xlim([stdt endt]);ylabel('WS (ms^{-1})');
title(sprintf('%s.  5-min surface meteorological data',cruise),'FontWeight','Bold','Interpreter','none')
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
subplot(6,1,2);plot(jdy,dir,'b.');
axis([stdt endt 0 360]);ylabel('WD (deg)');
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
subplot(6,1,3);plot(jdy,ta,'b',jdy,tsnk,'r','LineWidth',2);
% axis([stdt endt 20 35]);legend('Ta', 'Ts','Orientation','horizontal','location','NorthWest');ylabel('Temp (degC)');
% xlim([stdt endt]);
axis([stdt endt 24 35])
legend('Ta', 'Ts','Orientation','horizontal','location','NorthWest');ylabel('Temp (degC)');
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
subplot(6,1,4);plot(jdy,psp,'b',jdy,rl,'r','LineWidth',2);
axis([stdt endt -5 1100]);%legend('SW' , 'LW','Orientation','horizontal','location','NorthWest');
ylabel('Radiation (Wm^{-2})');
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
subplot(6,1,5);plot(jdy,RH,'b','LineWidth',2);
axis([stdt endt 50 100]);ylabel('RH (%)');
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
subplot(6,1,6);plot(jdy,press,'b.');
axis([stdt endt median(press(~isnan(press)))-10 median(press(~isnan(press)))+10]);ylabel('BP (mbar)');
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
xlabel('Day (UTC)');
orient tall
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
if prtit; print(graphdevice,[way_images_flux cruise '_5min_Surface_MetData' graphformat]); end;


figure;
set(gcf,'defaultaxesfonts',12)
set(gca,'fontname','time','fontsize',12);
subplot(6,1,1);plot(jdy,taub,'b','LineWidth',2);
ylabel('Stress (Nm^{-2})');
title(sprintf('%s.  5-min BULK flux data',cruise),'FontWeight','Bold','Interpreter','none')
% axis([stdt endt -.05 0.5]);
axis([stdt endt 0 0.4]);
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
subplot(6,1,2);plot(jdy,hnet,'b','LineWidth',2);
ylabel('Hnet (Wm^-^2)');axis([stdt endt -400 1200])
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
subplot(6,1,3);plot(jdy,0.95*rs,'b','LineWidth',2);
ylabel('Rs (Wm^-^2)');axis([stdt endt -5 1250])
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
subplot(6,1,4);plot(jdy,rlnet,'b','LineWidth',2);
ylabel('Rnl (Wm^-^2)');axis([stdt endt -100 0])
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
subplot(6,1,5);plot(jdy,hsb,'b','LineWidth',2);
ylabel('Hsb (Wm^-^2)');axis([stdt endt -10 100])
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
subplot(6,1,6);plot(jdy,hlb,'b','LineWidth',2);
ylabel('Hlb (Wm^-^2)');axis([stdt endt -5 500])
set(gca(gcf),'XTick',stdt:xs:endt);datetick('x','mm/dd','keepticks','keeplimits');
grid on;%grid minor
xlabel('Day (UTC)');
orient tall
annotation(gcf,'textbox',[0.007154 0.01077 0.4498 0.02462],'String',{'NOAA/ESRL/PSD/Weather & Climate Physics'},'FontSize',6,'FitBoxToText','off','LineStyle','none'); 
if prtit; print(graphdevice,[way_images_flux cruise '_5min_BULK_FluxData' graphformat]); end;


% % % [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(press>1025);press(tt)=press(tt-1);
% % % tt=find(press<970);press(tt)=press(tt-1);
% % % [mean(press(~isnan(press))) std(press(~isnan(press))) max(press(~isnan(press))) min(press(~isnan(press)))]
% % % 
% % % [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')
