% ENRR_sonde_profiles.m

clear all
close all
printit = false;

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

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

tdk = 273.16;
Rgas = 287.1;
cpa = 1004.67; % air specific heat

% set desired flight numbers and loop
for ii=1 %ii=[1,2,3,4,5,7,8,10,14,15,16] % all wrap, skinny wrap and coathanger flights
    
    path_plots = ['/Users/byronb/Raw_Data_Files.noindex/ENRR_2016/Fairall_processing/plots/RF',sprintf('%02.0f',ii)];
    mkdir(path_plots);
    
    % read one flight, drops interpolated to fixed altitude vector
    fname = ['*_RF',sprintf('%02.0f',ii),'_zinterp.mat'];
    file = dir(fname);
    % structure with fields: jd, lon, lat, z, U, udir, Un, Uw, ta, rh, P, qa
    load(fullfile(path_data,file.name));
    N = length(snd_zinterp.jd); % number of drops on this flight
    
    % set northern-most limit for convective analysis based of flight number
    switch ii
        case {1,2,8,10,15}
            latlim = 6;
        case 16
            latlim = 9;
        case {3,4,5,7}
            latlim = 11;
        case 14
            latlim = 18;
    end
    
    % read variables for all drops.  Profiles interp to fixed alt grid.
    jd    = snd_zinterp.jd;     % 1xN array, time of drop
    lon   = snd_zinterp.lon;    % 1xN array, lon at start of drop
    lat   = snd_zinterp.lat;    % 1xN array, lat at start of drop
    z_ref = snd_zinterp.z;      % 1200x1 altitude vactor for all profiles, m
    U     = snd_zinterp.U;      % 1200xN profile arrays...
    udir  = snd_zinterp.udir;
    Un    = snd_zinterp.Un;
    Ue    = snd_zinterp.Ue;
    ta    = snd_zinterp.ta;
    rh    = snd_zinterp.rh;
    P     = snd_zinterp.P;
    qa    = snd_zinterp.qa;
    thasx = (ta + tdk).*(1010./P).^0.287 - tdk; % air theta
    
    stdtdol = datestr(jd(1) + datenum(2015,12,31)); % time of first drop
    
    disp(['Number of profiles= ' num2str(N)])
    
    latax = ones(1200,N).*lat;      % 1200xN lat array, using constant lat/lon for each drop
    lonax = ones(1200,N).*lon;
    zax = ones(1200,N).*z_ref/1e3;  % 120xN alt array in km, same z_ref for each drop
    
    [xi,yi] = meshgrid(-2.5:.5:20,0:.5:12);  % latitude : altitude grid
    jj = find(isfinite(rh) & isfinite(latax) & isfinite(zax));
    ZI = griddata(latax(jj),zax(jj),rh(jj),xi,yi);
    v = 0:5:105;
    figure;
    subplot(2,1,1); hold on;
    [c,h] = contourf(xi,yi,ZI,v); clabel(c,h), colorbar;
    title(['ENRR2016  NOAA G-4 RH (%) JD= ',sprintf('%03.2f',jd(1)),' - ',sprintf('%03.2f',jd(N))]);
    xlabel('Latitude (deg)');ylabel('Altitude (km)');
    hold off;
    
    [xi,yi] = meshgrid(-2.5:.5:20,0:.1:3);  % latitude : altitude grid
    jj = find(isfinite(thasx) & isfinite(latax) & isfinite(zax));
    ZI = griddata(latax(jj),zax(jj),thasx(jj),xi,yi);
    v = 15:2:50;
    subplot(2,1,2); hold on;
    [c,h] = contourf(xi,yi,ZI,v); clabel(c,h), colorbar;
    title(['ENRR2016  NOAA G-4 Theta (C) JD= ',sprintf('%03.2f',jd(1)),' - ',sprintf('%03.2f',jd(N))]);
    xlabel('Latitude (deg)');ylabel('Altitude (km)');
    hold off;
    if printit
        pname = fullfile(path_plots,['RH_Theata_LatSection_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    [xi,yi]=meshgrid(-2.5:.5:20,0:.5:13);
    jj = find(isfinite(Un) & isfinite(latax) & isfinite(zax));
    ZI = griddata(latax(jj),zax(jj),Un(jj),xi,yi);
    v = -40:2:40;
    figure;
    subplot(2,1,1); hold on;
    [c,h] = contourf(xi,yi,ZI,v); clabel(c,h), colorbar;
    title(['ENRR2016  NOAA G-4 Un (m/s) JD= ',sprintf('%03.2f',jd(1)),' - ',sprintf('%03.2f',jd(N))]);
    xlabel('Latitude (deg)');ylabel('Altitude (km)');
    hold off;

    [xi,yi]=meshgrid(-2.5:.5:20,0:.5:13);
    jj = find(isfinite(Ue) & isfinite(latax) & isfinite(zax));
    ZI = griddata(latax(jj),zax(jj),Ue(jj),xi,yi);
    v = -40:2:40;
    subplot(2,1,2); hold on;
    [c,h] = contourf(xi,yi,ZI,v); clabel(c,h), colorbar;
    title(['ENRR2016  NOAA G-4 Ue (m/s) JD= ',sprintf('%03.2f',jd(1)),' - ',sprintf('%03.2f',jd(N))]);
    xlabel('Latitude (deg)');ylabel('Altitude (km)');
    hold off;
    if printit
        pname = fullfile(path_plots,['Un_Ue_LatSection_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    % smooth profiles with savitsky-golay filter
    % unaxs = sgolayfilt(Un,1,25);
    % ueaxs = sgolayfilt(Ue,1,25);
    unaxs = Un;
    ueaxs = Ue;
    paxs = sgolayfilt(P,1,25);
    qaxs = sgolayfilt(qa,1,25);
    rhoax = paxs*100./(Rgas*(ta + tdk).*(1 + 0.61*qaxs/1e3));  % air density
    
    % bin mean profiles by latitude from -2° to 17°
    [~, Yp,  ~, ~, ~] = binave2(lat,paxs,2,-2,17);
    [~, Yq,  ~, ~, ~] = binave2(lat,qaxs,2,-2,17);
    [~, Yun, ~, ~, ~] = binave2(lat,ueaxs,2,-2,17);
    [~, Yue, ~, ~, ~] = binave2(lat,unaxs,2,-2,17);
    [~, Yh,  ~, ~, ~] = binave2(lat,rh,2,-2,17);
    [lat_bins, Yth, ~, ~, ~] = binave2(lat,thasx,2,-2,17);
    
    % latitude gradient plots
    pp = find(z_ref==20);    % pick 20m surface level
    figure;
    subplot(4,1,1); plot(lat_bins,Yp(:,pp)','o-','linewidth',2); grid;
    ylabel('P (mb)'); title(['Mean Surface Gradients: ENRR  NOAA G-4 ',stdtdol(1:12) ] );
    subplot(4,1,2); plot(lat_bins,Yq(:,pp)','o-','linewidth',2); grid;
    ylabel('q (g/kg)');
    subplot(4,1,3); plot(lat_bins,Yun(:,pp)','o-',lat_bins,Yue(:,pp)','o-','linewidth',2); grid;
    ylabel('U (m/s)'); legend('Un','Ue','location','best');
    subplot(4,1,4); plot(lat_bins,Yth(:,pp)','o-','linewidth',2); grid;
    ylabel('Theta (C)'); xlabel('Latitude');
    if printit
        pname = fullfile(path_plots,['Lat_Surface_Gradients_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    
    an = NaN(3,1200);
    phipx = paxs*0;
    ae=an; aq=an; ap=an; ath=an;
    sigpax = NaN(1200,1);
    paxm = sigpax;  ns = sigpax;  rhoaxm = sigpax; qaxm = sigpax;
    unaxm = sigpax; ueaxm = sigpax; thaxm = sigpax;
    
    % dX/dLon and dX/dLat at all altitudes for wind, q, P and theta
    for i=1:1200
        kk = isfinite(unaxs(i,:)) & latax(i,:)<latlim;  % pick finite values at altitude i where lat < latlim
        aa = [ones(size(latax(i,kk)))',lonax(i,kk)',latax(i,kk)'];
        cc = aa\unaxs(i,kk)';
        an(:,i) = cc;
        
        kk = isfinite(ueaxs(i,:)) & latax(i,:)<latlim;
        aa = [ones(size(latax(i,kk)))',lonax(i,kk)',latax(i,kk)'];
        cc = aa\ueaxs(i,kk)';
        ae(:,i) = cc;
        
        kk = isfinite(qaxs(i,:)) & latax(i,:)<latlim;
        aa = [ones(size(latax(i,kk)))',lonax(i,kk)',latax(i,kk)'];
        cc = aa\qaxs(i,kk)';
        aq(:,i) = cc;
        
        kk = isfinite(thasx(i,:)) & latax(i,:)<latlim;
        aa = [ones(size(latax(i,kk)))',lonax(i,kk)',latax(i,kk)'];
        cc = aa\thasx(i,kk)';
        ath(:,i) = cc;
        
        kk = find(isfinite(paxs(i,:)) & latax(i,:)<latlim);
        sigpax(i) = nanstd(paxs(i,kk)); paxm(i) = nanmedian(paxs(i,kk));
        mm = find(abs(paxs(i,kk) - paxm(i)) < 4*sigpax(i)); ns(i) = length(mm);
        aa = [ones(size(latax(i,kk(mm))))',lonax(i,kk(mm))',latax(i,kk(mm))'];
        cc = aa\paxs(i,kk(mm))';
        ap(:,i) = cc;
        
        % compute mean profiles
        kk = isfinite(qaxs(i,:)) & latax(i,:)<latlim;
        rhoaxm(i) = nanmean(rhoax(i,kk));
        paxm(i)   = nanmean(paxs(i,kk));
        qaxm(i)   = nanmean(qaxs(i,kk));
        unaxm(i)  = nanmean(unaxs(i,kk));
        ueaxm(i)  = nanmean(ueaxs(i,kk));
        thaxm(i)  = nanmean(thasx(i,kk));
        
%         if i<1200
%             kk = isfinite(paxs(i,:))& latax(i,:)<latlim;
%             phipx(i+1,kk) = phipx(i,kk) - 287/9.83*(ta(i,kk)+tdk).*(paxs(i+1,kk) - paxs(i,kk))./paxs(i,kk);
%         end
    end
    
%     [km,~] = size(Yp);
%     Ypm = NaN(size(Yp)); Yprhm = Ypm;
%     for i=1:km
%         Ypm(i,:) = Yp(i,:) - nanmean(Yp(i,:));  % demeaned pressure profiles in lat bins
%         Yprhm(i,:) = Ypm(i,:)./rhoaxm';         % divided by air density
%     end
    
    % dv = dUn/dlat + dUe/dlon
    dv = -(an(3,:) + ae(2,:))/1e5; % horizontal divergence - actually 1.1058e5 m per degree
    cl = -(an(2,:) - ae(3,:))/1e5; % vertical vorticity
    
    % compute vertical velocity by integrating horizontal divergence
    wi = zeros(1,1200);
    for i=2:1200
        % integral of divergence = vert. vel.
        wi(i) = wi(i-1) - dv(i)*10; % Vertical velocity, 10m spacing between levels.
    end
    
    % smoothed mean q and theta profiles over all valid drops
    qaxms = sgolayfilt(qaxm,1,55);
    thaxms = sgolayfilt(thaxm,1,25);
    
    % vertical gradient of specific humidty
    dqdz = diff(qaxms)/10; dqdz = [dqdz;0];
    
    % vertical gradient of theta
    dthdz = diff(thaxms)/10; dthdz = [dthdz;0];
    
    % vertical gradient of water vapor density
    aa = 0.1; dqdz = filter(aa,[1 aa-1],dqdz);
    dQdz = rhoaxm.*dqdz;
    
    % smooth mean wind component profiles
    unaxms = sgolayfilt(unaxm,1,25);
    ueaxms = sgolayfilt(ueaxm,1,25);
    
    % vertical derivative
    dundz = diff(unaxms)/10; dundz = [dundz;0];
    duedz = diff(ueaxms)/10; duedz = [duedz;0];
    
    % average of f*v profiles
    kk = latax(1,:)<latlim;
    fun = nanmean((unaxs(:,kk).*sin(pi/180*latax(:,kk))),2)*1.45e-4;
    fue = nanmean((ueaxs(:,kk).*sin(pi/180*latax(:,kk))),2)*1.45e-4;
    
    % plots
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(-(an(3,:)' + ae(2,:)')/1e5, paxm,'-','linewidth',2); grid;
    hold off; set(gca,'yDir','reverse'); ylim([200 1010]);
    ylabel('P (mb)'); xlabel('Divergence (s^{-1})'); title(['ENRR  NOAA G-4 ',stdtdol(1:12) ] );
    if printit
        pname = fullfile(path_plots,['Divergence_Profile_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot((ae(3,:)'-an(2,:)')/1e5, paxm,'-','linewidth',2); grid;
    hold off; set(gca,'yDir','reverse'); ylim([200 1010]);
    ylabel('P (mb)'); xlabel('Curl (s^{-1})'); title(['ENRR  NOAA G-4 ',stdtdol(1:12) ] );
    if printit
        pname = fullfile(path_plots,['Curl_Profile_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(aq(3,:)'/1e5,paxm,'-',aq(2,:)'/1e5,paxm,'-','linewidth',2); grid;
    hold off; set(gca,'yDir','reverse'); ylim([200 1010]);
    ylabel('P (mb)'); xlabel('Q gradients g/kg/m'); legend(h1,'North grad','East grad','location','best');
    title(['ENRR  NOAA G-4 ',stdtdol(1:12)]);
    if printit
        pname = fullfile(path_plots,['Q_Gradients_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    zz = isfinite(lon);
    figure; plot(lon(zz),lat(zz),'o-',[-180,-140],[0,0],'k--','linewidth',2); grid;
    xlim([min(lon(zz))-1,max(lon(zz))+1]);
    xlabel('Longitude (deg)'); ylabel('Latitude (deg)'); title(['ENRR  NOAA G-4 ',stdtdol(1:12)]);
    if printit
        pname = fullfile(path_plots,['Dropsonde_Map_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    wi_P = -9.8.*rhoaxm.*wi';  % vertical velocity in Pa/s, Wi_P = -g * rohoa * w
    figure;
    subplot(1,2,1); plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(wi*100,paxm,'-','linewidth',2); grid;
%     h1 = plot(wi_P,paxm,'-','linewidth',2); grid;
    hold off; set(gca,'yDir','reverse'); ylim([200 1010]);
    ylabel('P (mb)'); xlabel('Vertical Velocity (cm/s)');
    title(['ENRR  NOAA G-4 ',stdtdol(1:12)]);
    subplot(1,2,2); plot(qaxm,paxm,'-','linewidth',2); grid;
    set(gca,'yDir','reverse'); ylim([200 1010]);
    xlabel('Specific Humidity (g/kg)');
    title(['ENRR  NOAA G-4 ',stdtdol(1:12)]);
    if printit
        pname = fullfile(path_plots,['W_q_Profiles_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    % density weighted wind components
    rhoun = rhoaxm.*unaxm; rhoue = rhoaxm.*ueaxm;
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(aq(3,:)'/1e5.*rhoun,paxm,'-',aq(2,:)'/1e5.*rhoue,paxm,'-',-wi'.*dQdz,paxm,'-','linewidth',2); grid;
    hold off; set(gca,'yDir','reverse'); ylim([200 1010]); xlim([-1e-4 2e-4]);
    ylabel('P (mb)'); xlabel('dQ/dt  g/m^2/s');
    legend(h1,'North grad','East grad','W*dq/dz'); title(['ENRR  NOAA G-4 ',stdtdol(1:12)]);
    if printit
        pname = fullfile(path_plots,['dQ_dt_Profiles_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    Q2 = -(-aq(3,:)'.*unaxm/1e5 - aq(2,:)'.*ueaxm/1e5 + wi'.*dqdz)*2.5;
    Q1 = (-ath(3,:)'.*unaxm/1e5 - ath(2,:)'.*ueaxm/1e5 + wi'.*dthdz);
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(Q2*3600*24,paxm,'-',Q1*3600*24,paxm,'-','linewidth',2); grid;
    hold off; set(gca,'yDir','reverse'); ylim([200,1010]);
    ylabel('P (mb)'); xlabel('Heating Rate (C/day)');
    legend(h1,'Q2','Q1','location','best'); title(['ENRR  NOAA G-4 ',stdtdol(1:12)]);
    if printit
        pname = fullfile(path_plots,['Q1_Q2_Profiles_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    Qr = -1.4;
    Iq = zeros(1200,1);
    for i=2:1200
        Iq(i) = Iq(i-1)+10*(Q1(i-1)-Q2(i-1)-Qr/24/3600);
    end
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(135-rhoaxm(:).*cpa.*Iq, paxm,'-','linewidth',2); grid;
    hold off; set(gca,'yDir','reverse'); ylim([200,1010]);
    ylabel('P (mb)'); xlabel('Moist Static Energy Flux (W/m^2)');
    title(['ENRR  NOAA G-4 ',stdtdol(1:12)]);
    if printit
        pname = fullfile(path_plots,['MSE_Flux_Profile_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    figure;
    subplot(1,2,1); 
    plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(-ap(3,:)'./rhoaxm/1e5*1e2,paxm,'-',fue,paxm,'-',-an(3,:)'.*unaxm/1e5,paxm,'-',-an(2,:)'.*ueaxm/1e5,paxm,'-','linewidth',2);
    hold off; set(gca,'yDir','reverse'); ylim([200,1010]); xlim([-2 2]*1e-4); grid;
    xlabel('dUn/Dt terms'); ylabel('P (mb)');
    legend(h1,'-dP/dn/\rho','-f*ue','un*Dundn','ue*Dunde','location','best');
    subplot(1,2,2);
    plot([0,0],[200,1200],'k-'); hold on;
    h2 = plot(-ap(2,:)'./rhoaxm/1e5*1e2,paxm,'-',-fun,paxm,'-',-ae(3,:)'.*unaxm/1e5,paxm,'-',-ae(2,:)'.*ueaxm/1e5,paxm,'-','linewidth',2);
    hold off; set(gca,'yDir','reverse'); ylim([200,1010]); xlim([-2 2]*1e-4); grid;
    xlabel('dUe/Dt terms'); ylabel('P (mb)');
    legend(h2,'-dP/de/\rho','f*un','un*Duedn','ue*Dende','location','best');
    if printit
        pname = fullfile(path_plots,['dU_dt_Terms_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end

%     close all;
end