% DYNAMO_2011_da_red1.m - process output da routines


%% initialize
disp('DYNAMO_2011_da_red1_v3')
close all;
clear;

%% set run variables
prtit = true;     % set to 1 to save flux result files to disk
savefigs = true;  % set to save figures to disk

datasetv = 'v3';    % expanded da file and full decorr of Tsonic

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

% set start julian date for plots
stjd = 275;

% path to 'Processed' folder with allcruise da file
da_path = '/Users/byronb/Raw_Data_Files.noindex/DYNAMO_2011/Revelle/flux/Processed_new/leg2';       % for Mac file system
% path for saving da figures
fig_path = '/Users/byronb/Raw_Data_Files.noindex/DYNAMO_2011/Revelle/flux/Processed_Images_new/leg2/da_red_v3';
mkdir(fig_path);  % create directory if it doesn't exist

% read da, cr files with 10-min results
x = load(fullfile(da_path,['DYNAMO_2011_da_',datasetv,'_',dcname,'_c25.txt']));       % 164 (v3) column array
yxx = load(fullfile(da_path,['DYNAMO_2011_cr_',datasetv,'_',dcname,'_c25.txt']));
% yxx columns
%         COLUMN  VARIABLE        DESCRIPTION
%         1       doy             decimal day of year
%         2       H2O_Lic_all     mean H2O, all data, mmol/m3
%         3       H2O_Lic_good    mean H2O, good data, mmol/m3
%         4       CO2_Lic_all     mean CO2, all data, mmol/m3
%         5       CO2_Lic_good    mean CO2, good data, mmol/m3
%         6       P_Licor         atm pressure from Licor box, kPa
%         7       T_Licor         air temperature from Licor box, C
%         8       wq              water vapor flux
%         9       hlc             latent heat flux = wq*Rho*Cp
%         10      not_used1       CO2 covariance, not computed
%         11      not_used2       CO2 covariance, not computed
%         12      std_H2O_all     std H2O, all data
%         13      std_H2O_good    std H2O, good data
%         14      std_CO2_all     std CO2, all data
%         15      std_CO2_good    std CO2, good data
%         16      P_std           std Licor box pressure
%         17      Licor_AGC       mean diagnostic value

% load all cruise mean met file
MMF = load('/Users/byronb/Raw_Data_Files.noindex/DYNAMO_2011/Revelle/flux/Processed_old/Revelle10minutesLeg2_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;

% 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 =   15;         % SCS pressure sensor height

% PSD data
U       =   x(:,16);    % true wind from PSD sonic - relative to earth 10/31/2012 - LB
dir     =   x(:,41);    %
urel    =   x(:,38);    % relative wind speed, m/s
reldir  =   x(:,42);    % relative wind direction, deg clockwise rel ship's bow (PSD sonic)
reldir(reldir>180) = reldir(reldir>180)-360; %convert to +/- 180 format
tsnk    =   x(:,20);    % sea snake SST
ta      =   x(:,17);    % vaisala air temperature
rh      =   x(:,18);    % vaisala rh
rs      =   x(:,21);    % downward solar flux, W/m^2 (PSD units) - avg of psp1 & psp2
sig_rl  =   x(:,95);    % stdev downward solar flux, W/m^2 (PSD units)
rl      =   x(:,22);    % downward IR flux, W/m^2 (PSD units) - avg of pir1 & pir2
sig_rs  =   x(:,96);    % stdev downward IR flux, W/m^2 (PSD units)
org     =   x(:,24) - min(x(:,24));    % rainrate, mm/hr (org rain gauge, uncorrected)
zu_etl  =   ones(size(U_scs))*18.16;    % measurement heights
zt_etl  =   ones(size(U_scs))*16.5;
zq_etl  =   ones(size(U_scs))*16.5;
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

%------ Bulk Results from motcorr 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
Hrain   =   x(:,91);    % rain heat flux, W/m^2 as per Gosnell et. al.
ustr0   =   x(:,74);	% u-star from initial bulk run
tstr0   =   x(:,79);	% t-star
qstr0   =   x(:,80);	% q-star

%------ Navigation & data quality indicators
Lat     =   x(:,151);   % latitude, deg, SCS
Lon     =   x(:,150);   % longitude, deg, SCS
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
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
badSon  =   x(:,159);   % number of NaNs in sonic data in each 10-min interval
Ts_noise =  x(:,164);   % T sonic noise variance

% variances and covariances, turbulence data
wuj      =   x(:,117);   % wu cov
wvj      =   x(:,119);   % wv cov
wt_son   =   x(:,121);   % wTsonic cov
uuj      =   x(:,125);   % u var
vvj      =   x(:,126);   % v var
wwj      =   x(:,127);   % w var
tt       =   x(:,128);   % Tsonic var

sp1      =   x(:,129);   % ship motion stdev, u direction (m/s)
sp2      =   x(:,130);   % ship motion stdev, v direction (m/s)
sp3      =   x(:,131);   % ship motion stdev, w direction (m/s)

% select inertial dissipation results
cu       =   x(:,134);   % cu^2 U method a
cw       =   x(:,136);   % cw^2 W method a
ct       =   x(:,138);   % ct^2 temperature method a
cql      =   x(:,140);   % cq^2 H2O Licor method a
usi      =   x(:,142);   % inertial dissipation scaling paramters, all method a
tsi      =   x(:,144);   %
qsi      =   x(:,146);   %
% cu       =   x(:,135);   % cu^2 U method b
% cw       =   x(:,137);   % cw^2 W method b
% ct       =   x(:,139);   % ct^2 temperature method b
% cql      =   x(:,141);   % cq^2 H2O Licor method b
% usi      =   x(:,143);   % inertial dissipation scaling paramters, all method b
% tsi      =   x(:,145);   %
% qsi      =   x(:,147);   %

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

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

% ocean current components
iii = find(MMF.yday>=275 & MMF.yday<279);
CN = -MMF.cspd(iii)'.*cos((MMF.cdir(iii)')./180.*pi);  %'toward' North
CE = -MMF.cspd(iii)'.*sin((MMF.cdir(iii)')./180.*pi);  %'toward' East

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

% compute U and dir relative to water
UN = U.*cos(dir./180.*pi); % U North
UE = U.*sin(dir./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))/(pi/180)+360,360);  % Wind direction (from) relative to water

% figure; plot(jdy,U,'k',jdy,Ul,'r');
% figure; plot(jdy,dir,'k',jdy,dirl,'r');

Uearth = U;          % save copies U & dir in earth coords
direarth = dir;
U = Ul; dir = dirl;  % U and dir now relative to water

% save adjusted wind speeds and original wind speeds
x(:,164) = Uearth;   % U in earth coords appended to da array
x(:,165) = direarth; % dir in earth coords appended to da array
x(:,16)  = U;        % U is now relative to water
x(:,41)  = dir;      % dir is now relative to water

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

% best bulk met values specified here...
u_best = U;
ts_best = tsnk;
qs_best = qse;
zts = ones(length(jdy),1)*0.05;	% 0.05m depth of seasnake on Revelle
ta_best = ta;
rh_best = rh;
qa_best = qa;
rs_best = rs;
rl_best = rl;
P_best = press;
rain_best = org;
rhoair = P_best.*100./(Rgas*(ta_best+tdk).*(1+0.61*qa_best/1000));
Le = (2.501-.00237*nanmedian1(ts_best))*1e6;

% set measurement heights
zu = zu_etl(1);      % wind speed measurement height (m)
zt = zt_im(1);       % air T measurement height (m)
zq = zq_im(1);      % air q measurement height (m)

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

% scalar heat fluxes from measurements
wt = wt_son+0.51*(ta+tdk).*ustr0.*qstr0;  % correct sonic wt for humidity contribution
hsc = wt*cpa.*rhoair;
hsc(hsc<-100) = NaN;
hsc_son = wt_son*cpa.*rhoair;
hlc_lic = wq_lic*Le.*rhoair/1000;

% plot for map of cruise track: m_map folder must be in MATLAB path!
Lonmin = 68;   % longitudes are 0-360, going E
Lonmax = 110;
Latmin = -15;    % deg N is positive
Latmax = 10;
LatLonLim = [Lonmin Lonmax Latmin Latmax];
map_DYNAMO_2011(Lon,Lat,'DYNAMO_2011',LatLonLim);

% compute clear sky downwelling longwave flux
rlcler = 5.67e-8*(ta+273.16).^4.*(.52+.085*sqrt(qa));   % model clear sky
rlclr = (.52+.13/60*abs(Lat)+(.082-.03/60*abs(Lat)).*sqrt(qa)).*(5.67e-8*(ta+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;
DYNAMO_2011_Rsclear;     % runs SolarRadiance model

% COARE options
jcool = 1;
jwarm = 1;
icount = 1;

% housekeeping variables
qcol_ac = 0;
qjoule = 0;
tau_ac = 0;
jtime = 0;
jamset = 0;
tau_old = .06;
hs_old = 10;
hl_old = 100;
dsea = 0;
dt_wrm = 0;
tk_pwp = 0;
fxp = .7;
ts_depth = zts(1);  % bulk water temperature sensor depth

%% major read data loop
[nn,nx] = size(x');         % lines of 10min data
ijx = 1;
d1 = datenum('31-Dec-2010'); % t0 2011 julian dates
clear dt;

% major read loop for bulk model run
U10N = NaN(nx,1); % for saving 10m neutral values
T10N = U10N; Q10N = U10N;

xsave = NaN(nx,19);
for ibg = 1:nx
    % set variables not in data base
%   P=1010;     % air pressure
    us = 0;       % surface current
    zi = 600;     % inversion ht

    % decode date
    [iyr,mon,iday,ihr,imin,isec] = datevec(jdy(ibg)+d1);
    timex = [iyr mon iday ihr imin isec];
    jd = jdy(ibg);

    % decode bulk met data
    u = u_best(ibg);
    phirel = reldir(ibg);  % relative wind direction
    tsea = ts_best(ibg);   % bulk sea surface temp
    t = ta_best(ibg);      % air temp
    P = P_best(ibg);       % bp mb
    qs = qs_best(ibg);     % bulk sea surface humidity
    q = qa_best(ibg);      % air humidity
    Rs = rs_best(ibg);     % downward solar flux
    Rl = rl_best(ibg);     % downward IR flux
    rain = rain_best(ibg); % rain rate
    if sig_h(ibg)<5
       jt=0;               % flux data quality indicator
    else
       jt=1;
    end;
    tilt = tiltx(ibg);     % flow tilt, fixed frame
    jm = sig_u(ibg);       % ship manauver index
    lonn = Lon(ibg);       % longitude
    latx = Lat(ibg);       % latitiude

    % variables for warm layer
    ntime = 1e6*mon+1e4*iday+100*ihr+imin;
    time = (ihr*3600+imin*60)/24/3600;
    intime = time;
    loc = (lonn+7.5)/15;
    Rnl = .97*(5.67e-8*(tsea+273.16)^4-Rl);
    Rns = .945*Rs;

    % set condition dependent stuff
    Le = (2.501-.00237*tsea)*1e6;
    cpv = cpa*(1+0.84*q/1000);
    rhoa = P*100/(Rgas*(t+tdk)*(1+0.61*q/1000));
    visa = 1.325e-5*(1+6.542e-3*t+8.301e-6*t*t-4.8e-9*t*t*t);
    Al = 2.1e-5*(tsea+3.2)^0.79;

    %% apply warm layer
    %  can exclude this section by setting jwarm = 0 above...
    if jwarm == 1;
        chktime = loc+intime*24;
        newtime = (chktime-24*fix(chktime/24))*3600;
        if icount>1                     % not on first time thru
            if newtime>21600 && jump==1
                %goto 16
            else
                if newtime < jtime      %re-zero at midnight
                    jamset = 0;
                    fxp = 0.5;
                    tk_pwp = 19;
                    tau_ac = 0;
                    qcol_ac = 0;
                    dt_wrm = 0;
                    jump = 0;           % goto 16
                else                    % set warm layer constants
                    rich = 0.65;        %crit rich
                    ctd1 = sqrt(2*rich*cpw/(Al*grav*rhow));
                    ctd2 = sqrt(2*Al*grav/(rich*rhow))/(cpw^1.5);

                    dtime = newtime-jtime;              % delta time for integrals
                    qr_out = Rnl+hs_old+hl_old+RF_old;  % total cooling at surface
                    q_pwp = fxp*Rns-qr_out;             % tot heat abs in warm layer
                    if q_pwp<50 && jamset==0            % check for threshold
                        %goto 16
                    else
                        jamset = 1;                     % indicates threshold crossed
                        tau_ac = tau_ac+tau_old*dtime;  % momentum integral
                        if qcol_ac+q_pwp*dtime>0        % check threshold for warm layer existence
                            for i=1:5                   % loop 5 times for fxp
                                fxp = 1-(0.28*0.014*(1-exp(-tk_pwp/0.014))+0.25*0.357*(1-exp(-tk_pwp/0.357))+0.47*12.82*(1-exp(-tk_pwp/12.82)))/tk_pwp;
                                % fg = fpaul(tk_pwp); fxp = fg(1);
                                qjoule = (fxp*Rns-qr_out)*dtime;
                                if qcol_ac+qjoule>0
                                    tk_pwp = min(19,ctd1*tau_ac/sqrt(qcol_ac+qjoule));
                                end;
                            end; % end i loop
                        else                            % warm layer wiped out
                            fxp = 0.75;
                            tk_pwp = 19;
                        end;                            % end sign check on qcol_ac
                        qcol_ac = qcol_ac+qjoule;       % heat integral
                        % compute dt_warm
                        if qcol_ac>0
                            dt_wrm = ctd2*(qcol_ac)^1.5/tau_ac;
                        else
                            dt_wrm = 0;
                        end;
                    end;            % end threshold check
                end;                % end midnight reset
                if tk_pwp<ts_depth
                    dsea = dt_wrm;
                else
                    dsea = dt_wrm*ts_depth/tk_pwp;
                end;
            end;                    % end 6am start first time thru
        else
            jump=1;                 % indicates data start period
        end;                        % end first time thru check
        jtime=newtime;
    end;                            % end warm layer model appl check

    %% run COARE for 10-min bulk model estimates
    ts = tsea + dsea; % adjust SST for warm layer
    qs = qsea(ts);    % compute qs at new SST
    xx = [u us ts t qs q Rs Rl rain zi  P zu zt zq latx jcool 0 0 0];
    
    % COARE 3.0
    [y, u10n, t10n, q10n] = cor30a(xx);

    % COARE 3.5 - vectorized, but here running point-by-pioint within loop
    %    to conform with output from cor30a above...
% 	[y, u10n, t10n, q10n] = coare35vn(u,zu,t,zt,rh_best(ibg),zq,P,ts,Rs,Rl,latx,zi,rain,NaN,NaN);

    %% output from COARE
    hsb    =    y(1);       % sensible heat flux W/m/m
    hlb    =    y(2);       % latent heat flux W/m/m
    tau    =    y(3);       % stress N/m/m
    zo     =    y(4);       % vel roughness length
    zot    =    y(5);       % temp "
    zoq    =    y(6);       % hum  "
    L      =    y(7);       % Obukhov Length
    usr    =    y(8);       % ustar
    tsr    =    y(9);       % tstar
    qsr    =    y(10);      % qstar  [g/g]
    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);      %
    % zax is expanded output array
    zax(1) =    jd;         % julian day
    zax(2:10)=  xx(1:9);    %
    zax(4) =    ts;         % Tsea [no cool skin]
    zax(11:32)= y;          %
    zax(33:35)= [dt_wrm tk_pwp ts];  % warm layer deltaT, thickness, corrected Tsea
    % previous values from cwf hp basic code
    hsx    =    Hsb(ibg);
    hlx    =    Hlb(ibg);
    taux   =    Taub(ibg);

    %% save various parts of data
    dt(ibg,:) = zax;
    % dt description: this array used for inertial diss processing later
    % [jd  u  us  Tsea[no cool skin]  t qs q Rs  Rl rain 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 dt_wrm tk_pwp  ts]
    %   1  2   3          4           5 6  7  8  9  10   11  12  13  14 15  16  17 18  19  20  21   22   23  24 25   26  27  28   29     30     31   32   33      34   35

    U10N(ibg,1) = u10n; % save 10m bulk model output and adjusted qs
    T10N(ibg,1) = t10n;
    Q10N(ibg,1) = q10n;

    eyes(ibg,:) = [jt tilt jm phirel];
    e = [num2str(mon),'  ',num2str(iday),'  ',num2str(ihr),'   ',num2str(dter),'  ',num2str(dt_wrm),'  ',num2str(tk_pwp)];
%     disp(e)
    hs_old = hsb;
    hl_old = hlb;
    RF_old = RF;
    tau_old = tau;
    icount = icount+1;
    if phirel<90 && phirel>-90  % save subset with good wind direction
        dtx(ijx,:) = zax;
        ijx = ijx+1;
    end;
end; % end major read loop


%% 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);
ii = find(isfinite(tsnk));
Le = (2.501-.00237*median(tsnk(ii)))*1e6;

% 11/01/2012 - LB - recompute hsc here as u* and q* have been updated by new coare run
ustr0 = dt(:,18);
tstr0 = dt(:,19);
qstr0 = dt(:,20);
wt = wt_son+0.51*(ta+tdk).*ustr0.*qstr0;  % correct sonic wt for humidity contribution
hsc = wt*cpa.*rhoair;
hsc(hsc<-100) = NaN;
hsc_son = wt_son*cpa.*rhoair;

% compute hl_webb correction from wbar_webb (dt(:,25))
hl_webb = dt(:,25).*rhoair*Le/1000.*qa;

% select best covariance latent heat
% ii = find(sgq_lic>.7 | hlc_lic<-20); hlc_lic(ii) = NaN;
% hlc_bl = hlc_oph; % (hlc_oph+hlc_lic)/2;
hlc_bl = hlc_lic; %  Licor only since 2006.

%% heat flux QC plots
lapse = grav/cpa; Theta10 = T10N + lapse*10; % Pot temp @ 10m
qs_adj = qs_best - dt(:,22)*1000; % qs - dqer
ts_adj = ts_best - dt(:,21); % ts - dter

wt_U = wt./U10N;
dTs = ts_adj-Theta10;
idx = (isfinite(wt_U) & isfinite(dTs) & dTs<3 & dTs>-2); %wT plot criteria
minx = min(dTs(idx)); maxx = max(dTs(idx));
coef = polyfit(dTs(idx),wt_U(idx),1); liney = coef(1)*[minx,maxx] + coef(2); xIntcpt = coef(2)/-coef(1);
figure; plot(dTs(idx),wt_U(idx),'r.',[minx,maxx],liney,'b'); ylim([-0.01,0.01]); grid;
hold on; plot([0,0],[-0.01,0.01],'k',[minx,maxx],[0,0],'k'); xlim([minx,maxx]);
ylabel('wT/U_{10n}'); xlabel('(Ts - \Theta_{10})');
title(['Ch10N x 1000 = ',num2str(coef(1)*1000,'%6.5f'),'   y-Intcpt = ',num2str(coef(2),'%6.5f'),'   x-Intcpt = ',num2str(xIntcpt,'%6.5f')]);

wq_U = (wq_lic-(dt(:,25).*dt(:,7)))./U10N; % (wq - wbar*qbar)/u10n
dQs = qs_adj-Q10N;
idx = (isfinite(wq_U) & isfinite(dQs)); %wq plot criteria
minx = min(dQs(idx)); maxx = max(dQs(idx));
coef = polyfit(dQs(idx),wq_U(idx),1); liney = coef(1)*[minx,maxx] + coef(2); xIntcpt = coef(2)/-coef(1);
figure; plot(dQs(idx),wq_U(idx),'r.',[minx,maxx],liney,'b'); ylim([-0.004,0.01]); grid;
hold on; plot([0,0],[-0.01,0.01],'k',[minx,maxx],[0,0],'k'); xlim([minx,maxx]);
ylabel('wq/U_{10n}'); xlabel('(Qs - Q_{10N})');
title(['Ce10N x 1000 = ',num2str(coef(1)*1000,'%6.5f'),'   y-Intcpt = ',num2str(coef(2),'%6.5f'),'   x-Intcpt = ',num2str(xIntcpt,'%6.5f')]);

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

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

%% ocean warm layer profiles
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-hlc_bl-hl_webb-HRF;       % net heat using cov fluxes hsc & hlc
hnet_id = 0.955*rs_best-rnl-hsib-hlib_lic-hl_webb-HRF;   % net heat using id fluxes hsib & hlib_lic

%% compute daily averages, do cloud plots

DYNAMO_2011_daily;

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

%% now compute 1-hr averages & bulk model
DYNAMO_2011_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
DYNAMO_2011_hr_bulk_c30;
% DYNAMO_2011_hr_bulk_c35;

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

% hourly avg stress components
taucx = ytauc;          % covariance streamwise
taucy = ytaucy;			% covariance x-stream
tauib = ytauib;			% inertial dissipation
taub = ytaub;			% coare model result

% checking best median values of stress from measurements and bulk
med_tau_plt_DYNAMO_2011;

%% Hourly net heat flux - using 'best' rl & rs
% net longwave relative to sea from warm-layer-corrected Ts (up flux is negative)
yrl = y4(:,9);
rnl_hr = 0.97*(5.67e-8*(yts + 273.15).^4 - yrl);
yrs = y4(:,8);
% hnet from COARE 3.5 using ship bulk met. y1(:,8) is ship PIR.
hnetb_hr = 0.955*yrs - rnl_hr - yhsb - yhlb - yhrain;
% net heat using covariance hl
hnet_hr = 0.955*yrs - rnl_hr - yhsc - (hlc_licbar + hl_webar) - yhrain;

%% save hourly flux to file if desired
if prtit
    prt_DYNAMO_2011_flux_hr;
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
jjj = find((reldir<90 & reldir>-90) & sig_h<5 & sig_u<0.6 & sp2<0.8 & rain_best<5 & sqrt(wwj)./ugw<2+.0015*u_best.^2 & sqrt(vvj)./ugu<2 & badSon<50);

% various plots of 10-min and hourly results
z_line = jdy*0;
figure; 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 (2014)'); ylabel('Stremwise stress (N/m2)');
grid; xlim([stjd (jdy(end)+1)]); ylim([-0.2 0.8]);
title([jdollars,' Stress: Filtered 10-min Cov, ID and COARE']);

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

figure; plot(jdy(jjj),hlc_bl(jjj),'bo',jdy(jjj),hlib_lic(jjj),'m+',jdy(jjj),hlbb(jjj),'g.',jdy,z_line,'k-')
legend('Cov','ID','Bulk','Location','Best');
grid; xlim([stjd (jdy(end)+1)]); ylim([-50 400]);
xlabel('DOY (2014)'); ylabel('Latent heat flux (W/m2)');
title([jdollars,' Latent Heat: Filtered 10-min Cov, ID and COARE']);

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

z_line2 = jda*0;
figure;
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 (2014)'); ylabel('Heat Flux Component (W/m^2)'); xlim([stjd (jda(end)+1)]);
ylim([-400 1000]); legend(h2, 'net IR','Solar','Bulk Hs','Bulk Hl','Bulk RF','Location','Best'); grid;
title([jdollars,' Net Hourly Heat Flux Components from COARE']);

figure;
h3 = plot([-400 400],[-400 400],'k',[-400 400],[0 0],'k',[0 0],[-400 400],'k',yhsb,yhsc,'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([-80 400 -80 400]); legend(h3,'yhsc','yhlib-lic + hl-webar','Location','northwest'); grid;
title([jdollars,' Hourly Sensible and Latent Heat: Turb vs Bulk']);

figure;
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 (2014)'); ylabel('Hourly Net Heat Flux (W/m^2)'); xlim([stjd (jda(end)+1)]);
ylim([-400 800]); legend(h4, 'Bulk','Turb(cov)','Location','Best'); grid;
title([jdollars,' Hourly Net Heat Flux w/Respect to Sea']);

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

