% ENRR_sonde_profiles.m

clear all
close all
printit = true;

% 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/data2';
cd(path_data);

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

% set desired flight numbers and loop
for ii=[1:5,7:16,18:19] % 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),'_BWB'];
    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
    
    % 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;     % Un and Ue are 'from' N and E
    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)])
    
     % set flight specific variables
    switch ii
        case 1
            conv_drops = [6:16,18:21,23:27]'; % good drops for convective analysis
        case 2
            conv_drops = [5:7,9,11:16];
        case 3
            conv_drops = 5:29;
        case 4
            conv_drops = 4:29;
        case 5
            conv_drops = 4:16;
        case 7
            conv_drops = 4:22;
        case 8
            conv_drops = 4:27;
        case 9
            conv_drops = 4:26;
        case 10
            conv_drops = 5:24;
        case 11
            conv_drops = 5:27;
        case 12
            conv_drops = [3:23,25:30];
        case 13
            conv_drops = 3:21;
        case 14
            conv_drops = 1:N;
        case 15
            conv_drops = [6:14,16:21];
        case 16
            conv_drops = 6:23;
        case 18
            conv_drops = 5:27;
        case 19
            conv_drops = [3:11,13:17];
    end
    
    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 (%) ','RF',sprintf('%02.0f',ii),' 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) ','RF',sprintf('%02.0f',ii),' 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) ','RF',sprintf('%02.0f',ii),' 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) ','RF',sprintf('%02.0f',ii),' 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 for convective analysis 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); Yp = Yp';
    [~, Yq,  ~, ~, ~] = binave2(lat,qaxs,2,-2,17); Yq = Yq';
    [~, Yun, ~, ~, ~] = binave2(lat,ueaxs,2,-2,17); Yun = Yun';
    [~, Yue, ~, ~, ~] = binave2(lat,unaxs,2,-2,17); Yue = Yue';
    [~, Yh,  ~, ~, ~] = binave2(lat,rh,2,-2,17); Yh = Yh';
    [lat_bins, Yth, ~, ~, ~] = binave2(lat,thasx,2,-2,17); Yth = Yth';
    
    % 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 ','RF',sprintf('%02.0f',ii),' ',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
    
    % set up arrays for convective analysis
    % first, select only valid drops for convective analysis
    unaxscd = unaxs(:,conv_drops);
    ueaxscd = ueaxs(:,conv_drops);
    qaxscd = qaxs(:,conv_drops);
    thasxcd = thasx(:,conv_drops);
    paxscd = paxs(:,conv_drops);
    rhoaxcd = rhoax(:,conv_drops);
    lataxcd = latax(:,conv_drops);
    lonaxcd = lonax(:,conv_drops);
    
    % convert lat/lon to enu grid points on N/E axis in meters
    % select origin as mean of lat/lon for convective drops
    lon0 = mean(lon(conv_drops));
    lat0 = mean(lat(conv_drops));
    wgs84 = wgs84Ellipsoid('meters');  % specifies ellipsoid and units...
    [xE,xN,~] = geodetic2enu(lat(conv_drops),lon(conv_drops),0,lat0,lon0,0,wgs84);
    
    % compute dX/dE and dX/dN at all levels for Ue, Un, q, P and theta
%     an = NaN(1200,2); ae=an; aq=an; ap=an; ath=an; % arrays for [dX/dE(ii),dX/dN(ii)]
%     paxm = NaN(1200,1); qaxm = paxm; unaxm = paxm; ueaxm = paxm; thaxm = paxm; rhoaxm = paxm;
%     for dd=1:1200
%         gg = isfinite(unaxscd(dd,:));
%         pcoefE = polyfit(xE(gg),unaxscd(dd,gg),1);
%         pcoefN = polyfit(xN(gg),unaxscd(dd,gg),1);
%         an(dd,:) = [pcoefE(1),pcoefN(1)];
%         gg = isfinite(ueaxscd(dd,:));
%         pcoefE = polyfit(xE(gg),ueaxscd(dd,gg),1);
%         pcoefN = polyfit(xN(gg),ueaxscd(dd,gg),1);
%         ae(dd,:) = [pcoefE(1),pcoefN(1)];
%         gg = isfinite(qaxscd(dd,:));
%         pcoefE = polyfit(xE(gg),qaxscd(dd,gg),1);
%         pcoefN = polyfit(xN(gg),qaxscd(dd,gg),1);
%         aq(dd,:) = [pcoefE(1),pcoefN(1)];
%         gg = isfinite(thasxcd(dd,:));
%         pcoefE = polyfit(xE(gg),thasxcd(dd,gg),1);
%         pcoefN = polyfit(xN(gg),thasxcd(dd,gg),1);
%         ath(dd,:) = [pcoefE(1),pcoefN(1)];
%         gg = isfinite(paxscd(dd,:));
%         pcoefE = polyfit(xE(gg),paxscd(dd,gg),1);
%         pcoefN = polyfit(xN(gg),paxscd(dd,gg),1);
%         ap(dd,:) = [pcoefE(1),pcoefN(1)];
%         
%         % compute means for this level
%         rhoaxm(dd) = nanmean(rhoaxcd(dd,:));
%         paxm(dd)   = nanmean(paxscd(dd,:));
%         qaxm(dd)   = nanmean(qaxscd(dd,:));
%         unaxm(dd)  = nanmean(unaxscd(dd,:));
%         ueaxm(dd)  = nanmean(ueaxscd(dd,:));
%         thaxm(dd)  = nanmean(thasxcd(dd,:));
%     end
    
    % compute regression coefficients w/respect to xE/xN coordinates in meters
    an = NaN(1200,3); ae=an; aq=an; ap=an; ath=an; % arrays for coefficients
    paxm = NaN(1200,1); qaxm = paxm; unaxm = paxm; ueaxm = paxm; thaxm = paxm; rhoaxm = paxm; % mean profiles
    for dd=1:1200
        kk = isfinite(unaxscd(dd,:));  % finite Un values at level dd
        aa = [ones(size(lataxcd(dd,kk)))',xE(kk)',xN(kk)'];
        cc = aa\unaxscd(dd,kk)';
        an(dd,:) = cc;
        kk = isfinite(ueaxscd(dd,:));  % Ue
        aa = [ones(size(lataxcd(dd,kk)))',xE(kk)',xN(kk)'];
        cc = aa\ueaxscd(dd,kk)';
        ae(dd,:) = cc;
        kk = isfinite(qaxscd(dd,:));  % q
        aa = [ones(size(lataxcd(dd,kk)))',xE(kk)',xN(kk)'];
        cc = aa\qaxscd(dd,kk)';
        aq(dd,:) = cc;
        kk = isfinite(thasxcd(dd,:)); % Theta
        aa = [ones(size(lataxcd(dd,kk)))',xE(kk)',xN(kk)'];
        cc = aa\thasxcd(dd,kk)';
        ath(dd,:) = cc;
        kk = isfinite(paxscd(dd,:));  % P
        aa = [ones(size(lataxcd(dd,kk)))',xE(kk)',xN(kk)'];
        cc = aa\paxscd(dd,kk)';
        ap(dd,:) = cc;
        
        % compute means for this level
        rhoaxm(dd) = nanmean(rhoaxcd(dd,:));
        paxm(dd)   = nanmean(paxscd(dd,:));
        qaxm(dd)   = nanmean(qaxscd(dd,:));
        unaxm(dd)  = nanmean(unaxscd(dd,:));
        ueaxm(dd)  = nanmean(ueaxscd(dd,:));
        thaxm(dd)  = nanmean(thasxcd(dd,:));
    end 
    
    % dv = dUn/dN + dUe/dE, cv = dUn/dE - dUe/dN
    % winds are 'from' so negate
    dv = -(an(:,3) + ae(:,2));  % horizontal divergence
    cl = -(an(:,2) - ae(:,3));  % vertical vorticity
    
    % compute vertical velocity by integrating horizontal divergence
    wi = zeros(1200,1);
    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
    fun = nanmean((unaxscd.*sin(pi/180*lataxcd)),2)*1.45e-4;
    fue = nanmean((ueaxscd.*sin(pi/180*lataxcd)),2)*1.45e-4;
    
    % plots
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(-(an(:,3) + ae(:,2)), 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 ','RF',sprintf('%02.0f',ii),' ',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)), 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 ','RF',sprintf('%02.0f',ii),' ',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),paxm,'-',aq(:,2),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 ','RF',sprintf('%02.0f',ii),' ',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),'bo-',lon(conv_drops),lat(conv_drops),'rx',[-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 ','RF',sprintf('%02.0f',ii),' ',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 ','RF',sprintf('%02.0f',ii),' ',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).*rhoun,paxm,'-',aq(:,2).*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 ','RF',sprintf('%02.0f',ii),' ',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 - aq(:,2).*ueaxm + wi.*dqdz)*2.5;
    Q1 = (-ath(:,3).*unaxm - ath(:,2).*ueaxm + 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 ','RF',sprintf('%02.0f',ii),' ',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 ','RF',sprintf('%02.0f',ii),' ',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*1e2,paxm,'-',fue,paxm,'-',-an(:,3).*unaxm,paxm,'-',-an(:,2).*ueaxm,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');
    title(['ENRR  NOAA G-4 ']);
    subplot(1,2,2);
    plot([0,0],[200,1200],'k-'); hold on;
    h2 = plot(-ap(:,2)./rhoaxm*1e2,paxm,'-',-fun,paxm,'-',-ae(:,3).*unaxm,paxm,'-',-ae(:,2).*ueaxm,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');
    title(['RF',sprintf('%02.0f',ii),' ',stdtdol(1:12)]);
    if printit
        pname = fullfile(path_plots,['dU_dt_Terms_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end

%      close all;
end