% gfd_core_read.m
% reads core flight-level data in netCDF format using snctools nc_varget
% reads cloud-physics flight-level data in netCDF format
% reads in calibrated total water from Lyman Alpha
% saves a 1Hz data set and a smaller 32 Hz data set in Matlab format
% Note: variables have dimension: n x 1 (or n x 32, for 32 Hz)
% Written by Ian Renfrew and Dave Sproson, May 2007

% *** set up data set here; Also change filename to be saved at end ***
idata = 274; data_name = 'B274';   
if (idata == 268) 
    filea = 'C:\data\gfd\FAAM\B268\core_faam_20070221_r16_b268.nc';
    fileb = 'C:\data\gfd\FAAM\B268\core-cloud-phy_faam_20070221_r0_b268.nc';
    filec = 'C:\data\gfd\FAAM\B268\core_faam_20070221_v002_r16_b268_tw.nc';    
    dir_out = 'C:\data\gfd\FAAM\B268\';
    cd(dir_out)
    disp(['Reading ' filea])
    sd_midnight = datenum(2007, 02, 21);
elseif (idata == 274)
    filea = 'U:\NTProfile\Documents\gfd\FAAM\B274\core_faam_20070302_r15_b274.nc';
    fileb = 'U:\NTProfile\Documents\gfd\FAAM\B274\core-cloud-phy_faam_20070302_r0_b274.nc';
    filec = 'U:\NTProfile\Documents\gfd\FAAM\B274\core_faam_20070302_v002_r16_b274_tw.nc';
    filed = 'U:\NTProfile\Documents\gfd\FAAM\B274\B274_postcorr_winds.xls';
    dir_out = 'C:\data\gfd\FAAM\B274\';
    cd(dir_out)
    disp(['Reading ' filea])  
    sd_midnight = datenum(2007, 03, 02);
elseif (idata == 276)
    filea = 'C:\data\gfd\FAAM\B276\core_faam_20070305_r16_b276.nc';
    fileb = 'C:\data\gfd\FAAM\B276\core-cloud-phy_faam_20070305_r0_b276.nc';
    filec = 'C:\data\gfd\FAAM\B276\core_faam_20070305_v002_r16_b276_tw.nc';
    dir_out = 'C:\data\gfd\FAAM\B276\';
    cd(dir_out)
    disp(['Reading ' filea])  
    sd_midnight = datenum(2007, 03, 05);
elseif (idata == 278)
    filea = 'C:\data\gfd\FAAM\B278\core_faam_20070309_r16_b278.nc';
    fileb = 'C:\data\gfd\FAAM\B278\core-cloud-phy_faam_20070309_r0_b278.nc';
    filec = 'C:\data\gfd\FAAM\B278\core_faam_20070309_v002_r16_b278_tw.nc';    
    dir_out = 'C:\data\gfd\FAAM\B278\';
    cd(dir_out)
    disp(['Reading ' filea])  
    sd_midnight = datenum(2007, 03, 09);
end;

% Read in time (secs after midnight) and define serial date
% Note: added interplotation over NAN with time variable
time = nc_varget(filea, 'PARA0515');   
time = naninterp(time);
sd = sd_midnight + time./86400;
[yyyy, mm, dd, hh, mi, ss] = datevec(sd);

% Read in the lat and lon arrays and interpolate over bad data
lat = nc_varget(filea, 'PARA0730');
lon = nc_varget(filea, 'PARA0731');
lat(find(nc_varget(filea, 'PARA0730FLAG')~=0)) = NaN;
lon(find(nc_varget(filea, 'PARA0731FLAG')~=0)) = NaN;
lat = naninterp(lat);
lon = naninterp(lon);

%-------------------------------------------------------------------
% core flight-level variables, if FLAG is not 0 then set to NaN
%-------------------------------------------------------------------
air_temp_32hz = nc_varget(filea, 'PARA0525');
air_temp_32hz(find(nc_varget(filea, 'PARA0525FLAG')~=0)) = NaN;
for i = 1:size(air_temp_32hz,1) air_temp(i) = nanmean(air_temp_32hz(i,:)); end

deiced_air_temp_32hz = nc_varget(filea, 'PARA0520');
deiced_air_temp_32hz(find(nc_varget(filea, 'PARA0520FLAG')~=0)) = NaN;
for i = 1:size(deiced_air_temp_32hz,1) deiced_air_temp(i) = nanmean(deiced_air_temp_32hz(i,:)); end

dewpoint_ge_4hz = nc_varget(filea, 'PARA0529');
dewpoint_ge_4hz(find(nc_varget(filea, 'PARA0529FLAG')~=0)) = NaN;
for i = 1:size(dewpoint_ge_4hz,1) dewpoint_ge(i) = nanmean(dewpoint_ge_4hz(i,:)); end

lwc_jw_4hz = nc_varget(filea, 'PARA0535');
lwc_jw_4hz(find(nc_varget(filea, 'PARA0535FLAG')~=0)) = NaN;
for i = 1:size(lwc_jw_4hz,1) lwc_jw(i) = nanmean(lwc_jw_4hz(i,:)); end

heading_inu_32hz = nc_varget(filea, 'PARA0562');
heading_inu_32hz(find(nc_varget(filea, 'PARA0562FLAG')~=0)) = NaN;
for i = 1:size(heading_inu_32hz,1) heading_inu(i) = nanmean(heading_inu_32hz(i,:)); end

ozone = nc_varget(filea, 'PARA0574');
ozone(find(nc_varget(filea, 'PARA0574FLAG')~=0)) = NaN;

radar_alt_2hz = nc_varget(filea, 'PARA0575');
radar_alt_2hz(find(nc_varget(filea, 'PARA0575FLAG')~=0)) = NaN;
for i = 1:size(radar_alt_2hz,1) radar_alt(i) = nanmean(radar_alt_2hz(i,:)); end

static_pressure_32hz = nc_varget(filea, 'PARA0576');
static_pressure_32hz(find(nc_varget(filea, 'PARA0576FLAG')~=0)) = NaN;
for i = 1:size(static_pressure_32hz,1) static_pressure(i) = nanmean(static_pressure_32hz(i,:)); end

static_pressure_s9_32hz = nc_varget(filea, 'PARA0778');
static_pressure_s9_32hz(find(nc_varget(filea, 'PARA0778FLAG')~=0)) = NaN;
for i = 1:size(static_pressure_s9_32hz,1) static_pressure_s9(i) = nanmean(static_pressure_32hz(i,:)); end

pressure_alt_32hz = nc_varget(filea, 'PARA0578');
pressure_alt_32hz(find(nc_varget(filea, 'PARA0578FLAG')~=0)) = NaN;
for i = 1:size(pressure_alt_32hz,1) pressure_alt(i) = nanmean(pressure_alt_32hz(i,:)); end

gps_alt = nc_varget(filea, 'PARA0582');
gps_alt(find(nc_varget(filea, 'PARA0582FLAG')~=0)) = NaN;

lwc_nev_8hz = nc_varget(filea, 'PARA0602');
lwc_nev_8hz(find(nc_varget(filea, 'PARA0602FLAG')~=0)) = NaN;
for i = 1:size(lwc_nev_8hz,1) lwc_nev(i) = nanmean(lwc_nev_8hz(i,:)); end

twc_nev_8hz = nc_varget(filea, 'PARA0605');
twc_nev_8hz(find(nc_varget(filea, 'PARA0605FLAG')~=0)) = NaN;
for i = 1:size(twc_nev_8hz,1) twc_nev(i) = nanmean(twc_nev_8hz(i,:)); end

twc_lyman_64hz = nc_varget(filea, 'PARA0664');
twc_lyman_64hz(find(nc_varget(filea, 'PARA0664FLAG')~=0)) = NaN;
for i = 1:size(twc_lyman_64hz,1) twc_lyman(i) = nanmean(twc_lyman_64hz(i,:)); end

v_32hz = nc_varget(filea, 'PARA0714');
v_32hz(find(nc_varget(filea, 'PARA0714FLAG')~=0)) = NaN;
for i = 1:size(v_32hz,1) v(i) = nanmean(v_32hz(i,:))'; end

u_32hz = nc_varget(filea, 'PARA0715');
u_32hz(find(nc_varget(filea, 'PARA0715FLAG')~=0)) = NaN;
for i = 1:size(u_32hz,1) u(i) = nanmean(u_32hz(i,:))'; end

w_32hz = nc_varget(filea, 'PARA0716');
w_32hz(find(nc_varget(filea, 'PARA0716FLAG')~=0)) = NaN;
for i = 1:size(w_32hz,1) w(i) = nanmean(w_32hz(i,:))'; end

if (nc_isvar(filea, 'PARA0723'))                % check variable exists
    sfc_temp_4hz = nc_varget(filea, 'PARA0723');
    sfc_temp_4hz(find(nc_varget(filea, 'PARA0723FLAG')~=0)) = NaN;
    for i = 1:size(sfc_temp_4hz,1) sfc_temp(i) = nanmean(sfc_temp_4hz(i,:)); end;
else
    sfc_temp = [];
    sfc_temp_4hz = [];
end

total_scat_blue_neph = nc_varget(filea, 'PARA0762');
total_scat_blue_neph(find(nc_varget(filea, 'PARA0762FLAG')~=0)) = NaN;

total_scat_green_neph = nc_varget(filea, 'PARA0763');
total_scat_green_neph(find(nc_varget(filea, 'PARA0763FLAG')~=0)) = NaN;

total_scat_red_neph = nc_varget(filea, 'PARA0764');
total_scat_red_neph(find(nc_varget(filea, 'PARA0764FLAG')~=0)) = NaN;

swd_clear = nc_varget(filea, 'PARA1019');
swd_clear(find(nc_varget(filea, 'PARA1019FLAG')~=0)) = NaN;

lwd = nc_varget(filea, 'PARA1021');
lwd(find(nc_varget(filea, 'PARA1021FLAG')~=3)) = NaN;
disp('Note: lwu flag of 3 is used')

swu_clear = nc_varget(filea, 'PARA1022');
swu_clear(find(nc_varget(filea, 'PARA1022FLAG')~=0)) = NaN;

lwu = nc_varget(filea, 'PARA1024');
lwu(find(nc_varget(filea, 'PARA1024FLAG')==3)) = NaN;
disp('Note: lwu flag of 3 is used')

%-----------------------------------------
% Cloud Physics variables
%-----------------------------------------
disp(['Reading ' fileb])

part_conc_2dc = nc_varget(fileb, 'PARA1301');
part_conc_2dc(find(nc_varget(fileb, 'PARA1301FLAG')~=0)) = NaN;

lwc_2dc = nc_varget(fileb, 'PARA1302');
lwc_2dc(find(nc_varget(fileb, 'PARA1302FLAG')~=0)) = NaN;

part_conc_2dp = nc_varget(fileb, 'PARA1311');
part_conc_2dp(find(nc_varget(fileb, 'PARA1311FLAG')~=0)) = NaN;

lwc_2dp = nc_varget(fileb, 'PARA1312');
lwc_2dp(find(nc_varget(fileb, 'PARA1312FLAG')~=0)) = NaN;

part_conc_pcasp = nc_varget(fileb, 'PARA1550');
part_conc_pcasp(find(nc_varget(fileb, 'PARA1550FLAG')~=0)) = NaN;

mass_content_pcasp = nc_varget(fileb, 'PARA1551');
mass_content_pcasp(find(nc_varget(fileb, 'PARA1551FLAG')~=0)) = NaN;

%------------------------------------
% Lyman alpha total water variable
%------------------------------------
disp(['Reading ' filec])
q_tw_64hz = nc_varget(filec, 'TWC_EVAP');
for i = 1:size(q_tw_64hz,1) q_tw(i) = nanmean(q_tw_64hz(i,:)); end

%----------------------------------------------------------------------
% INU-derived winds (use these for parts of B274, 276, 277 & 278) 
%--------------------------------------------------------------------
disp(['Reading ' filed])


%-----------------------------------------------------------------
% Save 1hz data to a matlab file
%----------------------------------------------------------------
disp('saving flightlevel 1Hz data')
save  B274_flightlevel_1hz_r16.mat ...
	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 static_pressure_s9...
	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 ...
    q_tw

 disp('saving flightlevel 32hz data')
 save B274_flightlevel_32hz_r16.mat ...
 	time sd hh mi ss lat lon heading_inu ...
     air_temp_32hz dewpoint_ge_4hz twc_lyman_64hz static_pressure_32hz static_pressure_s9_32hz sfc_temp_4hz ...
     radar_alt_2hz pressure_alt_32hz gps_alt ...
     u_32hz v_32hz w_32hz q_tw_64hz