% ENRR_sonde_profiles.m
% convective analysis computed with vertical pressure grid

%% initialize
clear
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
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 mean radiative flux divergence, -0.55 °C/day for convective conditions, per Johnson et al, DYNAMO
Qr = -0.55;

%% set desired flight numbers and loop
for ii=[1:5,7:16,18:19] % all wrap, skinny wrap and coathanger flights
    
    %% create output directory and read data
    disp(['Flight RF',sprintf('%02i',ii)]);
    path_plots = ['/Users/byronb/Raw_Data_Files.noindex/ENRR_2016/Fairall_processing/plots/RF',...
                  sprintf('%02.0f',ii),'_BWB_Pgrid'];
    mkdir(path_plots);
    
    % read one flight, drops interpolated to fixed pressure levels
    fname = ['*_RF',sprintf('%02.0f',ii),'_pinterp.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_pinterp.jd); % number of drops on this flight
    
    % read variables for all drops.  Profiles interp to fixed alt grid.
    jd    = snd_pinterp.jd;     % 1xN array, time of drop
    lon   = snd_pinterp.lon;    % 1xN array, lon at start of drop
    lat   = snd_pinterp.lat;    % 1xN array, lat at start of drop
    P_ref = snd_pinterp.P;      % 41x1 array, pressure levels
    zax   = snd_pinterp.z;      % 41xN altitude array, m
    U     = snd_pinterp.U;      % 41xN profile arrays...
    udir  = snd_pinterp.udir;
    Un    = snd_pinterp.Un;     % Un and Ue are 'from' N and E
    Ue    = snd_pinterp.Ue;
    ta    = snd_pinterp.ta;
    rh    = snd_pinterp.rh;
    qa    = snd_pinterp.qa;
%     thasx = (ta + tdk).*(1010./P).^0.287; % air theta, K
    thasx = ta + 0.0098*zax;    % Johnson definition of Theta
    pax   = ones(size(zax)).*P_ref;
    
    stdtdol = datestr(jd(1) + datenum(2015,12,31)); % time of first drop
    
    disp(['Number of profiles= ' num2str(N)])
    
    % set flight specific variables
    switch ii  % set day of year and indices for drops around study area
        case 1
            stjd=22;
            conv_drops = [6:16,18:21,23:27]; % good drops for convective analysis
        case 2
            stjd=26;
            conv_drops = [5:7,9,11:16];
        case 3
            stjd=27;
            conv_drops = 5:29;
        case 4
            stjd=30;
            conv_drops = 4:29;
        case 5
            stjd=31;
            conv_drops = 4:16;
        case 7
            stjd=35;
            conv_drops = 4:22;
        case 8
            stjd=44;
            conv_drops = 4:27;
        case 9
            stjd=46;
            conv_drops = 4:26;
        case 10
            stjd=47;
            conv_drops = 5:24;
        case 11
            stjd=48;
            conv_drops = 5:27;
        case 12
            stjd=49;
            conv_drops = [3:23,25:30];
        case 13
            stjd=50;
            conv_drops = 3:21;
        case 14
            stjd=53;
            conv_drops = 1:N;
        case 15
            stjd=58;
            conv_drops = [6:14,16:21];
        case 16
            stjd=59;
            conv_drops = 6:23;
        case 18
            stjd=61;
            conv_drops = 5:27;
        case 19
            stjd=63;
            conv_drops = [3:11,13:17];
    end
    
    %% latitude gradient contour plots
    latax = ones(41,N).*lat;      % 41xN lat array, using constant lat/lon for each drop
    lonax = ones(41,N).*lon;
    
    [xi,yi] = meshgrid(-2.5:.5:20,1000:-50:200);  % latitude : pressure altitude grid
    jj = isfinite(rh) & isfinite(latax);
    ZI = griddata(latax(jj),pax(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; set(gca,'yDir','reverse');
    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('P (mb)');
    hold off;
    
    jj = find(isfinite(thasx) & isfinite(latax));
    ZI = griddata(latax(jj),pax(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; set(gca,'yDir','reverse');
    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('P (mb)');
    hold off;
    if printit
        pname = fullfile(path_plots,['RH_Theata_LatSection_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    jj = find(isfinite(Un) & isfinite(latax));
    ZI = griddata(latax(jj),pax(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; set(gca,'yDir','reverse');
    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('P (mb)');
    hold off;

    jj = find(isfinite(Ue) & isfinite(latax));
    ZI = griddata(latax(jj),pax(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; set(gca,'yDir','reverse');
    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('P (mb)');
    hold off;
    if printit
        pname = fullfile(path_plots,['Un_Ue_LatSection_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    %% latitude gradients
    
    unaxs = Un;
    ueaxs = Ue;
    qaxs = qa;
    rhoax = pax*100./(Rgas*(ta + tdk).*(1 + 0.61*qaxs/1e3));  % air density
    
    % bin mean profiles by latitude from -2° to 17°
    [~, 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(P_ref==1000);    % pick 1000 mb level
    figure;
    subplot(3,1,1); plot(lat_bins,Yq(pp,:),'o-','linewidth',2); grid;
    ylabel('q (g/kg)'); title(['Mean Surface Gradients: ENRR  NOAA G-4 RF',sprintf('%02.0f',ii)] );
    subplot(3,1,2); plot(lat_bins,Yun(pp,:),'o-',lat_bins,Yue(pp,:),'o-','linewidth',2); grid;
    ylabel('U (m/s)'); legend('Un','Ue','location','best');
    subplot(3,1,3); 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
    
    %% profile analysis regression coefficients
    % 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 = pax(:,conv_drops);
    rhoaxcd = rhoax(:,conv_drops);
    lataxcd = latax(:,conv_drops);
    lonaxcd = lonax(:,conv_drops);
    zaxcd = zax(:,conv_drops);
    
    % convert lat/lon to enu grid points on N/E axes in meter units
    % select axes origin as mean of lat/lon for all 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 regression coefficients w/respect to xE/xN coordinates in meters
    an = NaN(41,3); ae=an; aq=an; ap=an; ath=an; % arrays for coefficients
    zaxm = NaN(41,1); qaxm = zaxm; unaxm = zaxm; ueaxm = zaxm; thaxm = zaxm; rhoaxm = zaxm; % mean profiles
    for dd=1:41
        kk = isfinite(unaxscd(dd,:));  % select 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,:));  % z
        aa = [ones(size(lataxcd(dd,kk)))',xE(kk)',xN(kk)'];
        cc = aa\zaxcd(dd,kk)';
        ap(dd,:) = cc;
        
        % compute means for this level
        rhoaxm(dd) = nanmean(rhoaxcd(dd,:));
        zaxm(dd)   = nanmean(zaxcd(dd,:));
        qaxm(dd)   = nanmean(qaxscd(dd,:));
        unaxm(dd)  = nanmean(unaxscd(dd,:));
        ueaxm(dd)  = nanmean(ueaxscd(dd,:));
        thaxm(dd)  = nanmean(thasxcd(dd,:));
    end 
    
    %% various profiles
    
    % Divergence and Vorticity
    % 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
    dpr = -100*20;                      % Delta press in Pa
    dzr = diff(zaxm); dzr = [155;dzr];  % Delta z in m
    wi2 = zeros(41,1); dwi = wi2;
    for i = 2:41
        wi2(i) = wi2(i-1) - dpr*dv(i); % Vertical velocity, omega, Pa/s
    end
    % wi2 = wi2./rhoaxm/9.83;  % pressure weighted w m/s

    % taper wi to zero at 100mb level (16.2 km)
    wim2 = wi2(41);
    dwi = wim2./100.^2.*P_ref.^2; % taper to zero at 100 mb
    wi2b = wi2 - dwi;

    wix = wi2;  % use omega
    
    % vertical gradient of specific humidty, g/kg/mb
    dqdp = diff(qaxm)/dpr; dqdp = [dqdp;0];

    % vertical gradient of theta deg/mb
    dthdp = diff(thaxm)/dpr; dthdp = [dthdp;0];

    % vertical gradient of water vapor density, g/m3/mb
    dQdp = rhoaxm.*dqdp;

    % vertical Un, Ue derivatives, m/s/mb
    dundp = diff(unaxm)/dpr; dundp = [dundp;0];
    duedp = diff(ueaxm)/dpr; duedp = [duedp;0];

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

    % here Q1 is actually Q1/Cp & Q2 isQ2/Le*Le/Cp, Q/Cp is deg/s units,  2.5 = Le/Cp
    Q1 = (-ath(:,3).*unaxm - ath(:,2).*ueaxm + wix.*dthdp);
    Q2 = -(-aq(:,3).*unaxm - aq(:,2).*ueaxm + wix.*dqdp)*2.5;
    
    %% total integrals
    Q1b = -cpa/9.85*nansum(dpr*Q1);
    Qr = -1 + 1.35e-3*Q1b;  % from CF analysis of Johnson, DYNAMO paper
    Q2b = -cpa/9.85*nansum(dpr*Q2);
    Qrb = -cpa/9.85*41*dpr*Qr/24/3600;
    Iq = -cpa/9.85*nansum(dpr*(Q1-Q2-Qr/24/3600)); % Q1-Q2-Qr = Hs+Hl
    wib = -nansum(dpr*wix)/800/100/100*3600;       % mean omega in mb/hr
    Prec = (Q1b-Qrb-10)/2500/1e3*24*3600;          % precipitation, mm/day

    Ints = [Q1b;Q2b;Qrb;Iq;wib;Prec];

    disp('Q1b      Q2b      Qrb      Iq       wi2       Precp');
    fprintf('%06.2f   %06.2f   %06.2f   %06.2f   %06.5f  %06.2f\n', Ints);
    
    %% compute moist static energy
    % first, compute mean variables for study area from SURFA model fields
    % out = [u,ta,rh,P,ts,rs,rl,qa,hs,hl]
    surfa_nc_path = '/Users/byronb/Raw_Data_Files.noindex/ENRR_2016/SURFA_nc';
    out = SURFA_avg(stjd,surfa_nc_path,lon,lat,conv_drops);
    
    % Integral of DeltaQ, W/m2
    I_MSE = zeros(41,1);
    for i = 2:41
        %  W/m2 = W/m2 + kg/m/s2 * W*s/kg/K * s2/m * K/s 
        I_MSE(i) = I_MSE(i-1) - dpr*cpa/9.85*(Q1(i-1)-Q2(i-1)-Qr/24/3600);
    end
    MSE0 = out(9) + out(10);  % surface MSE = hs0 + hl0 from SURFA
    
    %% compare SURFA heat fluxes to COARE 3.5
    % A=[usr tau hsb hlb hbb hsbb hlwebb tsr qsr zot zoq Cd Ch Ce  L zet dter dqer tkt Urf Trf Qrf RHrf UrfN Rnl Le rhoa UN U10 U10N Cdn_10 Chn_10 Cen_10 RF Qs Evap T10 Q10 RH10];
    %     1   2   3   4   5   6    7      8   9  10  11  12 13 14 15 16   17   18   19  20  21  22  23   24   25 26  27  28  29  30    31     32     33   34 35  36  37  38   39
    C35 = coare35vn(out(1),10,out(2),2,out(3),2,out(4),out(5),out(6),out(7),mean(lat(conv_drops)),2000,0,NaN,NaN);
    disp(['Sensible Heat: SURFA = ',sprintf('%2.2f',out(9)),'   COARE = ',sprintf('%2.2f',C35(3))]);
    disp(['Latent Heat:   SURFA = ',sprintf('%3.2f',out(10)),'  COARE = ',sprintf('%3.2f',C35(4))]);
    
    Ints = [Ints;out(9);out(10);C35(3);C35(4)];
    save(fullfile(path_plots,['Integrals_RF',sprintf('%02.0f',ii),'.mat']),'Ints');
    
    %% plots
    % divergence profile
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(-(an(:,3) + ae(:,2)), P_ref,'-','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)]);
    if printit
        pname = fullfile(path_plots,['Divergence_Profile_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end

    % curl profile
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot((ae(:,3) - an(:,2)), P_ref,'-','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)]);
    if printit
        pname = fullfile(path_plots,['Curl_Profile_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end

    % N/E q gradient profiles
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(aq(:,3),P_ref,'-',aq(:,2),P_ref,'-','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)]);
    if printit
        pname = fullfile(path_plots,['Q_Gradients_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end

    % maps of drop sonde locations
    zz = isfinite(lon);
    figure; plot(lon(zz),lat(zz),'bo-',lon,lat,'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)]);
    if printit
        pname = fullfile(path_plots,['Dropsonde_Map_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    figure; plot(xE/1000,xN/1000,'bd','linewidth',2); hold on;
    plot([-1000,1000],[0,0],'k-',[0,0],[-1000,1000],'k-'); grid;
    axis([min(xE/1000)-100,max(xE/1000)+100,min(xN/1000)-100,max(xN/1000)+100]);
    xlabel('Kilometers East'); ylabel('Kilometers North');
    title(['ENRR  NOAA G-4 RF',sprintf('%02.0f',ii)]);
    if printit
        pname = fullfile(path_plots,['Dropsonde_xGrid_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    % W and q profiles
    figure;
    subplot(1,2,1); plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(wi2/100*3600,P_ref,'b-','linewidth',2); grid;
    hold off; set(gca,'yDir','reverse'); ylim([200 1010]);
    ylabel('P (mb)'); xlabel('Omega mb/hr');
    title('ENRR  NOAA G-4');
    subplot(1,2,2); plot(qaxm,P_ref,'-','linewidth',2); grid;
    set(gca,'yDir','reverse'); ylim([200 1010]);
    xlabel('Specific Humidity (g/kg)');
    title(['RF',sprintf('%02.0f',ii)]);
    if printit
        pname = fullfile(path_plots,['W_q_Profiles_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    % Q1/Q2 profiles
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(Q2*3600*24,P_ref,'-',Q1*3600*24,P_ref,'-','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)]);
    if printit
        pname = fullfile(path_plots,['Q1_Q2_Profiles_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    % MSE profile =  MSE0 - I_MSE
    figure; plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(MSE0-I_MSE, P_ref,'-','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)]);
    if printit
        pname = fullfile(path_plots,['MSE_Flux_Profile_RF',sprintf('%02.0f',ii),'.png']);
        print(pname,'-dpng');
    end
    
    % dU/dt terms
    figure;
    subplot(1,2,1); 
    plot([0,0],[200,1200],'k-'); hold on;
    h1 = plot(fue,P_ref,'-',-an(:,3).*unaxm,P_ref,'-',-an(:,2).*ueaxm,P_ref,'-','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,'-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(-fun,P_ref,'-',-ae(:,3).*unaxm,P_ref,'-',-ae(:,2).*ueaxm,P_ref,'-','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,'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