% InterpSondes_P_modelGrid.m - reads sondes for each flight and bin-averages to fixed
% model pressure grid.

% float lv_ISBL6(lv_ISBL6) ;
% 		lv_ISBL6:units = "Pa" ;
% 		lv_ISBL6:long_name = "Isobaric surface"
% 
% lv_ISBL6 = 10000, 15000, 20000, 25000, 30000, 35000, 40000, 45000, 50000, 
%     55000, 60000, 65000, 70000, 75000, 80000, 85000, 90000, 92500, 95000, 
%     97500, 100000

% matlab script path
path_prog = '/Users/byronb/Documents/MATLAB/2016_ENRR';
cd(path_prog);                  % cd to matlab directory
restoredefaultpath              % reset matlab path
addpath(genpath(path_prog));	% add scripts to path
rehash;

path_data = '/Users/byronb/Raw_Data_Files.noindex/ENRR_2016/Fairall_processing/data2/';
cd(path_data);

% model grid bin edges
p_bin = [75;125;175;225;275;325;375;425;475;525;575;625;675;725;775;825;875;912.5;937.5;962.5;987.5;1015];
p_ref = [1000;975;950;925;900;850;800;750;700;650;600;550;500;450;400;350;300;250;200;150;100];

for ii=1:22 % iterate through 22 flights
    % read raw profile data
    fname = ['*_RF',sprintf('%02.0f',ii),'.mat'];
    file = dir(fname);
    disp(['Processing ',file.name]);
    snd_raw = load(fullfile(path_data,file.name), 'jdxh', 'zhh','Uhh','rhahh','tahh','lathh','lonhh','Phh','qahh','udirhh');
    N = length(snd_raw.jdxh); % number of drops on this flight
    z_raw = snd_raw.zhh;      % altitude vectors, bottom up
    P_filt = sgolayfilt(snd_raw.Phh,1,25);
    
    % wind speed
    U_bin = NaN(21,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(snd_raw.Uhh(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        U_bin(:,kk) = flipud(temp(:,2));
    end
    
    % wind direction 'from'
    dir_raw = snd_raw.udirhh;
    dir_raw = dir_raw*pi/180;
    dir_raw = unwrap(dir_raw);
    dir_bin = NaN(21,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(dir_raw(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        dir_bin(:,kk) = flipud(temp(:,2));
        dir_bin(:,kk) = dir_bin(:,kk)*180/pi;
        dir_bin(:,kk) = mod(dir_bin(:,kk)+8*360, 360);
    end
    
    % compute Un/Ue components, where Un is positive S to N and Ue is
    % positive for E to W (right handed convention).
    Un = U_bin .* cos(dir_bin*pi/180);
    Ue = U_bin .* sin(dir_bin*pi/180);
    
    % air temp
    ta_bin = NaN(21,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(snd_raw.tahh(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        ta_bin(:,kk) = flipud(temp(:,2));
    end
    
    % rh
    rh_raw = snd_raw.rhahh;
    rh_raw(isnan(rh_raw) & z_raw<12000) = 1;
    rh_bin = NaN(21,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(rh_raw(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        rh_bin(:,kk) = flipud(temp(:,2));
    end
    
    % qair
    qa_bin = NaN(21,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(snd_raw.qahh(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        qa_bin(:,kk) = flipud(temp(:,2));
    end
    
    % z
    z_bin = NaN(21,N);
    for kk = 1:N
        p1 = flipud(P_filt(:,kk));
        x1 = flipud(z_raw(:,kk));
        temp = interval_avg(p1,x1,p_bin);
        z_bin(:,kk) = flipud(temp(:,2));
    end
    
    % save interpolated profiles to new structure
    sname = [file.name(1:28),'_pinterp_modelGrid.mat'];
    snd_pinterp_modelGrid = struct('jd',snd_raw.jdxh,'lon',snd_raw.lonhh,'lat',snd_raw.lathh,...
        'P',p_ref,'z',z_bin,'U',U_bin,'udir',dir_bin,'Un',Un,'Ue',Ue,'ta',ta_bin,...
        'rh',rh_bin,'qa',qa_bin);
    save(sname,'snd_pinterp_modelGrid');
end



