% HiWinGS2013_da_red1.m - process output da routines

disp('HiWinGS2013_da_red1')
% comments on COARE version 3.0
%{
toga coare bulk flux model version 3.0
***************************************
uses following matlab subroutines:
	julday.m
	cor30a.m
	psiu_30.m
	psit_30.m
	grv.m

*********** basic specifications  *****
	zu=			height of wind measurement
	zt=			height of air temperature measurement
	zq=			height of air humidity measurement
	ts_depth	depth of water temperature measurement
	jwarm=		0=no warm layer calc, 1 =do warm layer calc
	jcool=		0=no cool skin calc, 1=do cool skin calc

***********   input data **************
	jdy=		date in decimal julian day format
	lon=		longitude (deg e=+)
	u=			wind speed (m/s), height zu
	us=			surface current (m/s)
	ts=			bulk surface sea temp (cent)
	t=			air temp (cent), height zt
	qs=			sea surface sat specific humidity (g/kg)
	q=			air specific humidity (g/kg), height zq
	Rs=			downward solar flux (w/m^2)
	Rl=			downward IR flux (w/m^2)
	zi=			inversion height (m)
	P=			air pressure (mb)
	rain=		rain rate (mm/hr)
	lat=		lat

********** output data  ***************
	hsb=        sensible heat flux (w/m^2)
	hlb=		latent heat flux (w/m^2) inlcuding webb correction
	RF=			rain heat flux(w/m^2)
	W_webb=		mean webb velocity (m/s)
	tau=		stress (nt/m^2)
	zo=			velocity roughness length (m)
	zot			temperature roughness length (m)
	zoq=		moisture roughness length (m)
	L=			Monin_Obukhov stability length
	usr=		turbulent friction velocity (m/s), including gustiness
	tsr			temperature scaling parameter (K)
	qsr			humidity scaling parameter (g/g)
	dter=		cool skin temperature depression (K)
	dqer=		cool skin humidity depression (g/g)
	tkt=		cool skin thickness (m)
	Cd=			velocity drag coefficient at zu, referenced to u
	Ch=			heat transfer coefficient at zt
	Ce=			moisture transfer coefficient at zq
	Cdn_10=		10-m velocity drag coeeficient, including gustiness
	Chn_10=		10-m heat transfer coeeficient, including gustiness
	Cen_10=		10-m humidity transfer coeeficient, including gustiness
%}

%% initialize
close('all');
clear all;
didread=0;  % set to 1 to skip file load
prtit=1;

%selects the desired dataset version
datasetv ='v3';
decorr = 1;
dcname = [];
if decorr==1
    dcname ='vdecorr';
end;

% set local path to procturb folder with concatenated da text file
% and QA/Qc'd 10 minute mean met data (mat file).
procturb_path = '/Volumes/User/bblomquist/Documents/Cruises/HiWinGS_2013/Analysis_BWB/procturb_v3_decorr_Oct2014/total/';
met_path = '/Volumes/User/bblomquist/Documents/Cruises/HiWinGS_2013/MeanNBulkfluxes/Matlab_format/';

%% load data  HiWinGS_2013_da_v3.txt
if didread~=1;
    % read file with 10 min average da data
    x = load([procturb_path,'HiWinGS_2013_da_',datasetv,'_',dcname,'.txt']);
    % read files with licor data
    yxx = load([procturb_path,'HiWinGS_2013_cr_',datasetv,'_',dcname,'.txt']);
    yco = load([procturb_path,'HiWinGS_2013_co_',datasetv,'_',dcname,'.txt']);
    % QC mean met data - will use this as MET inputs for COARE
    MMF = load([met_path,'HiWinGS_2013_10minMeans_r2.mat']);
    % da file using UH wind/motion
% 	xUH=load([procturb_path,'HiWinGS_2013_daUH' datasetv dcname '.txt']);
      
%% define variables   
    jdy = x(:,1); % decimal julian day at start of 10-min averages
    jdy = round(jdy*1440)/1440;
    MMF.Yday = MMF.Yday - 5/1440;
    MMF.Yday = round(MMF.Yday*1440)/1440;    % same format as jdy now
    
    % find time indices common to both data sets
    % c vector is array of common times, c=jdy(ia) & c=MMF.yday10m(ib)
    [c, ia, ib] = intersect(jdy,MMF.Yday);
    
    % ships 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(:,4);     % tsg water temperature, C (5 m depth)
    ta_im   =   x(:,5);     % imet air temperature, C 
    qs_tsg  =   x(:,6);     % imet sea surface specific humidity, g/kg 
    qa_im   =   x(:,7);     % imet air specific humidity, g/kg
    rs_im   =   x(:,8);     % downward solar flux, W/m^2 (imet units)
    rl_im   =   x(:,9);     % downward ir flux, W/m^2 (imet units)
    rain_im =   x(:,10);    % pricip, mm/hr
    zu_im   =   x(:,13);    % measurement heights, m
    zt_im   =   x(:,14);    % 
    zq_im   =   x(:,15);    % 

	% PSD data
    U       =   x(:,16);    % true wind from PSD sonic - relative to earth 10/31/2012 - LB
    dir     =   x(:,41);    %
    tsnk    =   x(:,18);    % sea snake temperature, C (0.05 m depth)
% 	ta      =   x(:,19);    % air temperature, PSD   
    ta      =   MMF.Tair(ib); % use QC data
    qse     =   x(:,20);    % sea surface specific humidity, g/kg (from snake)
%   qa      =   x(:,21);    % air specific humidity, g/kg
    qa      =   MMF.Qair(ib); % use QC data
    rh      =   MMF.RH(ib);  % use QC data
    rs      =   x(:,22);    % downward solar flux, W/m^2 (PSD units)
    sig_rl  =   x(:,95);    % stdev downward solar flux, W/m^2 (PSD units)
% 	rl      =   x(:,23);    % downward IR flux, W/m^2 (PSD units)
	rl      =   MMF.IR(ib); % use QC data
	sig_rs  =   x(:,96);    % stdev downward IR flux, W/m^2 (PSD units)
    
    org     =   x(:,24);    % rainrate, mm/hr (ETL STI optical rain gauge, uncorrected)
    org(org<0) = NaN;       % remove bad rain data
    zu_etl  =   x(:,27);    % measurement heights
    zt_etl  =   x(:,28);    %
    zq_etl  =   x(:,29);    %
   
%   press   =   x(:,157);   % pressure (mbar)
    press   =   MMF.Pair(ib); %  use QC data
    
   %------ Bulk Fluxes
    Hsb     =   x(:,74);    % bulk sensible heat flux, W/m^2
    Hlb     =   x(:,75);    % bulk latent heat flux, W/m^2 (includes Webb et al. correction?)
    Taub    =   x(:,76);    % bulk wind stress along mean wind, N/m^2
      
    ustr0   =   x(:,81);
    tstr0   =   x(:,82);
    qstr0   =   x(:,83);
    Hrain   =   x(:,87);    % rain heat flux, W/m^2 as per Gosnell et al
   
   %------ Navigation & data quality indicators
    Lat     =   x(:,151);   % latitude, deg  (not SCS pcode but GPS ETL)
    Lon     =   x(:,150);   % longitude, deg (not SCS pcode but GPS ETL)
    urel    =   x(:,38);    % relative wind speed, m/s
    reldir  =   x(:,42);    % relative wind direction (from), deg clockwise rel ship's bow(ETL sonic)
    tiltx   =   x(:,44)*180/pi; % flow tilt (beta), deg
    head    =   x(:,112);   % ship heading, deg clockwise rel north, ship's gyro
    sig_h   =   x(:,113);   % std dev of heading, >5 indicates ship manauver 
    ushp_im =   x(:,97);    % ship speed, m/s (SCS GPS)
    sig_uim =   x(:,98);    % std dev of ship speed, m/s (>.2 indicates manauver during average)
    x(:,153)=   ushp_im;    % odec not functional for HiWinGS, so overwrite with gps sog
    x(:,154)=   sig_uim;
    sig_him =   x(:,102);   % std dev of ship cog (SCS)
    ushp    =   x(:,153);   % doppler log (odec)  -->replaced with SCS SOG as ODEC not available
    sig_u   =   x(:,154);   % -->replaced with SCS SOG as ODEC not available
   
%% comparison plots
%     figure;plot(jdy,U,'r.',jdy,U_scs,'.');title(' True wind speed (m/s), blue=ship, red=PSD');
%     figure;plot(jdy,dir,'r.',jdy,dir_scs,'.');title(' True wind direction (m/s), blue=ship, red=PSD');
%     figure;plot(jdy,tsnk,'r.',jdy,tsg,'.');title(' SST(degC), blue=ship tsg, red=PSD seasnake');
%     figure;plot(jdy,ta,'r.',jdy,ta_im,'.');title('Air temperature (degC), blue=ship, red=PSD');
%     figure;plot(jdy,qa,'r.',jdy,qa_im,'.');title('Air specific humidity (g/kg), blue=ship, red=PSD');
%     figure;plot(jdy,rs,'r.',jdy,rs_im,'.');title('downward solar flux, W/m^2, blue=ship, red=PSD');
%     figure;plot(jdy,rl,'r.',jdy,rl_im,'.');title('downward ir flux, W/m^2, blue=ship, red=PSD');
%     figure;plot(jdy,org,'r.',jdy,rain_im,'.');title('rainrate, mm/hr, blue=ship, red=PSD');
% 
%     ll=find(reldir>180);
%     figure;plot(reldir,tiltx,'b.',reldir(ll)-360,tiltx(ll),'b.',[-90 -90],[-20 20],'r',[90 90],[-20 20],'r');axis([-180 180 -20 20]);grid
%     title('tilt angle, degrees');
    
%% variances and covariances, turbulence data
   wu       =   x(:,116);   % hare
   wuj      =   x(:,117);   % edson
   wv       =   x(:,118);   % hare
   wvj      =   x(:,119);   % edson
   wt_son   =   x(:,121);
   uu       =   x(:,122);   % hare
   uuj      =   x(:,125);   % edson
   vvj      =   x(:,126);   % edson
   vv       =   x(:,123);   % hare
   wwj      =   x(:,127);   % edson
   ww       =   x(:,124);   % hare
   tt       =   x(:,128);
    
   sp1      =   x(:,129);   % motion stdev, u (m/s)
   sp2      =   x(:,130);   % motion stdev, v (m/s)
   sp3      =   x(:,131);   % motion stdev, w (m/s)
   
   % idiss results
%  cu       =   x(:,134);   % cu^2 method 1
%  cw       =   x(:,136);   % cw^2 method 1
   cu       =   x(:,135);   % cu^2 method 2
   cw       =   x(:,137);   % cw^2 method 2
   ct       =   x(:,138);   % ct^2 method 2
   cql      =   yco(:,243); % cq^2 Licor method 2
   usi      =   x(:,143);   % inertial dissipation scaling paramters, all method 2
   tsi      =   x(:,145);   %  
   qsi      =   x(:,147);   % 
end;

%% make modifications, calibrations, limit tests, etc
eval(['HiWinGS2013_fixit_da1'])

%% correct U for currents
% load current velocity and direction from MeansNFluxes file
MMF = load([met_path,'HiWinGS_2013_10minMeansNFluxes_r1.mat']);  
% compute U and dir relative to water using current at 14m depth
U_relative2water;
Uearth = U;          % save U & dir in earth coords
direarth = dir;
U = Ul; dir = dirl;  % U and dir now relative to water
x(:,160) = Uearth;   % U in earth coords added to da array
x(:,161) = direarth; % dir in earth coords added to da array
x(:,16)  = U;        % U is now relative to water
x(:,41)  = dir;      % dir is now relative to water

%% setup constants and other variables
tdk = 273.16;
grav = 9.72;
Rgas = 287.1;
cpa = 1004.67;
rhoair = press.*100./(Rgas*(ta+tdk).*(1+0.61*qa/1000));
ii = find(isfinite(tsnk));
Le = (2.501-.00237*median(tsnk(ii)))*1e6;

be = 0.026;
cpw = 4000;
rhow = 1022;
visw = 1e-6;
tcw = 0.6;

% seasnake deployment
tm_in = [282];      % doy deployed into water
tm_ot = [365];      % doy removed from water    
[n m] = size(jdy);
zsnk = ones(n,1)*.05;   % depth of seasnake
ztsg = ones(n,1)*5;     % depth of tsg

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

% substitute tsg when snake is out of water
tsnk_m = tsg; zts = ztsg;
[n,m] = size(tm_in);
jg=[];
for k=1:m
   jf = find(jdy'>tm_in(k) & jdy'<tm_ot(k)); % find when snake is in water
   jg = [jg jf];
end;
tsnk_m(jg) = tsnk(jg);   % SST temperature when snake is in
zts(jg) = zsnk(jg);      % depth of water temp measurement
Ts = tsnk_m;
% ii=find(isnan(tsnk));Ts(ii)=tsg(ii);
% qse=qsea(Ts);

% get map of cruise track
lat = Lat;
lon = Lon;
lond = lon;
latd = lat;
figure; map_HiWinGS2013;

% longwave clear sky radiative 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); % as f(lat)
rlcler = rlclr; % use f(lat)

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

% set instrument heights and bulk model parameters
zu = zu_etl(1);
% zt=zt_etl(1);
% zq=zq_etl(1);
zt = 15;
zq = 15;

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

% housekeep 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 = .05;  % bulk water temperature sensor depth

%% major read data loop
[nn,nx] = size(x');  % lines of 10min data  
ijx = 1;
d1 = datenum('31-Dec-2012'); % t0 for day of year
clear dt;

for ibg = 1:nx  % major read loop
	%% 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(ibg);   
    phirel = reldir(ibg);  % relative wind direction
	tsea = Ts(ibg);        % bulk sea surface temp
	t = ta(ibg);           % air temp
    P = press(ibg);        % 03/05/09 --> added pressure in the loop
	qs = qse(ibg);         % bulk sea surface humidity
	q = qa(ibg);           % air humidity
	Rs = rs(ibg);          % downward solar flux
	Rl = rl(ibg);          % downward IR flux
    rain = org(ibg);       % rain rate
    if sig_h<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
    ts_depth = zts(ibg);   % .05m, bulk water temperature sensor depth
  
    % 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
	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 = .5;
					tk_pwp = 19;
					tau_ac = 0;
					qcol_ac = 0;
					dt_wrm = 0;
					jump = 0;           % goto 16
                else                    % set warm layer constants
					rich = .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;					%idicates data start period
		end;                        % end first time thru check
		jtime=newtime;
	end;                            % end warm layer model appl check
    
    %% run COARE
    ts = tsea+dsea;
    qs = qsea(ts);
	xx = [u us ts t qs q Rs Rl rain zi  P zu zt zq latx jcool 0 0 0];  % set data for basic flux alogithm
  	y = cor30a(xx);
%   A = coare35vn(u,zu,t,zt,rh(ibg),zq,P,ts,Rs,Rl,latx,zi,rain,NaN,NaN);  % newest vectorized version
   %A=[usr tau hsb hlb hbb hsbb hlwebb tsr qsr zot zoq Cd Ch Ce  L zet dter dqer tkt Urf Trf Qrf RHrf UrfN Rnl Le rhoa UN U10 U10N Cdn_10 Chn_10 Cen_10 RF Qs Evap T10 Q10 RH10];
   %    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  36  37  38   39
    % configure to match old C3 ouput format.
%   y = [A(3) A(4) A(2) A(42) A(10) A(11) A(15) A(1) A(8) A(9)/1000 A(17) A(18) A(19) A(34) A(41) A(12) A(13) A(14) A(31)/1000 A(32)/1000 A(33)/1000 A(40) ];
    
    %% output from routine
    hsb    =    y(1);       % sensible heat flux W/m/m
    hlb    =    y(2);       % latent
    tau    =    y(3);       % stress
    zo     =    y(4);       % vel roughness
    zot    =    y(5);       % temp "
    zoq    =    y(6);       % hum  "
    L      =    y(7);       % Ob 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 @ 10 [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:
    % [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
    
    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>270  % save subset with good wind direction
        dtx(ijx,:) = zax;
        ijx = ijx+1;
    end;
end; % end major read loop

%% further processing
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 recopmputes hsc here as u* and q* have changed (for the good).
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; 

% ii = find(sgq_lic>.7 | hlc_lic<-20); hlc_lic(ii) = NaN;

hl_webb = dt(:,25).*rhoair*Le/1000.*qa;

% create best covariance latent heat
% hlc_bl = hlc_oph; % (hlc_oph+hlc_lic)/2;
hlc_bl = hlc_lic; %  Licor only since 2006.

jdollars = 'HiWinGS2013';
HiWinGS2013_da_id2;

hlib_lic(hlib_lic>1000) = NaN;

% 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

rnl = 0.97*(5.67e-8*(Ts+273.15).^4-rl);
hnet = 0.95*rs-rnl-Hsbb-Hlbb-HRF;

latmin = 52;latmax = 64;
HiWinGS2013_daily; % perform daily averages, do plots

if prtit
%   hsc(hsc<-100)=NaN;
    prt_HiWinGS2013_flux_10;      %to create an ascii file
end;

HiWinGS2013_hr_ave; % convert to one hr averages

hr_bulk_HiWinGS2013; % recompute bulk, ID, and covariance fluxes at 1-hr
med_hl_plt_HiWinGS2013;
taucx = ytauc;
taucy = ytaucy;
tauib = ytauib;
taub = ytaub;
reldir = mod(thet1rel+360,360);

med_tau_plt_HiWinGS2013;

if prtit
    prt_HiWinGS2013_flux_hr;  % to create an ascii file
end;

ii = find( isfinite(hsbb) & isfinite(hlbb));
mean([ 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));
mean([ yhlb(ii)' yhlc_bl(ii)' hlc_licbar(ii)' hl_webar(ii)]);
ii = find( isfinite(yhlc_bl) & isfinite(yhlb) & isfinite(hlc_licbar));
median([ yhlb(ii)' yhlc_bl(ii)' hlc_licbar(ii)' hl_webar(ii)]);

rhobar = yP.*100./(Rgas*(y1(:,19)+tdk).*(1+0.61*y1(:,21)/1000));
wbar_t = (1+1.61*y1(:,21)/1e3).*yhsb'./rhobar./cpa./(y(:,19)+tdk);
wbar_w = 1.61*yhlb'/Le;

reldir = x(:,42);  % relative wind direction (from), deg clockwise rel ship's bow (PSD sonic)
jj = find((reldir<90 | reldir>270) & sig_h<5 & sig_u<1.4 & sp2<1.2  & org<5 & sqrt(wwj)./ugw<2+.0015*U.^2 & sqrt(vvj)./ugu<2); %& isledist>3.5 );

figure; plot(jdy(jj),-(rhoair(jj).*wuj(jj)),'bo',jdy(jj),(rhoair(jj).*usib(jj).^2),'m+',jdy(jj),taubb(jj),'g-')
legend('Cov','ID','Bulk'); grid
xlabel('DOY (2013)'); ylabel('Stremwise stress (N/m2)');

figure; plot(jdy(jj),hsc(jj),'bo',jdy(jj),hsib(jj),'m+',jdy(jj),hsbb(jj),'g-')
legend('Cov','ID','Bulk'); grid; %ylim([-50 100])
xlabel('DOY (2013)'); ylabel('Sensible heat flux (W/m2)');

figure; plot(jdy(jj),hlc_bl(jj),'bo',jdy(jj),hlib_lic(jj),'m+',jdy(jj),hlbb(jj),'g-')
legend('Cov','ID','Bulk'); grid; %ylim([-50 300])
xlabel('DOY (2013)'); ylabel('Latent heat flux (W/m2)');
