% plot SURFA fields & G4 track for specified flight
% compute mean SURFA data for convective study area defined
%   by polygon of dropsondes surrounding study area.

%% initialize

close all
clearvars

printit = false;

% set flight number
RF = 1;

switch RF  % 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


% 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 to surfa nc files
surfa_nc_path = '/Users/byronb/Raw_Data_Files.noindex/ENRR_2016/SURFA_nc';

path_plots = ['/Users/byronb/Raw_Data_Files.noindex/ENRR_2016/Fairall_processing/plots/RF',sprintf('%02.0f',RF),'_BWB'];
mkdir(path_plots);  % create directory for figures if it doesn't exist

% dropsonde data path
path_G4_data = '/Users/byronb/Raw_Data_Files.noindex/ENRR_2016/Fairall_processing/data2';

graphformat = '.png';
graphdevice = '-dpng';

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

%% load data

% loads structure with 3-hr model fields for requested date range
% path to data is defined in read_surfa_ecmwf.m
% fields: dfl,jd,lat,lon,evap,lowcld,midcld,wcloud,P,ts,tdew,ta,ux,uy,u,tau,hs,hl,rs,rl,rs0,rl0,rsn,rln
sfa = ENRR_read_surfa_ecmwf(surfa_nc_path,stjd,stjd);  % one day only

% read one flight, drops interpolated to fixed altitude vector
fname = ['*_RF',sprintf('%02.0f',RF),'_zinterp.mat'];
G4file = dir(fullfile(path_G4_data,fname));
% structure with fields: jd, lon, lat, z, U, udir, Un, Uw, ta, rh, P, qa
load(fullfile(path_G4_data,G4file.name));

% 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
lon = lon+360;

%% plot images and flight track

% define lat/lon bounds for complete flight track and
% find SURFA indices corresponding to these limits
lat_lo = -2;        lat_lo_idx = find(sfa.lat==lat_lo);
lat_hi = 20;        lat_hi_idx = find(sfa.lat==lat_hi);
lon_lo = floor(min(lon))-1;   lon_lo_idx = find(sfa.lon==lon_lo);
lon_hi = ceil(max(lon))+1;   lon_hi_idx = find(sfa.lon==lon_hi);

% image limit arrays
x_lim = [lon_lo,lon_hi];
y_lim = [lat_hi,lat_lo];  % image is plotted from upper right corner...

N = length(sfa.jd);

figure;
dstr = sfa.dfl{1}(end-15:end-3);
imagesc(x_lim,y_lim,sfa.midcld(lon_lo_idx:lon_hi_idx,lat_hi_idx:lat_lo_idx,1)');
axis xy; grid; hold on;
c = colorbar; c.Label.String = 'Mid Cloud Fraction, %';
caxis([0,100]); c.Ticks = 0:20:100;
xxx = isfinite(lat) & isfinite(lon);
plot(lon(xxx),lat(xxx),'ro-',lon(conv_drops),lat(conv_drops),'yx','linewidth',2); hold off;
ylabel('latitiude'); xlabel('longitude');
title(['DOY ',sprintf('%07.3f',sfa.jd(1)),'  SURFA-ECMWF Mid Cloud, G4-RF',sprintf('%02.0f',RF),'   ',dstr(1:4),...
    ' ',dstr(5:6),' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
if printit
    pname = fullfile(path_plots,['G4_track_SURFA_midCld_RF',sprintf('%02.0f',RF),'.png']);
    print(pname,'-dpng');
end

figure;
dstr = sfa.dfl{1}(end-15:end-3);
imagesc(x_lim,y_lim,sfa.wcloud(lon_lo_idx:lon_hi_idx,lat_hi_idx:lat_lo_idx,1)');
axis xy; grid; hold on;
c = colorbar; c.Label.String = 'Cloud Water, g/m^2';
caxis([0,600]); c.Ticks = 0:100:600;
xxx = isfinite(lat) & isfinite(lon);
plot(lon(xxx),lat(xxx),'ro-',lon(conv_drops),lat(conv_drops),'yx','linewidth',2); hold off;
ylabel('latitiude'); xlabel('longitude');
title(['DOY ',sprintf('%07.3f',sfa.jd(1)),'  SURFA-ECMWF Cloud Water, G4-RF',sprintf('%02.0f',RF),'   ',dstr(1:4),...
    ' ',dstr(5:6),' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
if printit
    pname = fullfile(path_plots,['G4_track_SURFA_CldWtr_RF',sprintf('%02.0f',RF),'.png']);
    print(pname,'-dpng');
end

figure;
dstr = sfa.dfl{1}(end-15:end-3);
imagesc(x_lim,y_lim,sfa.hs(lon_lo_idx:lon_hi_idx,lat_hi_idx:lat_lo_idx,1)');
axis xy; grid; hold on;
c = colorbar; c.Label.String = 'Sens Heat, W/m^2';
caxis([-100,200]); c.Ticks = -100:50:200;
xxx = isfinite(lat) & isfinite(lon);
plot(lon(xxx),lat(xxx),'ro-',lon(conv_drops),lat(conv_drops),'yx','linewidth',2); hold off;
ylabel('latitiude'); xlabel('longitude');
title(['DOY ',sprintf('%07.3f',sfa.jd(1)),'  SURFA-ECMWF Sensible Heat Flux, G4-RF',sprintf('%02.0f',RF),'   ',dstr(1:4),...
    ' ',dstr(5:6),' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
if printit
    pname = fullfile(path_plots,['G4_track_SURFA_hs_RF',sprintf('%02.0f',RF),'.png']);
    print(pname,'-dpng');
end

figure;
dstr = sfa.dfl{1}(end-15:end-3);
imagesc(x_lim,y_lim,sfa.hl(lon_lo_idx:lon_hi_idx,lat_hi_idx:lat_lo_idx,1)');
axis xy; grid; hold on;
c = colorbar; c.Label.String = 'Latent Heat, W/m^2';
caxis([-50,400]); c.Ticks = -50:50:400;
xxx = isfinite(lat) & isfinite(lon);
plot(lon(xxx),lat(xxx),'ro-',lon(conv_drops),lat(conv_drops),'yx','linewidth',2); hold off;
ylabel('latitiude'); xlabel('longitude');
title(['DOY ',sprintf('%07.3f',sfa.jd(1)),'  SURFA-ECMWF Latent Heat Flux, G4-RF',sprintf('%02.0f',RF),'   ',dstr(1:4),...
    ' ',dstr(5:6),' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
if printit
    pname = fullfile(path_plots,['G4_track_SURFA_hl_RF',sprintf('%02.0f',RF),'.png']);
    print(pname,'-dpng');
end

figure;
dstr = sfa.dfl{1}(end-15:end-3);
imagesc(x_lim,y_lim,sfa.tau(lon_lo_idx:lon_hi_idx,lat_hi_idx:lat_lo_idx,1)');
axis xy; grid; hold on;
c = colorbar; c.Label.String = 'Stress, N/m^2';
caxis([0,0.2]); c.Ticks = 0:0.05:0.2;
xxx = isfinite(lat) & isfinite(lon);
plot(lon(xxx),lat(xxx),'ro-',lon(conv_drops),lat(conv_drops),'yx','linewidth',2); hold off;
ylabel('latitiude'); xlabel('longitude');
title(['DOY ',sprintf('%07.3f',sfa.jd(1)),'  SURFA-ECMWF Wind Stress, G4-RF',sprintf('%02.0f',RF),'   ',dstr(1:4),...
    ' ',dstr(5:6),' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
if printit
    pname = fullfile(path_plots,['G4_track_SURFA_ta_RF',sprintf('%02.0f',RF),'.png']);
    print(pname,'-dpng');
end

figure;
dstr = sfa.dfl{1}(end-15:end-3);
imagesc(x_lim,y_lim,sfa.evp(lon_lo_idx:lon_hi_idx,lat_hi_idx:lat_lo_idx,1)');
axis xy; grid; hold on;
c = colorbar; c.Label.String = 'Evaporation, g/m^2';
caxis([0,0.14]); c.Ticks = 0:0.02:0.14;
xxx = isfinite(lat) & isfinite(lon);
plot(lon(xxx),lat(xxx),'ro-',lon(conv_drops),lat(conv_drops),'yx','linewidth',2); hold off;
ylabel('latitiude'); xlabel('longitude');
title(['DOY ',sprintf('%07.3f',sfa.jd(1)),'  SURFA-ECMWF Evaporation, G4-RF',sprintf('%02.0f',RF),'   ',dstr(1:4),...
    ' ',dstr(5:6),' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
if printit
    pname = fullfile(path_plots,['G4_track_SURFA_evap_RF',sprintf('%02.0f',RF),'.png']);
    print(pname,'-dpng');
end

Le = (2.501-.00237*25)*1e3;  % approximate Le

figure;
dstr = sfa.dfl{1}(end-15:end-3);
imagesc(x_lim,y_lim,sfa.hl(lon_lo_idx:lon_hi_idx,lat_hi_idx:lat_lo_idx,1)'/Le);
axis xy; grid; hold on;
c = colorbar; c.Label.String = 'Est Evaporation, g/m^2';
caxis([0,0.14]); c.Ticks = 0:0.02:0.14;
xxx = isfinite(lat) & isfinite(lon);
plot(lon(xxx),lat(xxx),'ro-',lon(conv_drops),lat(conv_drops),'yx','linewidth',2); hold off;
ylabel('latitiude'); xlabel('longitude');
title(['DOY ',sprintf('%07.3f',sfa.jd(1)),'  SURFA-ECMWF HL/Le, G4-RF',sprintf('%02.0f',RF),'   ',dstr(1:4),...
    ' ',dstr(5:6),' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
if printit
    pname = fullfile(path_plots,['G4_track_SURFA_EstEvap_RF',sprintf('%02.0f',RF),'.png']);
    print(pname,'-dpng');
end

%% compute mean SURFA field values within study area

% define vertices of polygon
Xlon = lon(conv_drops)'; Ylat = lat(conv_drops)';

% duplicate fisrt vertex at end of array to complete polygon
Xlon = [Xlon; Xlon(1)]; Ylat = [Ylat; Ylat(1)];

% build array of SURFA lat/lon indices
Nlon = length(sfa.lon(lon_lo_idx:lon_hi_idx)); Nlat = length(sfa.lat(lat_hi_idx:lat_lo_idx));
SFAlon = repmat(sfa.lon(lon_lo_idx:lon_hi_idx),Nlat,1);
SFAlat = NaN(Nlon*Nlat,1);
for ii = 1:Nlat
    SFAlat(Nlon*(ii-1)+1:Nlon*ii) = sfa.lat(lat_hi_idx+ii-1);
end

% find SURFA indices that fall within the study area polygon
idx = inpolygon(SFAlon,SFAlat,Xlon,Ylat);
lon_valid = SFAlon(idx); lat_valid = SFAlat(idx);
lon_idx = NaN(length(lon_valid),1); lat_idx = NaN(length(lat_valid),1);
for ii = 1:length(lon_valid)
    lon_idx(ii) = find(sfa.lon==lon_valid(ii));
    lat_idx(ii) = find(sfa.lat==lat_valid(ii));
end

% test plot to show SURFA grid cells included in means...
figure;
dstr = sfa.dfl{1}(end-15:end-3);
imagesc(x_lim,y_lim,sfa.midcld(lon_lo_idx:lon_hi_idx,lat_hi_idx:lat_lo_idx,1)');
axis xy; grid; hold on;
c = colorbar; c.Label.String = 'Mid Cloud Fraction, %';
caxis([0,100]); c.Ticks = 0:20:100;
xxx = isfinite(lat) & isfinite(lon);
plot(lon(xxx),lat(xxx),'ro-',lon_valid,lat_valid,'kx','linewidth',2); hold off;
ylabel('latitiude'); xlabel('longitude');
title(['DOY ',sprintf('%07.3f',sfa.jd(1)),'  SURFA-ECMWF Mean Grids, G4-RF',sprintf('%02.0f',RF),'   ',dstr(1:4),...
    ' ',dstr(5:6),' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
if printit
    pname = fullfile(path_plots,['G4_track_SURFA_MeanGrids_RF',sprintf('%02.0f',RF),'.png']);
    print(pname,'-dpng');
end


% compute means for SURFA field values within polygon
% mid-level cloud percentage
temp = NaN(length(lon_valid),1);
for ii = 1:length(lon_valid)
    temp(ii) = sfa.midcld(lon_idx(ii),lat_idx(ii),1);
end
MeanMidCld = mean(temp);

% cloud water content, g/m2
temp = NaN(length(lon_valid),1);
for ii = 1:length(lon_valid)
    temp(ii) = sfa.wcloud(lon_idx(ii),lat_idx(ii),1);
end
MeanCldWater = mean(temp);

% sensible heat, W/m2
temp = NaN(length(lon_valid),1);
for ii = 1:length(lon_valid)
    temp(ii) = sfa.hs(lon_idx(ii),lat_idx(ii),1);
end
MeanSensHeat = mean(temp);

% latent heat, W/m2
temp = NaN(length(lon_valid),1);
for ii = 1:length(lon_valid)
    temp(ii) = sfa.hl(lon_idx(ii),lat_idx(ii),1);
end
MeanLatHeat = mean(temp);

% wind stress, N/m2
temp = NaN(length(lon_valid),1);
for ii = 1:length(lon_valid)
    temp(ii) = sfa.tau(lon_idx(ii),lat_idx(ii),1);
end
MeanTau = mean(temp);

% evaporation, kg/m2
temp = NaN(length(lon_valid),1);
for ii = 1:length(lon_valid)
    temp(ii) = sfa.evp(lon_idx(ii),lat_idx(ii),1);
end
MeanEvap = mean(temp);


