% HiWinGS_2013_da_red1.m - process output da routines

% THIS VERSION USING REFINED BULK MET DATA FROM THE 10-MIN MEAN MET FILE
%       HiWinGS_2013_10minMeans_r2.mat
%

%% initialize
close all;
clear;

%% set run parameters
prtit = true;       % set true to save flux result files to disk
savefigs = true;    % set true to save plots as png files

% datasetv = 'v1';  % read da files from standard Hare / Edson motion correction
% datasetv = 'v2';  % read v2 version da fles with modified tilt correction
datasetv = 'v3';    % v2 with decorr of Tsonic and q-licor

% set flag for decorrelation
decorr = true;     % true for da version with decorr
if decorr
    dcname ='decorr';
else
    dcname = 'std';
end;

% set limits for filter variables to be used in all routines
rdir_lo     = -60;  % rel wind direction, deg from bow
rdir_hi     = 60;
rain_lim    = 5;    % rain rate, mm/hr
sp2_lim     = 0.8;  % vplat std dev  m/s
sigu_lim    = 0.6;  % ship speed std dev  m/s
sigh_lim    = 5;    % ship heading std dev, deg
badSon_lim  = 50;   % # bad sonic points per 10min interval (PSD only)
missing_lim = 50;   % # missing sonic points per 10min interval(PSD only)
ushp_lim    = 3.5;  % ship speed, m/s
tilt_lim    = 10;   % flow tilt angle, degrees
badlic_lim  = 50;   % # bad licor points per 10min interval
missLic_lim = 50;   % # missing licor points per 10min interval
missMot_lim = 50;   % # missing licor points per 10min interval

% choose ID method, 'a' or 'b' only
IDflag = 'a';

% for map plot...
Lonmin = 290;   % longitudes are 0-360, going E
Lonmax = 320;
Latmin = 40;    % deg N is positive
Latmax = 72;
LatLonLim = [Lonmin Lonmax Latmin Latmax];

%% read files with 10-min results

% path to 'Processed' folder da/cr files
da_path = '/Users/byronb/Raw_Data_Files.noindex/HiWinGS_2013/Knorr/flux/Processed_v3';       % for Mac file system
% path to folder for saving plots
fig_path = '/Users/byronb/Raw_Data_Files.noindex/HiWinGS_2013/Knorr/flux/Processed_images_v3/da_red';
mkdir(fig_path);  % create directory if it doesn't exist

x = load(fullfile(da_path,['HiWinGS_2013_da_',datasetv,'_',dcname,'.txt']));       % 219 column array
yxx = load(fullfile(da_path,['HiWinGS_2013_cr_',datasetv,'_',dcname,'.txt']));     % 17 column array

% load HiWinGS refined mean met file, 5135 x 22 array with column format:
% Yday Lat Lon SOG COG Heading Cspd Cdir Wspd Wdir WspdR WdirR rdir Pair RH Tair Qair Tsea Tsnake Solar IR  Pcp
%  1    2   3   4   5     6     7    8    9    10   11    12    13   14  15  16   17   18    19    20   21  22
load(fullfile(da_path,'HiWinGS_MM10_r2.mat'));

% define variables from loaded arrays
% see readme on da file or check run_motcorr_flux script for column descriptions
jdy = x(:,1);           % decimal julian day, 10-min timebase
stdt = floor(jdy(1)); endt = floor(jdy(end))+1;

%% define variables
% ship imet / scs data
U_scs   =   x(:,2);     % true wind speed, m/s; imet sonic anemometer
dir_scs =   x(:,34);    % true wind direction (from), deg (clockwise rel north)
tsg     =   x(:,6);     % tsg water temperature, C (5 m depth)
ta_im   =   x(:,3);     % imet air temperature, C
rs_im   =   x(:,7);     % downward solar flux, W/m^2 (imet units)
rl_im   =   x(:,8);     % downward ir flux, W/m^2 (imet units)
rain_im =   x(:,10);    % precip, mm/hr
zu_im   =   ones(size(U_scs))*15;    % measurement heights for ship sensors, m
zt_im   =   zu_im;
zq_im   =   zu_im;
rh_im   =   x(:,48);    % SCS RH
P_im    =   x(:,5);     % SCS BP mb
qs_tsg  =   qsea_p([tsg,P_im]);             % imet sea surface specific humidity, g/kg
qa_im   =   qair_p([ta_im,rh_im,P_im]);     % imet air specific humidity, g/kg
zP_ship =   zu_im;      % SCS pressure sensor height

% PSD data
U       =   x(:,16);    % true wind from PSD sonic - relative to earth
dir     =   x(:,41);    % true wind direction, PSD sonic
urel    =   x(:,38);    % relative wind speed, m/s
reldir  =   x(:,42);    % relative wind direction, +/- 180 deg rel ship's bow (PSD sonic)
tsnk    =   x(:,20);    % sea snake SST
ta      =   x(:,17);    % vaisala air temperature
rh      =   x(:,18);    % vaisala rh
rs      =   (x(:,21)+x(:,216))/2;    % downward solar flux, W/m^2 - avg of psp1 & psp2
sig_rl  =   x(:,95);    % stdev downward solar flux, W/m^2 (PSD units)
rl      =   (x(:,22)+x(:,217))/2;    % downward IR flux, W/m^2 - avg of pir1 & pir2
sig_rs  =   x(:,96);    % stdev downward IR flux, W/m^2 (PSD units)
org     =   x(:,24);    % rainrate, mm/hr (psd rain gauge, uncorrected)
zu_etl  =   ones(size(U_scs))*15.91;    % measurement heights
zt_etl  =   ones(size(U_scs))*14;
zq_etl  =   ones(size(U_scs))*14;
press   =   x(:,19);    % psd vaisala pressure (mbar)
qse     =   qsea_p([tsnk,press]);   % sea surface specific humidity, g/kg (from snake)
qa      =   qair_p([ta,rh,press]);  % air specific humidity, g/kg
zP      =   10;         % psd pressure sensor height, m

% UH wind data
U_uh       =  x(:,211);
dir_uh     =  x(:,212);
relspd_uh  =  x(:,207);
reldir_uh  =  x(:,208);

% 10-min results from initial C35 run
Hsb     =   x(:,76);    % bulk sensible heat flux, W/m^2
Hlb     =   x(:,77);    % bulk latent heat flux, W/m^2
Taub    =   x(:,75);    % bulk wind stress along mean wind, N/m^2
ustr0   =   x(:,74);    % u-star
tstr0   =   x(:,79);    % t-star
qstr0   =   x(:,80);    % q-star
Hrain   =   x(:,91);    % rain heat flux, W/m^2 as per Gosnell et. al.

% Navigation & data quality indicators
Lat     =   x(:,151);   % latitude, deg, SCS
Lon     =   x(:,150);   % longitude, deg, SCS, 0-360 deg going E
Lon(Lon<0) = Lon(Lon<0) + 360;  % reformat negative W longitudes
cog     =   x(:,105);   % cog, deg, PSD gps
cog(cog<0) = cog(cog<0)+360;
sog     =   x(:,109);   % sog, m/s, PSD gps
tiltx   =   x(:,44)*180/pi; % wind flow tilt, deg, PSD sonic
tilt_uh =   x(:,214)*180/pi; % tilt, UH sonic
head    =   x(:,112);   % ship heading, deg clockwise rel north, ship's gyro
sig_h   =   x(:,113);   % std dev of heading, >5 indicates ship maneuver
sog_im  =   x(:,97);    % sog, m/s (SCS GPS)
sig_uim =   x(:,98);    % std dev of ship speed, m/s (>.2 indicates maneuver)
cog_im  =   x(:,101);   % scs cog
cog_im(cog_im<0) = cog_im(cog_im<0)+360;
sig_him =   x(:,102);   % std dev of ship cog (SCS)
ushp    =   sog_im;     % doppler speed log (odec), ship - NOT WORKING, USE SOG
sig_u   =   sig_uim;    % std dev of speed log - NOT WORKING USE STD SOG
missing =   x(:,158);   % missing lines of sonic data in each 10-min interval
missing(isnan(missing)) = 0;
badSon  =   x(:,159);   % number of NaNs in sonic data in each 10-min interval
badSon(isnan(badSon)) = 0;
Ts_noise =  x(:,218);   % psd Tsonic noise variance
Ts_noise_u = x(:,219);  % uh Tsonic noise variance
missing_lic = x(:,220); % missing lines of licor data in each 10-min interval
badLic  = x(:,221);     % number of NaNs in licor H2O data in each 10-min interval
missing_mot = x(:,222); % missing lines of motion data in each 10-min interval

% variances and covariances, turbulence data
wuj      =   x(:,117);   % wu cov
wu_uh    =   x(:,185);
wvj      =   x(:,119);   % wv cov
wv_uh    =   x(:,186);
wt_son   =   x(:,121);   % wTs cov
wt_son_uh =  x(:,187);
uuj      =   x(:,125);   % u var
uu_uh    =   x(:,188);
vvj      =   x(:,126);   % v var
vv_uh    =   x(:,189);
wwj      =   x(:,127);   % w var
ww_uh    =   x(:,190);
tt       =   x(:,128);   % sonic T var
tt_uh    =   x(:,191);

sp1      =   sqrt(x(:,129));   % ship motion stdev, u direction (m/s)
sp1_uh   =   sqrt(x(:,192));
sp2      =   sqrt(x(:,130));   % ship motion stdev, v direction (m/s)
sp2_uh   =   sqrt(x(:,193));
sp3      =   sqrt(x(:,131));   % ship motion stdev, w direction (m/s)
sp3_uh   =   sqrt(x(:,194));

% select inertial dissipation results
switch IDflag
    case 'a'
        cu       =   x(:,134);   % cu^2 U method a
        cu_uh    =   x(:,195);
        cw       =   x(:,136);   % cw^2 W method a
        cw_uh    =   x(:,197);
        ct       =   x(:,138);   % ct^2 temperature method a
        ct_uh    =   x(:,199);
        cql      =   x(:,140);   % cq^2 H2O Licor method a
        usi      =   x(:,142);   % inertial dissipation scaling paramters, all method a
        usi_uh   =   x(:,201);
        tsi      =   x(:,144);
        tsi_uh   =   x(:,203);
        qsi      =   x(:,146);
    case 'b'
        cu       =   x(:,135);   % cu^2 U method b
        cu_uh    =   x(:,196);
        cw       =   x(:,137);   % cw^2 W method b
        cw_uh    =   x(:,198);
        ct       =   x(:,139);   % ct^2 temperature method b
        ct_uh    =   x(:,200);
        cql      =   x(:,141);   % cq^2 H2O Licor method b
        usi      =   x(:,143);   % inertial dissipation scaling paramters, all method b
        usi_uh   =   x(:,202);
        tsi      =   x(:,145);
        tsi_uh   =   x(:,204);
        qsi      =   x(:,147);
end;

%% make adjustments, calibrations, etc.
% apply fixes to data and select best quality measurements
% for use in bulk flux calculations.  Edit this script to
% make adjustments for 'best' bulk model input dataset
HiWinGS_2013_fixit_da1

%% correct U for ocean surface currents
% otherwise current data must be computed from speed log or ADCP...
% assume current is specified in oceanographic convention - direction 'to'
%%% set current to zero for now - HiWinGS u_best already current corrected %%%

rdcon = pi/180;
% CNa = -cspd_a.*cos(cdir_a./180.*pi);  % current 'toward' North (meteorological convention)
% CEa = -cspd_a.*sin(cdir_a./180.*pi);  % current 'toward' East
% CN = intavg(jd_scs_a,CNa,jdy);        % average to jdy 10-min timebase
% CE = intavg(jd_scs_a,CEa,jdy);
CN = zeros(length(jdy),1);              % set current to zero for now
CE = zeros(length(jdy),1);

% compute U and dir relative to water
UN = u_best.*cos(dir_best./180.*pi); % U North
UE = u_best.*sin(dir_best./180.*pi); % U East
Ul = sqrt((UN+CN).^2+(UE+CE).^2);                   % Wind speed (m/s) relative to water
dirl  = mod(atan2((UE+CE),(UN+CN))/rdcon+360,360);  % Wind direction (from) relative to water

Uearth = u_best; direarth = dir_best;     % save copies U & dir in earth coords
u_best = Ul; dir_best = dirl;             % use current relative winds for U/dir

% append adjusted wind speeds and best wind speeds to x array
x(:,220) = Uearth;      % U in earth coords appended to da array
x(:,221) = direarth;    % dir in earth coords appended to da array
x(:,222) = u_best;      % ocean surface relative winds
x(:,223) = dir_best;

%% constants and other variables for the bulk flux calculations
tdk = 273.16;
grav = 9.72;
Rgas = 287.1;
cpa = 1004.67;
be = 0.026;
cpw = 4000;
rhow = 1022;
visw = 1e-6;
tcw = 0.6;
ztsg = ones(length(jdy),1)*5;	% 5m depth of tsg

rhoair = P_best.*100./(Rgas*(ta_best+tdk).*(1+0.61*qa_best/1000));
Le = (2.501-.00237*nanmedian1(ts_best))*1e6; % cruise mean Le

% average cu and cw for smoother structure function
cs = (cu+.75*cw)/2;
ii = find(isnan(cu));
cs(ii) = .75*cw(ii);
kk = find(isnan(cw));
cs(kk) = cu(kk);

cs_uh = (cu_uh+.75*cw_uh)/2;
ii = find(isnan(cu_uh));
cs_uh(ii) = .75*cw_uh(ii);
kk = find(isnan(cw_uh));
cs_uh(kk) = cu_uh(kk);

% plot for map of cruise track: m_map folder must be in MATLAB path!
% using Albers equal-area projection for high lat map...
map_HiWinGS_2013(Lon,Lat,'HiWinGS_2013',LatLonLim);

% compute clear sky downwelling longwave flux
rlcler = 5.67e-8*(ta_best+273.16).^4.*(.52+.085*sqrt(qa_best));   % model clear sky
rlclr = (.52+.13/60*abs(Lat)+(.082-.03/60*abs(Lat)).*sqrt(qa_best)).*(5.67e-8*(ta_best+tdk).^4); % model as f(lat)
rlcler = rlclr; % use f(lat) model

% shortwave clear sky radiative flux
js = stdt; je = endt;
iv = ones(length(jdy),1)*NaN;
min_set = 5;
HiWinGS_2013_Rsclear;     % runs SolarRadiance model

%% run COARE for 10-min bulk model estimates

zi = 600;     % inversion ht default
y = coare35vn(u_best,zu,ta_best,zt,rh_best,zq,P_best,ts_best,rs_best,rl_best,Lat,zi,rain_best,NaN,NaN);
% output of C35 set to match old C30 format + U10N
% y = [hsb hlb tau zo zot zoq L usr tsr qsr dter dqer tkt RF wbar Cd Ch Ce Cdn_10/1000 Chn_10/1000 Cen_10/1000 ug U10N];
%       1   2   3   4  5   6  7  8   9  10   11   12  13  14  15  16 17 18    19           20          21      22  23
hsb    =    y(:,1);       % sensible heat flux W/m2
hlb    =    y(:,2);       % latent heat flux W/m2
tau    =    y(:,3);       % stress N/m2
zo     =    y(:,4);       % vel roughness length
zot    =    y(:,5);       % temp "
zoq    =    y(:,6);       % hum  "
L      =    y(:,7);       % Obukhov Length, m
usr    =    y(:,8);       % ustar
tsr    =    y(:,9);       % tstar
qsr    =    y(:,10);      % qstar (g/kg)
dter   =    y(:,11);      % cool skin delta T
dqer   =    y(:,12);      %  "   "      "   q
tkt    =    y(:,13);      % thickness of cool skin
RF     =    y(:,14);      % rain heat flux
Wbar_webb = y(:,15);      % webb mean w
Cd     =    y(:,16);      % drag @ zu
Ch     =    y(:,17);      %
Ce     =    y(:,18);      % Dalton
Cdn_10 =    y(:,19);      % neutral drag @ 10m [includes gustiness]
Chn_10 =    y(:,20);      %
Cen_10 =    y(:,21);      %
ug     =    y(:,22);      % gustiness
u10n   =    y(:,23);      % 10m neutral windspeed, m/s

% dt is expanded output array
dt = NaN(length(jdy),35);
dt(:,1)    = jdy;    % julian day
dt(:,2:10) =  [u_best,sog_im,ts_best,ta_best,qs_best,qa_best,rs_best,rl_best,rain_best];
dt(:,11:32)  = y(:,1:22);    % standard COARE output
dt(:,33:34) = 0;     % no warm layer results; set COLS 33:35 to zero
dt(:,35) = ts_best;

%% further processing
% select results from dt
pst = psitg(zq./dt(:,17));
[np,a] = size(dt);
dtw = dt(:,33);
zw = dt(:,34);
Hsbb = dt(:,11);
Hlbb = dt(:,12);
HRF = dt(:,24);
Taubb = dt(:,13);

% compute hsc - u* and q* have been updated by new coare run
wt = wt_son + 0.51*(ta_best+tdk).*usr.*qsr/1000;  % correct sonic wt for humidity contribution
wt_uh = wt_son_uh + 0.51*(ta_best+tdk).*usr.*qsr/1000;
hsc = wt*cpa.*rhoair;
hsc_uh = wt_uh*cpa.*rhoair;
hsc_son = wt_son*cpa.*rhoair;
hsc_son_uh = wt_son_uh*cpa.*rhoair;

% compute 'best' 10-min hsc & tau as avg of filtered psd and uh results
% find periods of bad data...
jjj = find((reldir_best>rdir_hi | reldir_best<rdir_lo) | sig_h>sigh_lim | ...
            sp2_best>sp2_lim | sig_u>sigu_lim | badSon>badSon_lim);
hsc_filt = hsc; hsc_filt(jjj) = NaN;
tauc_filt = -(rhoair.*wuj); tauc_filt(jjj) = NaN;
taux_filt = -(rhoair.*wvj); taux_filt(jjj) = NaN;

jjj = find((reldir_best>rdir_hi|reldir_best<rdir_lo) | sig_h>sigh_lim | sig_u>sigu_lim);
hsc_uh_filt = hsc_uh; hsc_uh_filt(jjj) = NaN;
tauc_uh_filt = -(rhoair.*wu_uh); tauc_uh_filt(jjj) = NaN;
taux_uh_filt = -(rhoair.*wv_uh); taux_uh_filt(jjj) = NaN;

hsc_best = nanmean1([hsc_filt';hsc_uh_filt']); hsc_best = hsc_best';
tauc_best = nanmean1([tauc_filt';tauc_uh_filt']); tauc_best = tauc_best';
taux_best = nanmean1([taux_filt';taux_uh_filt']); taux_best = taux_best';

% compute hl_webb correction from C35 Wbar_webb
hl_webb = Wbar_webb.*rhoair*Le/1000.*qa;

% select best covariance latent heat
hlc_lic = wq_lic*Le.*rhoair/1000;
hlc_bl = hlc_lic; %  Licor only since 2006.

% file name & plot title text string used in subsequent calls
jdollars = 'HiWinGS-2013';

%% inertial dissipation processing
HiWinGS_2013_da_id2;
hlib_lic(hlib_lic>1000) = NaN;	% remove bad id hl values

%% ocean warm layer profiles - does nothing for HiWinGS since warm layer not computed
dtz = zeros(np,6);
zz=0:5;
for i=1:6
   ii = find(zw>zz(i));
   dtz(ii,i) = dtw(ii).*(1-zz(i)./zw(ii));  % warming from mixed layer to z
end;
dt5 = dtz(:,6);
tmxd = tsg-dtz(:,6);  % mixed layer temp: tsg with with 5-m warm layer removed

%% 10-minute net heat flux
rnl = 0.97*(5.67e-8*(ts_best+273.15).^4 - rl_best);                 % net longwave heat flux
hnet = 0.955*rs_best - rnl - Hsbb - Hlbb - HRF;                     % net heat from coare
hnet_c = 0.955*rs_best - rnl - hsc_best - hlc_bl - hl_webb - HRF;   % net heat using cov fluxes
hnet_id = 0.955*rs_best - rnl - hsib - hlib_lic - hl_webb - HRF;    % net heat using id fluxes

%% compute daily avgs & save results

HiWinGS_2013_daily; % cloud plots

%% save ascii file of 10-min flux results if desired
if prtit
    prt_HiWinGS_2013_flux_10;
end;

%% now compute 1-hr averages & rerun bulk model
HiWinGS_2013_hr_ave;

% recomputing bulk values with hourly avg input variables using COARE 3.5
% output: yhsb, yhlb, ytaub, yLb, yusb, ytsb, yqsb, ydter, yhrain, etc...
% also plots updated bulk values with best licor turb hl
HiWinGS_2013_hr_bulk;

%% median hourly hl & stress results vs bulk and vs U
med_hl_plt_HiWinGS_2013;

med_tau_plt_HiWinGS_2013; % stress plots

% Hourly net heat flux
% net longwave relative to sea from warm-layer-corrected Ts (up flux is negative)
rnl_hr = 0.97*(5.67e-8*(yts + 273.15).^4 - y1(:,8));
% hnet from COARE 3.5 using ship bulk met. y1(:,8) is ship PIR.
hnetb_hr = 0.955*y1(:,8) - rnl_hr - yhsb - yhlb - yhrain;
% net heat using hsib and  hlib only due to problems with hlc
hnet_hr = 0.955*y1(:,8) - rnl_hr - yhsib - (yhlib_lic + hl_webar) - yhrain;

if prtit
    prt_HiWinGS_2013_flux_hr;     % save hourly flux to file if desired
end;

% ii = find( isfinite(hsbb) & isfinite(hlbb));
% disp('    hsbb      hlbb       rs       Rscl       rl       rlclr     rnl        hnet')
% nanmean1([ hsbb(ii) hlbb(ii) rs(ii) Rscl(ii) rl(ii) rlclr(ii) rnl(ii) hnet(ii)])
% ii = find( isfinite(yhlc_bl) & isfinite(yhlb) & isfinite(hlc_licbar));
% disp('    yhlb     yhlc_bl  hlc_licbar hl_webar')
% nanmean1([ yhlb(ii)' yhlc_bl(ii)' hlc_licbar(ii)' hl_webar(ii)])
% ii = find( isfinite(yhlc_bl) & isfinite(yhlb) & isfinite(hlc_licbar));
% disp('    yhlb     yhlc_bl  hlc_licbar hl_webar   -  medians')
% nanmedian1([ yhlb(ii)' yhlc_bl(ii)' hlc_licbar(ii)' hl_webar(ii)])

% webb correction from hourly bulk hs and hl
rhobar = yP.*100./(Rgas*(yta+tdk).*(1+0.61*yqa/1000));     % hr mean air density
wbar_t = (1+1.61*yqa/1000).*yhsb./rhobar./cpa./(yta+tdk);  % webb correction
wbar_w = 1.61*yhlb/Le;

%% define selection criteria for plots of 10-min data - showing PSD results only
jjj = find((reldir_best<rdir_hi & reldir_best>rdir_lo) & sig_h<sigh_lim & sig_u<sigu_lim & ...
            sp2_best<sp2_lim & (sqrt(ww_best)./ugw < 2+.0015*u_best.^2) & (sqrt(vv_best)./ugu < 2) & missing<missing_lim);

% various plots of 10-min and hourly results
z_line = jdy*0;

figure;
plot(jdy,missing,'b.',jdy,missing_mot,'r.',jdy,missing_lic,'g.',jdy,z_line,'k-')
legend('Sonic','Motion','Licor','Location','Best');
grid; xlim([stdt endt]);
xlabel('DOY (2014)'); ylabel('missing data points per 10min');
title([jdollars,' Check missing PSD 10Hz data']);

figure;
minx = min(taubb(jjj))-0.02; maxx = max(taubb(jjj))+0.02;
plot(jdy,z_line,'k-',jdy(jjj),-(rhoair(jjj).*wuj(jjj)),'bo',jdy(jjj),(rhoair(jjj).*usib(jjj).^2),'m+',jdy(jjj),taubb(jjj),'g.')
legend('','Cov','ID','Bulk','Location','Best');
xlabel('DOY'); ylabel('Stremwise stress (N/m2)');
grid; xlim([stdt endt]); ylim([minx maxx]);
title([jdollars,' Stress: Filtered 10-min Cov, ID and COARE']);

figure;
minx = min(hsib(jjj))-10; maxx = max(hsib(jjj))+10;
plot(jdy(jjj),hsc(jjj),'bo',jdy(jjj),hsib(jjj),'m+',jdy,hsbb,'g.',jdy,z_line,'k-')
legend('Cov','ID','Bulk','Location','Best');
grid; xlim([stdt endt]); ylim([minx maxx])
xlabel('DOY'); ylabel('Sensible heat flux (W/m2)');
title([jdollars,' Sensible Heat: Filtered 10-min Cov, ID and COARE']);

figure;
minx = min(hlbb)-30; maxx = max(hlbb)+30;
plot(jdy(jjj),hlc_bl(jjj),'bo',jdy(jjj),hlib_lic(jjj),'m+',jdy,hlbb,'g.',jdy,z_line,'k-')
legend('Cov','ID','Bulk','Location','Best');
grid; xlim([stdt endt]); ylim([minx maxx]);
xlabel('DOY'); ylabel('Latent heat flux (W/m2)');
title([jdollars,' Latent Heat: Filtered 10-min Cov, ID and COARE']);

figure;
minx = min(hnet)-30; maxx = max(hnet)+30;
h1 = plot(jdy,z_line,'k-',jdy(jjj),hnet_id(jjj),'mo',jdy,hnet,'g.');
h1 = h1(2:end); % select traces for legend, omit zero line
xlabel('DOY'); ylabel('Net heat flux (W/m^2)'); xlim([stdt endt]);
ylim([minx maxx]); legend(h1, 'ID Hl & Hs','Bulk Hl & Hs','Location','Best'); grid;
title([jdollars,' Filtered 10-min Net Heat Flux: Hs & Hl from ID & COARE']);

% hourly plots
z_line2 = jda*0;
figure;
minx = min(-yhlb)-30; maxx = max((0.95*y1(:,7)))+30;
h2 = plot(jda,z_line2,'k',jda,-(rnl_hr),'r',jda,(0.95*y1(:,7)),'b',jda,-yhsb,'g',jda,-yhlb,'m', jda,-yhrain,'c');
h2 = h2(2:end); % select traces for legend, omit zero line
xlabel('DOY'); ylabel('Heat Flux Component (W/m^2)'); xlim([stdt endt]);
ylim([minx maxx]); legend(h2, 'net IR','Solar','Bulk Hs','Bulk Hl','Bulk RF','Location','Best'); grid;
title([jdollars,' Net Hourly Heat Flux Components from COARE']);

figure;
minx = min(yhsib)-10; maxx = max((yhlib_lic + hl_webar))+30;
h3 = plot([minx maxx],[minx maxx],'k',[minx maxx],[0 0],'k',[0 0],[minx maxx],'k',yhsb,yhsib,'or',yhlb,(yhlib_lic + hl_webar),'db');
h3 = h3(4:end); % select traces for legend, omit 1:1 line
xlabel('Bulk Heat Flux (W/m^2)'); ylabel('Turb Heat Flux (W/m^2)');
axis([minx maxx minx maxx]); legend(h3,'Sensible (ID)','Latent (ID)','Location','northwest'); grid;
title([jdollars,' Hourly Filtered Sensible and Latent Heat: Turb vs Bulk']);

figure;
minx = min(hnetb_hr)-30; maxx = max(hnetb_hr)+30;
h4 = plot(jda,z_line2,'k',jda,hnetb_hr,'r',jda,hnet_hr,'b');
h4 = h4(2:end); % select traces for legend, omit zero line
xlabel('DOY'); ylabel('Hourly Net Heat Flux (W/m^2)'); xlim([stdt endt]);
ylim([minx maxx]); legend(h4, 'Bulk','Turb','Location','Best'); grid;
title([jdollars,' Hourly Net Heat Flux w/Respect to Sea']);

figure;
plot(yusb,yusr_psd,'r.',yusb,yusr_uh,'b.',yusb,yusr_id_psd,'g.',yusb,yusr_id_uh,'m.');
hold on; plot(yusb,yusb,'k-'); grid; axis([0,1.3,0,1.3]);
xlabel('ustar, bulk, m/s'); ylabel('ustar, various, m/s');
legend('cov psd','cov uh','id psd','id uh','location','northwest')
title([jdollars,' Hourly Friction Velocity vs Bulk']);

figure;
plot(yU,A(:,19),'k.',yU,ycd10n_id_psd,'r.',yU,ycd10n_id_uh,'b.'); grid;
xlabel('U m/s'); ylabel('Cd10N');axis([0,26,0,0.003])
legend('bulk','id psd','id uh','location','northwest')
title([jdollars,' Hourly Cd10N vs U']);

%% save all figures

if savefigs
    graphformat = '.png';
    graphdevice = '-dpng';
    figs = findall(0,'type','figure');
    sorted_figs = sortfighandlearray(figs); % put in ascending order
    N = length(sorted_figs);
    if N>0
        for ii=1:N
            figure(sorted_figs(ii));
            ppath = fullfile(fig_path,['da_figure_',sprintf('%02i',ii),graphformat]);
            print(graphdevice,ppath);
        end;
    end;
end;

