% gfd_flightlevel_diagnostics
% calculates number of diagnostics for matlab-formatted core data
% use after gfd_core_read or load relevant .mat file
%
% matlab-formatted flightlevel variables:
% 	time sd hh mi ss lat lon heading_inu ...
% 	air_temp deiced_air_temp dewpoint_ge ...
% 	ozone ...
% 	radar_alt static_pressure pressure_alt gps_alt ...
% 	lwc_jw lwc_nev twc_nev twc_lyman ...
% 	u v w sfc_temp ...
% 	total_scat_blue_neph total_scat_green_neph total_scat_red_neph ...
% 	swd_clear lwd swu_clear lwu ...
% 	part_conc_2dc lwc_2dc part_conc_2dp lwc_2dp part_conc_pcasp mass_content_pcasp 
%     
disp('Calculating flightlevel diagnostics')

% wind speed, wind direction
[ws, wd] = uv_to_wswd(u, v);    

% calculate air density in units kg/m^3
rho = (static_pressure*100)./(287*air_temp);

% thermodynamic variables: non-deiced (take these as default)
 [air_theta, air_thetae, q_ge, rh_ge, rhi_ge, qsat, qsati] = thermo2_td(air_temp-273.15, static_pressure, dewpoint_ge-273.15);

% Calculate SST from upward LW, emissivity=0.97 (Curry and Webster p251) 
stefan_boltzmann = 5.67e-8; 
sfc_temp_lwu = NaN*ones(size(lwu));
ee = find(~isnan(lwu));
sfc_temp_lwu(ee) = ((0.97*lwu(ee)/stefan_boltzmann).^0.25);
clear ee;

% calculate surface bulk turbulent fluxes; only valid for low-level legs 
% radar alt indicated ~100 feet so about ~30 m
sst = 278*ones(size(ws));
z_flux_calc = 30*ones(1,4);
[shf_bulk, lhf_bulk, tau_bulk] = sfcflux_hf_robust(ws, air_temp-273.14, rh_ge, sst-273.14, static_pressure, z_flux_calc);

% save diagnostics to .mat file
idata = 268; data_name = 'B268';   
if (idata == 268) 
    dir_out = 'C:\data\gfd\FAAM\B268\';
    cd(dir_out)
    disp('saving B268_flightlevel_diagnostics_1hz_v1')
    save B268_flightlevel_diagnostics_1hz_v1.mat ...
    ws wd ...
    air_theta air_thetae q_ge rh_ge rhi_ge qsat qsati ...
    sfc_temp_lwu shf_bulk lhf_bulk tau_bulk
end;