% epicz_read
% BWB for CAPRICORN 2016
% Reads flux_5min files & does plots and other stuff

%% initialize
clear all
close all
fclose all;
warning('off','MATLAB:MKDIR:DirectoryExists');

dv = clock;  % gets current date and time as date vector.
Vdate = datevec(datenum(dv(1),dv(2),dv(3)-1));
    
cruise = 'CAPRICORN_2016';  % Cruise name, e.g. 'STRATUS_2009', etc...
yyyy = '2016';              % year string
ship = 'Investigator';      % Research vessel
data_drive = 'D:\Data';     % define base path - DAS system

disp(['epicz_read: ' cruise]);

% paths relative to data_drive
path_proc_data = fullfile(data_drive,cruise,ship,'flux','Processed_Raw');
path_raw_images = fullfile(data_drive,cruise,ship,'flux','Raw_Images','epicz');
mkdir(path_raw_images);  % create directory if necessary

% matlab script path
path_prog = fullfile('c:\Users\Sergio\Documents\MATLAB\2016_CAPRICORN\flux');
cd(path_prog);                  % cd to matlab directory
restoredefaultpath              % reset matlab path
addpath(genpath(path_prog));	% add scripts to path
savepath; rehash;

% Adjustments & postcalibrations for PSD temperatures, all deg C
% tsea_adj = -0.5;  % sea snake high by 0.5 deg - 3/13/16 test
% ta_adj = 0;
% td1_adj = 0;
% tc1_adj = 0;
% td2_adj = 0;
% tc2_adj = 0;
% adj = [tsea_adj,ta_adj,td1_adj,tc1_adj,td2_adj,tc2_adj];

% flux configuration
% fsonic = 10;    % frequency (Hz)
% flicor = 10;    % frequency (Hz)

% Clear sky model configuration
k1 = 0.05;	% aerosol optical depth, band 1
k2 = 0.05;	% aerosol optical depth, band 2
oz = 0.2;	% column ozone
iv = 1.25;	% column water vapor (cm), if not calculated from obs.

% Track map limits
Lonmin = 130;
Lonmax = 160;
Latmin = -56;
Latmax = -40;

% sensor heights
zu_etl = 23-2;	% wind speed measurement height (m)
zt_etl = 21.5-2;	% air T measurement height (m)
zq_etl = 21.5-2;	% air q measurement height (m)
zu_ship = 25-2;	% scs wind speed measurement height (m)
zt_ship = 31.4-4;	% scs air T measurement height (m)
zq_ship = 31.4-4;	% scs air q measurement height (m)
dtsg = 7;   % depth of ship tsgm, tssm

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

%% concat 5-min flux files
filecat(fullfile(path_proc_data,[cruise,'_PSD_flux_5min_*.txt']),fullfile(path_proc_data,[cruise,'_PSD_flux_5min_all.txt']),'Overwrite');
filecat(fullfile(path_proc_data,[cruise,'_PSD_flux_30min_*.txt']),fullfile(path_proc_data,[cruise,'_PSD_flux_30min_all.txt']),'Overwrite');

%% read 5 min flux file with following format
%{
1  jd_5min        % Decimal julian day at beginning of line
2  wspd           % psd true wind speed at 23 m ht, m/s
3  wdir           % psd true wind direction, deg
4  Tsea           % psd seasnake T, C
5  tsgm           % ship theromsalinograph T, C
6  tssm           % ship theromsalinograph salinity, psu
7  Tvais          % psd air T at 18 m, C
8  qsea           % psd air specific humidity at sea surface, g/kg
9  qair           % psd air specific humidity at 18 m, g/kg
10  psp           % psd solar flux, w/m^2
11  pir           % psd IR flux, w/m^2
12  org           % psd optical raingage precip rate, mm/hr
13  sogm          % ship sog from gps, m/s
14  lrgm          % ship heading from gyrocompass, deg
15  rwspd         % rel wind speed, m/s
16  rwdir         % rel wind direction, deg
17  latm          % decimal latitude, deg
18  lonm          % decimal longtude, deg
19  ztsg          % depth of SST sensor used in heat flux calc, m
20  sogm_std      % standard deviation of ship speed, m/s
21  taub          % wind stress, COARE 3.5, N/m^2
22  hsb           % sensible heat flux, COARE 3.5, w/m^2
23  hlb           % latent heat flux, COARE 3.5, w/m^2
24  rf            % rain heat flux, w/m^2
25  imta          % IMET air temp, C (31.4 m)
26  imqa          % IMET air specific humidity, g/kg (31.4 m)
27  imum          % IMET true wind speed, m/s (25 m)
28  imdm          % IMET true wind direction, deg (25 m)
29  imsol         % IMET solar flux, w/m^2
30  imir          % IMET IR flux, w/m^2
31  impress       % BP, mb @ 20 m
32  imrh          % RH (%) (31.4 m)
33  U10N          % 10m neutral windspeed, COARE 3.5
34  Cdn_10        % 10m neutral drag coefficient, COARE 3.5
%}

x = load(fullfile(path_proc_data,[cruise,'_PSD_flux_5min_all.txt']));

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)
% remove data when snake is out of water during CAPRICORN
bbb = find(jdy<75 | (jdy>75.8&jdy<76.6) | (jdy>78.83&jdy<79.26) | ...
    (jdy>79.825&jdy<80.28) | (jdy>82.040&jdy<82.93) | (jdy>85.84&jdy<88.345) | ...
    (jdy>88.87&jdy<90.06));
tsnk(bbb) = NaN;

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_etl = x(:,8);    % psd sea surface specific humidity,from snake (g/kg)
qa = x(:,9);         % air specific humidity, PSD (g/kg) (z = 15 m)
psp = x(:,10);       % 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 sog (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
qa_im = x(:,26);     % IMET air specific humidity, g/kg
s_ship = x(:,27);    % IMET true wind speed, m/s
dir_ship = x(:,28);  % IMET true wind direction, deg
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, %

% choose selected values for later use
sog = x(:,13);       % same as ushp
cog = x(:,14);       % use heading
% org = zeros(length(org1),1);
org = org1;
rs = psp_ship;
% tsea = tsnk;
tsea = tsg;
qse = qsea(tsea);

% set date range as desired
stdt = floor(jdy(1));
endt = ceil(jdy(end)); 

%% plot track map
figure;
m_proj('mercator','longitudes',[Lonmin Lonmax],'latitudes',[Latmin Latmax]);
m_grid;
m_coast('linewidth',2,'color',0.5+[0 0 0]);
m_coast('patch',[.7 .7 .7],'edgecolor','none');
m_line(lon,lat,'color','r','linestyle','none','marker','.','markersize',6);
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');
if prtit;
    ppath = fullfile(path_raw_images,[cruise,'_ship_track_all',graphformat]);
    print(graphdevice,ppath);
end;

%% run bulk model

% run bulk model with 5 min data
% using vectorized COARE 3.5 with output adjusted to match historical cor30b
% y=[hsb hlb tau zo zot zoq L usr tsr qsr dter dqer tkt RF wbar Cd Ch Ce Cdn_10 Chn_10 Cen_10 ug ];
%     1   2   3   4  5   6  7  8   9  10   11   12  13  14  15  16 17 18    19      20    21  22
% coare35vn adjusted to return U10N as separate output array, 3/2016 BWB

[yy,U10_et] = coare35vn(U,zu_etl,ta,zt_etl,RH,zq_etl,press,tsea,rs,rl,lat,600,org,NaN,NaN);

zou = yy(:,4);
zoq = yy(:,6);
ust = yy(:,8);
cdn10 = yy(:,19);
cen10 = yy(:,21);
hlb = yy(:,2);
taub = yy(:,3);
hsb = yy(:,1);
dter = yy(:,11);
dqer = yy(:,12)*1e3; % in g/kg
L = yy(:,7);
qsr = yy(:,10);
tsr = yy(:,9);
wbar = yy(:,15);

U10_im = s_ship+ust/.4.*(log(10./zu_ship)-Psit_26(10./L)+Psit_26(zu_ship./L));

%% constants
tdk = 273.16;
grav = 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 = isfinite(tsea);
Le = (2.501-.00237*median(tsea(ii)))*1e6;

%% radiation models & heat budget
% Downwelling IR
rlclr = (0.52+0.13/60*abs(lat)+(0.082-0.03/60.*abs(lat)).*sqrt(qa)).*(5.67e-8*(ta+tdk).^4);

% Downwelling Solar
latd = nanmean1(lat);
lond = nanmean1(lon);
watvap = iv;
jdx = jdy + double(2.5/1440); % 5-min jd bin centers
tutc = (jdx - floor(jdx)).*24;    % decimal hour of the day, 0-24
clear sw sz saz dirs sky;
[Rscl,sz,saz,dirs,sky] = SolarRadiance(latd,lond,jdx,tutc,watvap,nanmean1(press),k1,k2,oz);

rsclr = Rscl;
rlnet = 0.97*(rl-5.67e-8*(tsea+tdk).^4);
rlnetc = 0.97*(rlclr-5.67e-8*(tsea+tdk).^4);
hnet = 0.95*rs+rlnet-hsb-hlb-hrain;

%% further processing
cu = .00087; bu = .00018; cx = .00085; bx = .00015;
dts = tsea-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(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;

hl_webb = rhoair.*wbar.*qa.*Le/1000;

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;
%ii = find(dq10n<-5);dq10n(ii) = qse(ii)-qp(ii);

%% plots
xtickspace = 2;

% Temperatures
figure; plot(jdy,tsea,'b.',jdy,ta,'g'); xlim([stdt endt]);
set(gca,'fontsize',12);
xlabel('Day (UTC)'); ylabel('Temperature (degC)');grid on; %grid minor;
title(sprintf('%s.  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;
    ppath = fullfile(path_raw_images,[cruise,'_SST_Tair_comparison',graphformat]);
    print(graphdevice,ppath);
end;

% Sp Hum
figure;plot(jdy,qa,'b');
set(gca,'fontsize',12); xlim([stdt,endt]);
xlabel('Day (UTC)');grid on; grid minor;
ylabel('Specific Humidty (g/kg)');
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;
    ppath = fullfile(path_raw_images,[cruise,'_SpecRH',graphformat]);
    print(graphdevice,ppath);
end;

% Sonic Winds
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('Wspd (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} (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} (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;
    ppath = fullfile(path_raw_images,[cruise,'_PSD_Sonic_speeds',graphformat]);
    print(graphdevice,ppath);
end;

% True Wspd/Wdir
figure;
subplot(2,1,1);set(gca,'fontsize',12);
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');
xlim([stdt endt]);
ylabel('Speed (ms^{-1})');

subplot(2,1,2);set(gca,'fontsize',12);
plot(jdy,dir,'.'); axis([stdt endt 0 360]);
grid on; %grid minor;
set(gca(gcf),'XTick',stdt:xtickspace:endt);datetick('x','mm/dd','keepticks','keeplimits');
set(gca(gcf),'YTick',0:90:360)
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;
    ppath = fullfile(path_raw_images,[cruise,'_PSD_Sonic_SpdDir',graphformat]);
    print(graphdevice,ppath);
end;

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

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

jdq = jdy;%

eval('epicz_daily_CAPRICORN_2016');

%% epics diurnal
jdq = jdy;

eval('epicz_diurnal_CAPRICORN_2016');

%% lat/lon plots
stlon = min(lon)-1; enlon = max(lon)+1;
stlat = min(lat)-1; enlat = max(lat)+1;
ii = find(jdy>=stdt & jdy<endt);

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

subplot(3,1,2); set(gca,'fontsize',10);
plot(lon(ii),tsnk(ii),'k.',lon(ii),ta(ii),'r.');
xlim([stlon enlon]);ylabel('Temp (C)');
grid on;%grid minor
legend('SST','Ta','location','Best')

subplot(3,1,3); set(gca,'fontsize',10);
plot(lon(ii),hsb(ii),'k.',lon(ii),hlb(ii),'r.');
axis([stlon enlon -20 400]);xlabel('Long (deg)');ylabel('Heat (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;
    ppath = fullfile(path_raw_images,[cruise,'_U_TaSST_HsHl_lon',graphformat]);
    print(graphdevice,ppath);
end;

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

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

subplot(3,1,3); set(gca,'fontsize',10);
plot(lat(ii),hsb(ii),'k.',lat(ii),hlb(ii),'r.');
xlim([stlat,enlat]); xlabel('Lat (deg)'); ylabel('Heat (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;
    ppath = fullfile(path_raw_images,[cruise,'_U_TaSST_HsHl_lat',graphformat]);
    print(graphdevice,ppath);
end;

%% grid plots

% %[xi,yi] = meshgrid(-71:.25:-66,41.5:.25:44.5); %EPICS?
% %[xi,yi] = meshgrid(-60:.25:-35,-55:.25:-48); %GASEXIII
% [xi,yi] = meshgrid(stlon:0.25:enlon, stlat:0.25:enlat);
% 
% ia = find(isfinite(lon) & isfinite(U));
% ZI = griddata(lon(ia),lat(ia),U(ia),xi,yi);
% figure;[c,h] = contour(xi,yi,ZI); clabel(c,h), colorbar;  grid;
% ylabel('Latitude'); xlabel('Longitude'); title('Wind Speed (m/s)');
% 
% ia = find(isfinite(lon) & isfinite(ta-tsea));
% ZI = griddata(lon(ia),lat(ia),ta(ia)-tsea(ia),xi,yi);
% figure;[c,h] = contour(xi,yi,ZI); clabel(c,h), colorbar; grid;
% ylabel('Latitude'); xlabel('Longitude'); title('Air-Sea T difference (C)');
% 
% ia = find(isfinite(lon) & isfinite(tsea));
% ZI = griddata(lon(ia),lat(ia),tsea(ia),xi,yi);
% figure;[c,h] = contour(xi,yi,ZI); clabel(c,h), colorbar; grid;
% ylabel('Latitude'); xlabel('Longitude'); title('Sea T C)');
% 
% ia = find(isfinite(lon) & isfinite(stabeff));
% ZI = griddata(lon(ia),lat(ia),stabeff(ia),xi,yi);
% figure;[c,h] = contour(xi,yi,ZI); clabel(c,h), colorbar; grid;
% ylabel('Latitude'); xlabel('Longitude'); title('Stability Effect on Fluxes at 10-m');
% 
% %% ASSMAN intercomparison
% if ass
%    % 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;

%% big met plots
xs = xtickspace; %step for x-axis tick

h1 = figure;
set(gca,'fontsize',10);
subplot(6,1,1);plot(jdy,U,'b','LineWidth',2);
xlim([stdt endt]); ylabel('U');
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('Dir');
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,tsea,'r','LineWidth',2);
xlim([stdt endt]); ylabel('T');
legend('Ta', 'Ts','Orientation','horizontal','location','best');
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);
xlim([stdt endt]); ylim([-5,max(psp)+20]);
ylabel('Rl/Rs');
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.');
xlim([stdt endt]); ylabel('Pmb');
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'); 
ppath = fullfile(path_raw_images,[cruise,'_5min_Surface_MetData.fig']);
saveas(h1,ppath);


h2 = figure;
set(gca,'fontsize',10);
subplot(6,1,1); plot(jdy,taub,'b','LineWidth',2);
ylabel('taub'); xlim([stdt,endt]);
title(sprintf('%s.  5-min BULK flux 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,hnet,'b','LineWidth',2);
ylabel('Hnet'); xlim([stdt,endt]); ylim([-500,800]);
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'); xlim([stdt,endt]); ylim([-10,800]);
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') ; xlim([stdt,endt]); ylim([-150,50]);
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'); xlim([stdt,endt]);
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'); xlim([stdt,endt]);
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'); 
ppath = fullfile(path_raw_images,[cruise,'_5min_Bulk_FluxData.fig']);
saveas(h2,ppath);


% % % [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')
