% plot clouds & G4 track for specified flight

%% initialize

close all
clearvars

printit = true;

% set flight number
RF = 18;

switch RF
    case 1
        stjd=22;
    case 2
        stjd=26;
    case 3
        stjd=27;
    case 4
        stjd=30;
    case 5
        stjd=31;
    case 7
        stjd=35;
    case 8
        stjd=44;
    case 9
        stjd=46;
    case 10
        stjd=47;
    case 11
        stjd=48;
    case 12
        stjd=49;
    case 13
        stjd=50;
    case 14
        stjd=53;
    case 15
        stjd=58;
    case 16
        stjd=59;
    case 18
        stjd=61;
    case 19
        stjd=63;
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 image

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

x_lim = [lon_lo,lon_hi];
y_lim = [lat_hi,lat_lo];
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-','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-','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-','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-','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-','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_tau RF',sprintf('%02.0f',RF),'.png']);
    print(pname,'-dpng');
end


