% InterpSondes.m - reads sondes for each flight and interpolates to fixed
% vertical grid, filling NaNs except at start and end.

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

path_data = '/Users/byronb/Raw_Data_Files.noindex/ATOMIC_2020/drops/raw';
cd(path_data);

z_ref = 10:10:12000; z_ref = z_ref'; % fixed altitude vector 10-12000 m

for ii=1:11 % iterate through 11 flights
    % read raw profile data
    fname = ['*_RF',sprintf('%02.0f',ii),'.mat'];
    file = dir(fname);
    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
    
    % wind speed
    U_raw = snd_raw.Uhh;
    U_interp = NaN(1200,N);
    for kk = 1:N
        nn = find(isfinite(z_raw(:,kk)) & isfinite(U_raw(:,kk)));
        if isempty(nn); continue; end
        U_interp(:,kk) = interp1(z_raw(nn,kk),U_raw(nn,kk),z_ref);
    end
    
    % wind direction 'from'
    udir = snd_raw.udirhh;
    dir_raw = udir*pi/180;
    dir_raw = unwrap(dir_raw);
    dir_interp = NaN(1200,N);
    for kk = 1:N
        nn = find(isfinite(z_raw(:,kk)) & isfinite(dir_raw(:,kk)));
        if isempty(nn); continue; end
        dir_interp(:,kk) = interp1(z_raw(nn,kk),dir_raw(nn,kk),z_ref);
        dir_interp(:,kk) = dir_interp(:,kk)*180/pi;
        dir_interp(:,kk) = mod(dir_interp(:,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_interp .* cos(dir_interp*pi/180);
    Ue = U_interp .* sin(dir_interp*pi/180);
    
    % air temp
    ta_raw = snd_raw.tahh;
    ta_interp = NaN(1200,N);
    for kk = 1:N
        nn = find(isfinite(z_raw(:,kk)) & isfinite(ta_raw(:,kk)));
        if isempty(nn); continue; end
        ta_interp(:,kk) = interp1(z_raw(nn,kk),ta_raw(nn,kk),z_ref);
    end
    
    % rh
    rh_raw = snd_raw.rhahh;
    rh_raw(isnan(rh_raw) & z_raw<12000) = 1;
    rh_interp = NaN(1200,N);
    for kk = 1:N
        nn = find(isfinite(z_raw(:,kk)) & isfinite(rh_raw(:,kk)));
        if isempty(nn); continue; end
        rh_interp(:,kk) = interp1(z_raw(nn,kk),rh_raw(nn,kk),z_ref);
    end
    
    % qair
    qa_raw = snd_raw.qahh;
    qa_interp = NaN(1200,N);
    for kk = 1:N
        nn = find(isfinite(z_raw(:,kk)) & isfinite(qa_raw(:,kk)));
        if isempty(nn); continue; end
        qa_interp(:,kk) = interp1(z_raw(nn,kk),qa_raw(nn,kk),z_ref);
    end
    
    % P
    P_raw = snd_raw.Phh;
    P_interp = NaN(1200,N);
    for kk = 1:N
        nn = find(isfinite(z_raw(:,kk)) & isfinite(P_raw(:,kk)));
        if isempty(nn); continue; end
        P_interp(:,kk) = interp1(z_raw(nn,kk),P_raw(nn,kk),z_ref);
    end
    
    % save interpolated profiles to new structure
    sname = [file.name(1:24),'_zinterp.mat'];
    snd_z = struct('jd',snd_raw.jdxh,'lon',snd_raw.lonhh,'lat',snd_raw.lathh,...
        'z',z_ref,'U',U_interp,'udir',dir_interp,'Un',Un,'Ue',Ue,...
        'ta',ta_interp,'rh',rh_interp,'P',P_interp,'qa',qa_interp);
    save(sname,'snd_z');
end



