% STRATUS_2017_da_red1.m - process output da routines

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

% matlab script path
% path_prog = 'c:\Users\Sergio\Documents\MATLAB\2017_STRATUS\da_red_analysis';
path_prog = '/Users/byronb/Documents/MATLAB/2017_STRATUS/da_red_analysis';
cd(path_prog);                  % cd to matlab directory
restoredefaultpath              % reset matlab path
addpath(genpath(path_prog));	% add scripts to path
savepath; rehash;

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

datasetv = 'v4';    % full decorr of Tsonic

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

% set limits for various filter variables to be used in all routines
rdir_lo     = -90;	% rel wind direction, deg from bow
rdir_hi     = 90;
rain_lim    = 1;    % rain rate, mm/hr
sp2_lim     = 0.8;  % vplat std dev  m/s
sigu_lim    = 1.5;  % ship speed std dev  m/s
sigh_lim    = 5;    % ship heading std dev, deg
badSon_lim  = 50;   % # bad sonic points per 10min interval
missing_lim = 50;   % # missing sonic points per 10min interval
ushp_lim    = 3.5;  % ship speed, m/s
tilt_lim    = 10;   % flow tilt angle, degrees
sig_rwd_lim = 15;   % std dev rwdir

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

% Track map limits
Lonmin = 240;
Lonmax = 310;
Latmin = -25;
Latmax = 33;
LatLonLim = [Lonmin Lonmax Latmin Latmax];

% as measured on CalWater 01/13/2015: waterline to base of tower = 7.16 m.
zu = 18.16;     % wind speed measurement height (m): 7.16m + 11m
zt = 16.5;      % air T measurement height (m)
zq = 16.5;      % air q measurement height (m)
zP = 17;        % psd pressure sensor height
zWXT = 17;      % WXT sensor height (m)
zu_ship = 15;   % ship SCS sensor heights
zt_ship = 15;
zq_ship = 15;
ztsg = 5;	    % 5m depth of tsg

tdk = 273.15;

%% read files with 10-min results

% path to 'Processed' folder with allcruise da files
da_path = '/Users/byronb/Raw_Data_Files.noindex/STRATUS_2017/RHB/flux/Processed_STRATUS_only';	% for Mac file system
% da_path = 'D:\Data\STRATUS_2017\RHB\flux\Processed';
fig_path = '/Users/byronb/Raw_Data_Files.noindex/STRATUS_2017/RHB/flux/Processed_Images/da_red_plots';
% fig_path = 'D:\Data\STRATUS_2017\RHB\flux\Processed_Images\da_red_plots';
mkdir(fig_path);  % create directory if it doesn't exist

% read da, cr files with 10-min results, 3312 x 173
%{
da file format for STRATUS 2017
COL     VARIABLE        DESCRIPTION
    1       doy             decimal day of year
    2:10    xship           COARE 3.5 input, scs:[Utrue,Tair,RH,P_mb,Tsea,
                            Rs,Rl,lat,rain]
    11:15   no_data         To maintain historical column numbering...
    16:24   xpsd            COARE 3.5 input, PSD, as above
    25:29   no_data         To maintain historical column numbering...
    30      u_ship          ship mean true wind speed (m/s)
    31      u_ship_min      ship minimum true wind speed (m/s)
    32      u_ship_max      ship maximum true wind speed (m/s)
    33      u_ship_std      ship true wind speed stdev (m/s)
    34      wdir_ship       ship mean true wind direction
    35      wdir_ship_min   ship minimum true wind direction
    36      wdir_ship_max   ship maximum true wind direction
    37      wdir_ship_std   ship true wind direction stdev
    38      rwspd_psd       streamwise relative wind speed (m/s)from PSD sonic
    39      Un              N component, corrected true windspeed, PSD sonic
    40      Ue              E component, corrected true windspeed, PSD sonic
    41      wdir_psd        true wind direction from N,E components, PSD sonic
    42      rwdir_psd       relative wind direction from PSD sonic
    43      azm             azimuth angle from first rotation, radians
    44      tilt            tilt angle from first rotation, radians
    45      Tsea_psd        seasnake SST, degC
    46      Tair_psd        air temperature, PSD vaisala, degC
    47      RH_psd          RH, PSD vaisala
    48      RH_ship         RH, ship sensor
    49      Q_psd           specific humidity, PSD sensors, g/kg
    50      Q_ship          specific humidity, ship sensors, g/kg
    51      Qsea_psd        sea surface specific humidity from seasnake SST
    52      rhoa_psd        moist air density from PSD sensors (kg/m3)
    53:73   yship           COARE 3.5 output, ship data:[usr,tau,hsb,hlb,hlwebb,
                         tsr,qsr,zot,zoq,Cd,Ch,Ce,L,zet,dter,dqer,tkt,RF,
                         Cdn_10,Chn_10,Cen_10]
    74:94   ypsd            COARE 3.5 output, PSD data, as above
    95      Rl_std_psd      std PIR flux, PSD
    96      Rs_std_psd      std PSP flux, PSD
    97      sog_ship        mean sog, ship gps, m/s
    98      sog_ship_std    std sog, ship gps
    99      sog_ship_min    minimum sog, ship gps
    100     sog_ship_max    maximum sog, ship gps
    101     cog_ship        mean cog, deg, ship gps
    102     cog_ship_std    std cog, ship gps
    103     cog_ship_min    minimum cog, ship gps
    104     cog_ship_max    maximum cog, ship gps
    105     cog_psd         mean cog, deg, PSD gps
    106     cog_psd_std     std cog, PSD gps
    107     cog_psd_min     minimum cog, PSD gps
    108     cog_psd_max     maximum cog, PSD gps
    109     sog_psd         mean sog, m/s, PSD gps
    110     sog_psd_min     minimum sog, PSD gps
    111     sog_psd_max     maximum sog, PSD gps
    112     lrg             mean heading, deg, ship gyro
    113     lrg_std         std heading, ship gyro
    114     lrg_min         minimum heading, ship gyro
    115     lrg_max         maximum heading, ship gyro
    116     no_data         wu, Hare correction, not computed
    117     wu              wu covariance, Edson correction
    118     no_data         wv, Hare correction, not computed
    119     wv              wv covariance, Edson correction
    120     no_data         wT, Hare correction, not computed
    121     wt              wT covariance, Edson correction
    122     no_data         uvar, Hare correction, not computed
    123     no_data         vvar, Hare correction, not computed
    124     no_data         wvar, Hare correction, not computed
    125     uvar            u variance, Edson correction
    126     vvar            v variance, Edson correction
    127     wvar            w variance, Edson correction
    128     tvar            T variance
    129     Sp1             U platform velocity std, Edson correction
    130     Sp2             V platform velocity std, Edson correction
    131     Sp3             W platform velocity std, Edson correction
    132     hd_sin_std      stdev sine lrg heading
    133     hd_cos_std      stdev cosine lrg heading
    134     Cu2a            Cu structure func parameter, method a
    135     Cu2b            Cu structure func parameter, method b
    136     Cw2a            Cw structure func parameter, method a
    137     Cw2b            Cw structure func parameter, method b
    138     Ct2a            Ct structure func parameter, method a
    139     Ct2b            Ct structure func parameter, method b
    140     Cq2a            Ch2o structure func parameter, method a
    141     Cq2b            Ch2o structure func parameter, method b
    142     usid_a          ustar from idiss, method a
    143     usid_b          ustar from idiss, method b
    144     tsid_a          tstar from idiss, method a
    145     tsid_b          tstar from idiss, method b
    146     qsid_a          qstar from idiss, method a, kg/kg
    147     qsid_b          qstar from idiss, method b, kg/kg
    148     Lid_a           L from idiss, method a
    149     Lid_b           L from idiss, method b
    150     lon             longitude, scs gps
    151     lat             latitiude, scs gps
    152     Ubar            streamwise true windspeed from Edson correction
    153     odec            odec (ship speed through water), m/s
    154     odec_std        stdev odec
    155     odec_min        minimum odec
    156     odec_max        maximum odec
    157     P_mb            atmospheric pressure, PSD sensor, mb
    158     missing         missing sonic wind data lines per 10 min interval
    159     badPnts         # NaNs in sonic data per 10 min inteerval
    160     hs_cov          sensible heat flux from covariance
    161     hl_cov          latent heat flux from covariance, uncorrected
    162     hlc_cov         webb corrected hl_cov
    163     sal_scs         salinity, ship system, ppt
    164     Ts_noise        Sonic temperature noise variance
    165     rwdir_psd_std   std dev relative wind direction
    166     rwdir_wxt       rel wind direction ,0-360 deg from bow, WXT
    167     rwspd_wxt       wind speed, m/s, WXT
    168     Tair_wxt        air temperature, C, WXT
    169     RH_wxt          RH, %, WXT
    170     Pmb_wxt         air pressure, mb or hPa, WXT
    171     rain_wxt        rain rate, mm/hr, WXT
    172     spd_wxt         true wind speed, WXT
    173     dir_wxt         true wind dir, WXT
%}
x = load(fullfile(da_path,['STRATUS_2017_da_',dcname,'.txt']));
yxx = load(fullfile(da_path,['STRATUS_2017_cr_',dcname,'.txt']));

% 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);     % avg port/stbd solar flux, W/m^2 (imet units)
rl_im   =   x(:,8);     % avg port/stbd ir flux, W/m^2 (imet units)
rain_im =   x(:,10);    % precip, mm/hr
zu_im   =   ones(size(U_scs))*zu_ship;    % measurement heights for ship sensors, m
zt_im   =   ones(size(U_scs))*zt_ship;
zq_im   =   ones(size(U_scs))*zq_ship;
rh_im   =   x(:,4);     % SCS RH
P_im    =   x(:,5);                         % ship is already sea level corrected SCS pressure 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

% 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, deg clockwise rel ship's bow (PSD sonic)
tsnk    =   x(:,20);    % sea snake SST, but after jd 144 tsnk = tsg
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_rs  =   x(:,96);    % stdev downward IR flux, W/m^2 (PSD units)
rl      =   x(:,22);    % downward IR flux, W/m^2 (PSD units) - avg of pir1 & pir2
sig_rl  =   x(:,95);    % stdev downward solar flux, W/m^2 (PSD units)
org     =   x(:,24);    % rainrate, mm/hr (org rain gauge, uncorrected)
press   =   x(:,19)+0.125*zP;       % sea level corrected 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
rwdir_wxt = x(:,166);   % rel wind direction ,0-360 deg from bow, WXT
rwspd_wxt = x(:,167);   % wind speed, m/s, WXT
Tair_wxt  = x(:,168);   % air temperature, C, WXT
RH_wxt    = x(:,169);   % RH, WXT
Pmb_wxt   = x(:,170)+0.125*zWXT;   % sea level air pressure, mb or hPa, WXT
rain_wxt  = x(:,171);   % rain rate, mm/hr, WXT
spd_wxt   = x(:,172);   % true wind speed, WXT
dir_wxt   = x(:,173);   % true wind dir, WXT
qa_wxt    = qair_p([Tair_wxt,RH_wxt,Pmb_wxt]);  % air specific humidity, g/kg

% 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
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
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(:,113);   % std dev of ship heading (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(:,164);   % T sonic noise variance
std_rwdir = x(:,165);   % std dev rel wind direction

% 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
switch IDflag
    case 'a'
        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);
    case 'b'
        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);
end

%% make adjustments, calibrations, etc.
% apply fixes to data and select best quality measurements
% for use in bulk flux calculations.  Edit this script to
STRATUS_2017_fixit_da1
% best bulk met values specified here...
u_best = U;
dir_best = dir;
if jdy<140
    ts_best = tsnk; ts_best(jdy>135.8 & jdy<136.5) = tsg(jdy>135.8 & jdy<136.5) - 0.7;
else
    ts_best = tsg;
end
ta_best = ta;
rh_best = RH_wxt;
rs_best = rs;
rl_best = rl;
P_best = Pmb_wxt;
rain_best = org;

qa_best = qair_p([Tair_wxt,RH_wxt,Pmb_wxt]);  % recompute q  from WXT measurements;
qs_best = qsea_p([ts_best P_best]);
rhoair = P_best.*100./(Rgas*(ta_best+tdk).*(1+0.61*qa_best/1000));
Le = (2.501-.00237*nanmedian1(ts_best))*1e6;

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

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

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

%% constants and other variables for the bulk flux calculations
grav = 9.72;
Rgas = 287.1;
cpa = 1004.67;
be = 0.026;
cpw = 4000;
rhow = 1022;
visw = 1e-6;
tcw = 0.6;

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

% 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!
map_STRATUS_2017(Lon,Lat,'STRATUS_2017',LatLonLim);

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

% shortwave clear sky radiative flux
js = stdt; je = endt;
iv = NaN(length(jdy),1);    % Rsclear sets to qa/qrat
min_set = 5;
STRATUS_2017_Rsclear;     % runs SolarRadiance model

% COARE options
jcool = 1;
if jdy<140
    jwarm = 0;  % no run warm layer run for STRATUS
else
    jwarm = 1;
end
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;

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

% major read loop for bulk model run
U10N = NaN(nx,1); % for saving 10m neutral values
T10N = U10N; Q10N = U10N;
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
    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<ztsg
                    dsea = dt_wrm;
                else
                    dsea = dt_wrm*ztsg/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
    if jdy<140
        ts = tsea;  % not adjusting warm layer for STRATUS
    else
        ts = tsea + dsea; % adjust tsg for warm layer
        qs = qsea(ts);    % compute qs at new SST
    end

    % COARE 3.0
    xx = [u us ts t qs q Rs Rl rain zi  P zu zt zq latx jcool 0 0 0];
%     [y, u10n, t10n, q10n] = cor30a(xx);

    % COARE 3.5 - normally vectorized, but here running point-by-pioint within loop
    %    y output conforms with output from cor30a above...
	[y,u10n,t10n,q10n,u03n,t03n,rh03] = 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;
    U03N(ibg,1) = u03n; % save 10m bulk model output and adjusted qs
    T03N(ibg,1) = t03n;
    RH03(ibg,1) = rh03;

    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(ts_best));
Le = (2.501-.00237*median(ts_best(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_best+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))-0.5; maxx = max(dTs(idx))+0.5;
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(['Ch_{10N}x10^{-3} = ',num2str(coef(1)*1000,'%4.3f'),', y-Int = ',num2str(coef(2),'%3.2e'),', x-Int = ',num2str(xIntcpt,'%4.3f')]);

wq_U = (wq_lic-(dt(:,25).*dt(:,7)))./U10N; % webb corrected: (wq - wbar*qbar)/u10n
dQs = qs_adj-Q10N;
idx = isfinite(wq_U) & isfinite(dQs) & wq_U<4*nanstd1(wq_U); %wq plot criteria
minx = min(dQs(idx))-1; maxx = max(dQs(idx))+1;
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.01,0.015]); grid;
hold on; plot([0,0],[-0.015,0.015],'k',[minx,maxx],[0,0],'k'); xlim([minx,maxx]);
ylabel('wq/U_{10n}'); xlabel('(Qs - Q_{10N})');
title(['Ce_{10N}x10^{-3} = ',num2str(coef(1)*1000,'%4.3f'),', y-Int = ',num2str(coef(2),'%3.2e'),', x-Int = ',num2str(xIntcpt,'%4.3f')]);

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

%% inertial dissipation processing
STRATUS_2017_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 avgs & save results

STRATUS_2017_daily; % cloud plots

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

%% now compute 1-hr averages & rerun bulk model
STRATUS_2017_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
STRATUS_2017_hr_bulk_c35;

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

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

med_tau_plt_STRATUS_2017; % stress plots

% Hourly net heat flux
% 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 ID hl/hs
hnet_hr = 0.955*yrs - rnl_hr - yhsib - (yhlib_lic + hl_webar) - yhrain;

if prtit
    prt_STRATUS_2017_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<rdir_hi & reldir>rdir_lo) & sig_h<sigh_lim & sig_u<sigu_lim & std_rwdir<sig_rwd_lim &...
            sp2<sp2_lim & (sqrt(wwj)./ugw < 2+.0015*u_best.^2) & (sqrt(vvj)./ugu < 2) & missing<missing_lim);

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

figure;
slo = min([min(-(rhoair(jjj).*wuj(jjj))),min((rhoair(jjj).*usib(jjj).^2)),min(taubb(jjj))])-0.1;
shi = max([max(-(rhoair(jjj).*wuj(jjj))),max((rhoair(jjj).*usib(jjj).^2)),max(taubb(jjj))])+0.1;
h0 = 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.');
h0 = h0(2:end);   % omit zero line from legend
legend(h0,'Cov','ID','Bulk','Location','Best');
xlabel('DOY (2014)'); ylabel('Stremwise stress (N/m2)');
grid; xlim([stdt endt]); ylim([slo,shi]);
title([jdollars,' Stress: Filtered 10-min Cov, ID and COARE']);

figure;
hslo = min([min(hsc(jjj)),min(hsib(jjj)),min(hsbb(jjj))])-25;
hshi = max([max(hsc(jjj)),max(hsib(jjj)),max(hsbb(jjj))])+25;
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([stdt endt]); ylim([hslo,hshi])
xlabel('DOY (2014)'); ylabel('Sensible heat flux (W/m2)');
title([jdollars,' Sensible Heat: Filtered 10-min Cov, ID and COARE']);

hllo = min([min(hlc_bl(jjj)),min(hlib_lic(jjj)),min(hlbb(jjj))])-50;
hlhi = max([max(hlc_bl(jjj)),max(hlib_lic(jjj)),max(hlbb(jjj))])+50;
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([stdt endt]); ylim([hllo,hlhi]);
xlabel('DOY (2014)'); ylabel('Latent heat flux (W/m2)');
title([jdollars,' Latent Heat: Filtered 10-min Cov, ID and COARE']);

figure;
hnlo = min(min(hnet(jjj)),min(hnet_id(jjj)))-50;
hnhi = max(max(hnet(jjj)),max(hnet_id(jjj)))+50;
h1 = plot(jdy,z_line,'k-',jdy(jjj),hnet(jjj),'g.',jdy(jjj),hnet_id(jjj),'m+');
h1 = h1(2:end); % select traces for legend, omit zero line
xlabel('DOY (2014)'); ylabel('Net heat flux (W/m^2)'); xlim([stdt endt]);
ylim([hnlo,hnhi]); legend(h1, 'bulk hs&hl','ID hs&hl','Location','Best'); grid;
title([jdollars,' Filtered 10-min Net Heat Flux']);

% hourly plots
z_line2 = jda*0;
hflo = min([min(-rnl_hr),min(-yhsb),min(-yhlb)])-50;
hfhi = max((0.95*y1(:,7)))+50;
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([stdt endt]);
ylim([hflo,hfhi]); legend(h2, 'net IR','Solar','Bulk Hs','Bulk Hl','Bulk RF','location','north'); grid;
title([jdollars,' Net Hourly Heat Flux Components from COARE']);

figure;
hlo = min([min(yhsb),min((yhsib + yhsc)/2),min(yhlb),min((yhlib_lic + hl_webar))])-50;
hhi = max([max(yhsb),max((yhsib + yhsc)/2),max(yhlb),max((yhlib_lic + hl_webar))])+50;
h3 = plot([hlo,hhi],[hlo,hhi],'k',[hlo,hhi],[0 0],'k',[0 0],[hlo,hhi],'k',yhsb,(yhsib + yhsc)/2,'r.',yhlb,(yhlib_lic + hl_webar),'b.');
h3 = h3(4:end); % select traces for legend, omit zero & 1:1 lines
xlabel('Bulk Heat Flux (W/m^2)'); ylabel('Turb Heat Flux (W/m^2)');
axis([hlo,hhi,hlo,hhi]); legend(h3,'Sensible (cov+ID)/2','Latent (ID+webb)','Location','northwest'); grid;
title([jdollars,' Sensible and Latent Heat: Turb vs Bulk, Hourly Avg']);

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([stdt endt]);
ylim([hnlo,hnhi]); 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

