% interp_sondes_halo_nc_P.m - reads halo dropsondes from folder
% and interps to pressure grid with 20mb bin width and 41 bins
% from 1000 to 200 mb (bin centers).

%% initialize
% matlab script path
path_prog = '/Users/byronb/Documents/MATLAB/2020_ATOMIC';
cd(path_prog);                  % cd to matlab directory
addpath(genpath(path_prog));	% add scripts to path
rehash;

% set this path and run
path_data = '/Users/byronb/Raw_Data_Files.noindex/ATOMIC_2020/halo_drops/raw_nc/2020_02_13';
cd(path_data);

% for saving mat files
save_path = '/Users/byronb/Raw_Data_Files.noindex/ATOMIC_2020/halo_drops/mat';
mkdir(save_path);

% 39 level fixed pressure bin edges 250-1010 mb
p_bin = 250:20:1010; p_bin = p_bin';

%% load all variables
files = dir('*.nc'); NN = length(files);
jdxh = NaN(NN,1); lathh = jdxh; lonhh = jdxh;
U_bin = NaN(38,NN);
dir_bin = NaN(38,NN);
ta_bin = NaN(38,NN);
rh_bin = NaN(38,NN);
qa_bin = NaN(38,NN);
z_bin = NaN(38,NN);

for ii = 1:NN
    dfl = files(ii).name;
    snd = load_nc(dfl);
    yyyy = str2double(dfl(2:5)); mm = str2double(dfl(6:7)); dd = str2double(dfl(8:9));
    hh = str2double(dfl(11:12)); mn = str2double(dfl(13:14)); ss = str2double(dfl(15:16));
    jd0 = md2yd_vn(yyyy,mm,dd) + hh/24 + mn/1440 + ss/86400;
    jdxh(ii) = jd0;
    lathh(ii) = snd.reference_lat;
    lonhh(ii) = snd.reference_lon;
    zhh = snd.alt;
    udirhh = snd.wdir;
    Uhh = sqrt(snd.u_wind.^2 + snd.v_wind.^2);
    rhahh = snd.rh;
    tahh = snd.tdry;
    Phh = snd.pres;
    td = snd.dp;
    qahh = td2q(td,Phh);
    
    % wind speed
    p1 = flipud(Phh);
    x1 = flipud(Uhh);
    temp = interval_avg(p1,x1,p_bin);
    U_bin(:,ii) = flipud(temp(:,2));
    if ii==1
        p_ref = flipud(temp(:,1)) + 10;  % reference pressure bin centers
    end
    
    % wind direction 'from'
    dir_raw = udirhh*pi/180;
    dir_raw = unwrap(dir_raw);
    x1 = flipud(dir_raw);
    temp = interval_avg(p1,x1,p_bin);
    dir_bin(:,ii) = flipud(temp(:,2));
    dir_bin(:,ii) = dir_bin(:,ii)*180/pi;
    dir_bin(:,ii) = mod(dir_bin(:,ii)+8*360, 360);
    
    % air temp
    x1 = flipud(tahh);
    temp = interval_avg(p1,x1,p_bin);
    ta_bin(:,ii) = flipud(temp(:,2));

    % rh
    x1 = flipud(rhahh);
    temp = interval_avg(p1,x1,p_bin);
    rh_bin(:,ii) = flipud(temp(:,2));
  
    % qair
    x1 = flipud(qahh);
    temp = interval_avg(p1,x1,p_bin);
    qa_bin(:,ii) = flipud(temp(:,2));

    % z
    x1 = flipud(zhh);
    temp = interval_avg(p1,x1,p_bin);
    z_bin(:,ii) = flipud(temp(:,2));
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);

%% save to structure

sname = fullfile(save_path,['halo_',files(1).name(2:9),'_all_pinterp.mat']);
snd_p = struct('jd',jdxh','lon',lonhh','lat',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_p');



