% plot Calwater 2015 SURFA and ship data

%% initialize

close all
clearvars

% set DOY range
stjd = 21;
endjd = 21;

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

save_plots = true;  % set 'true' to save plots
run_movies = true;  % set 'true' to create mp4 movies of SURFA model fields
surfa_fig_path = '/Users/byronb/Raw_Data_Files.noindex/ENRR_2016/SURFA_plots';
mkdir(surfa_fig_path);  % create directory for figures if it doesn't exist

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

%% 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,endjd);

% load filtered hourly flux file from RHB - this file should have a header line defining variable names
% variables: jd,ushp,wspd,wdir,rwspd,rwdir,hed,ts,ta,qs,qa,hsc,hsib,hsb,hlc,hlib,hlb,tauc,taux,tauib,taub
%        rs,rl,rain,jplume,tilt,jmanuv,ct,cq,cu,cw,hrain,hwebb,lat,lon,ushp2,wspd_im,wdir_im
%        cog,tsg,ta_im,qs_im,qa_im,rs_im,rl_im,qlic,std_qlic,co2,std_co2,pmb,wspd2,wdir2,rlnet,hnetb,hnet_turb
%import_hourly_flux

% air constants for later use
% tdk = 273.15;
% Q = qa/1e3;     % kg/kg
% Beta = 1.2;
% von  = 0.4;     % von karman constnt
% grav = 9.83; 
% Rgas = 287.1;
% Le   = (2.501-.00237*ts)*1e6;  % latent heat vaporization of water
% cpa  = 1004.67;                % Specific heat dry air
% cpv  = cpa*(1+0.84*Q);         % Specific heat moist air
% rhoa = pmb*100./(Rgas*(ta+tdk).*(1+0.61*Q));  % air density
% visa = 1.326e-5*(1+6.542e-3.*ta+8.301e-6*ta.^2-4.84e-9*ta.^3); % air viscosity

%% extract 3-Hr ship and surfa time series at location of RHB

% define empty arrays for 3-hr surfa time series
% sfa_P = NaN(size(sfa.jd)); sfa_ts = sfa_P; sfa_ta = sfa_P;
% sfa_u = sfa_P; sfa_tau = sfa_P; sfa_hs = sfa_P;
% sfa_hl = sfa_P; sfa_rs = sfa_P; sfa_rl = sfa_P; sfa_qa = sfa_P;

% define empty arrays for 3-hr ship time series
% rhb_P = NaN(size(sfa.jd)); rhb_ts = rhb_P; rhb_ta = rhb_P; rhb_u = rhb_P;
% rhb_tauc = rhb_P; rhb_tauib = rhb_P; rhb_taub = rhb_P; rhb_qa = rhb_P;
% rhb_hsc = rhb_P; rhb_hsib = rhb_P; rhb_hsb = rhb_P;
% rhb_hlc = rhb_P; rhb_hlib = rhb_P; rhb_hlb = rhb_P;
% rhb_rs = rhb_P; rhb_rl = rhb_P; rhb_lat = rhb_P; rhb_lon = rhb_P;

% adjust timestamps to center of averaging intervals
%    sfa timestamp is the end of a 3-hr interval
%    psd timestamp is the start of a 1-hr interval
% sfa.jd = sfa.jd - 1.5/24;
% jd = jd + 0.5/24;

% for ii = 1:length(sfa.jd)
%     % ship data nearest 3-hr surfa forcast time - sfa.jd(ii)
%     jd_idx = find(abs(jd - sfa.jd(ii))<1.5/24);
%     % sfa latitude index for ship location at sfa.jd(ii)
%     [~,lat_idx] = min(abs(sfa.lat - nanmean1(lat(jd_idx))));
%     % sfa longitude index for ship location at sfa.jd(ii)
%     [~,lon_idx] = min(abs(sfa.lon - nanmean1(lon(jd_idx))));
%     
%     % extract surfa time series at ship location
%     sfa_P(ii) = sfa.P(lon_idx,lat_idx,ii);
%     sfa_ts(ii) = sfa.ts(lon_idx,lat_idx,ii);
%     sfa_ta(ii) = sfa.ta(lon_idx,lat_idx,ii);
%     sfa_u(ii) = sfa.u(lon_idx,lat_idx,ii);
%     sfa_tau(ii) = sfa.tau(lon_idx,lat_idx,ii);
%     sfa_hs(ii) = sfa.hs(lon_idx,lat_idx,ii);
%     sfa_hl(ii) = sfa.hl(lon_idx,lat_idx,ii);
%     sfa_rs(ii) = sfa.rs(lon_idx,lat_idx,ii);
%     sfa_rl(ii) = sfa.rl(lon_idx,lat_idx,ii);
%     sfa_qa(ii) = sfa.qa(lon_idx,lat_idx,ii);
%     
%     % compute mean ship result on surfa timebase
%     rhb_P(ii) = nanmean1(pmb(jd_idx));
%     rhb_ts(ii) = nanmean1(ts(jd_idx));
%     rhb_ta(ii) = nanmean1(ta(jd_idx));
%     rhb_u(ii) = nanmean1(wspd(jd_idx));
%     rhb_tauc(ii) = nanmean1(tauc(jd_idx));
%     rhb_tauib(ii) = nanmean1(tauib(jd_idx));
%     rhb_taub(ii) = nanmean1(taub(jd_idx));
%     rhb_hsc(ii) = nanmean1(hsc(jd_idx));
%     rhb_hsib(ii) = nanmean1(hsib(jd_idx));
%     rhb_hsb(ii) = nanmean1(hsb(jd_idx));
%     rhb_hlc(ii) = nanmean1(hlc(jd_idx));
%     rhb_hlib(ii) = nanmean1(hlib(jd_idx));
%     rhb_hlb(ii) = nanmean1(hlb(jd_idx));
%     rhb_rs(ii) = nanmean1(rs(jd_idx));
%     rhb_rl(ii) = nanmean1(rl(jd_idx));
%     rhb_lat(ii) = nanmean1(lat(jd_idx));
%     rhb_lon(ii) = nanmean1(lon(jd_idx));
%     rhb_qa(ii) = nanmean1(qa(jd_idx));
% end;

%% run comparison plots

% one2one = [-10000,10000];
% 
% % met comparison: qa, ta, ts, U
% figure();
% qamin = min(min(sfa_qa),min(rhb_qa))-1; qamax = max(max(sfa_qa),max(rhb_qa))+1;
% tamin = min(min(sfa_ta),min(rhb_ta))-1; tamax = max(max(sfa_ta),max(rhb_ta))+1;
% tsmin = min(min(sfa_ts),min(rhb_ts))-1; tsmax = max(max(sfa_ts),max(rhb_ts))+1;
% Umin = min(min(sfa_u),min(rhb_u))-2; Umax = max(max(sfa_u),max(rhb_u))+2;
% subplot(2,2,1);
%     plot(sfa_qa,rhb_qa,'bo',one2one,one2one,'k-'); grid;
%     axis([qamin,qamax,qamin,qamax]); xlabel('SURFA'); ylabel('RHB');
%     title('Specific Humidity, g/kg');
% subplot(2,2,2);
%     plot(sfa_ta,rhb_ta,'bo',one2one,one2one,'k-'); grid;
%     axis([tamin,tamax,tamin,tamax]); xlabel('SURFA'); ylabel('RHB');
%     title('Air Temp, C');
% subplot(2,2,3);
%     plot(sfa_ts,rhb_ts,'bo',one2one,one2one,'k-'); grid;
%     axis([tsmin,tsmax,tsmin,tsmax]); xlabel('SURFA'); ylabel('RHB');
%     title('Surface Temp, C');
% subplot(2,2,4);
%     plot(sfa_u,rhb_u,'bo',one2one,one2one,'k-'); grid;
%     axis([Umin,Umax,Umin,Umax]); xlabel('SURFA'); ylabel('RHB');
%     title('Wind Speed, m/s');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_met_comparison',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% % stress comparison
% figure();
% taumin = min(min(sfa_tau),min(rhb_taub))-0.2; taumax = max(max(sfa_tau),max(rhb_taub))+0.2;
% plot(sfa_tau,rhb_tauc,'bo',sfa_tau,rhb_tauib,'gd',sfa_tau,rhb_taub,'r*',one2one,one2one,'k-'); grid;
% axis([taumin,taumax,taumin,taumax]); xlabel('SURFA'); ylabel('RHB');
% legend('tauc','tauib','taub','location','northwest'); title('Wind Stress, N/m^2');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_stress_comparison',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% % latent heat comparison
% figure();
% hlmin = min(min(sfa_hl),min(rhb_hlb))-5; hlmax = max(max(sfa_hl),max(rhb_hlb))+5;
% plot(sfa_hl,rhb_hlc,'bo',sfa_hl,rhb_hlib,'gd',sfa_hl,rhb_hlb,'r*',one2one,one2one,'k-'); grid;
% axis([hlmin,hlmax,hlmin,hlmax]); xlabel('SURFA'); ylabel('RHB');
% legend('hlc','hlib','hlb','location','northwest'); title('Latent Heat Flux, W/m^2');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_hl_comparison',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% % sensible heat comparison
% figure();
% hsmin = min(min(sfa_hs),min(rhb_hsb))-5; hsmax = max(max(sfa_hs),max(rhb_hsb))+5;
% plot(sfa_hs,rhb_hsc,'bo',sfa_hs,rhb_hsib,'gd',sfa_hs,rhb_hsb,'r*',one2one,one2one,'k-'); grid;
% axis([hsmin,hsmax,hsmin,hsmax]); xlabel('SURFA'); ylabel('RHB');
% legend('hsc','hsib','hsb','location','northwest'); title('Sensible Heat Flux, W/m^2');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_hs_comparison',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% % downwelling solar comparison
% figure();
% rsmin = min(min(sfa_rs),min(rhb_rs))-20; rsmax = max(max(sfa_rs),max(rhb_rs))+20;
% plot(sfa_rs,rhb_rs,'bo',one2one,one2one,'k-'); grid;
% axis([rsmin,rsmax,rsmin,rsmax]); xlabel('SURFA'); ylabel('RHB');
% title('Downwelling Solar, W/m^2');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_rs_comparison',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% % downwelling IR comparison
% figure();
% rlmin = min(min(sfa_rl),min(rhb_rl))-10; rlmax = max(max(sfa_rl),max(rhb_rl))+10;
% plot(sfa_rl,rhb_rl,'bo',one2one,one2one,'k-'); grid;
% axis([rlmin,rlmax,rlmin,rlmax]); xlabel('SURFA'); ylabel('RHB');
% title('Downwelling IR, W/m^2');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_rl_comparison',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% close all;
% 
% %% daily time series plots
% 
% % met time series: qa, ta, ts, U
% figure;
% subplot(4,1,1);
%     plot(sfa.jd,rhb_qa,'ro',sfa.jd,sfa_qa,'k-'); grid;
%     axis([sfa.jd(1),sfa.jd(end),qamin,qamax]); ylabel('q air, g/kg');
%     legend('RHB','SURFA','location','northwest'); title('Specific Humidity');
% subplot(4,1,2);
%     plot(sfa.jd,rhb_ta,'ro',sfa.jd,sfa_ta,'k-'); grid;
%     axis([sfa.jd(1),sfa.jd(end),tamin,tamax]); ylabel('T air, C');
%     title('Sfc Air Temperature')
% subplot(4,1,3);
%     plot(sfa.jd,rhb_ts,'ro',sfa.jd,sfa_ts,'k-'); grid;
%     axis([sfa.jd(1),sfa.jd(end),tsmin,tsmax]); ylabel('SST, C');
%     title('Ocean Surface Temperature');
% subplot(4,1,4);
%     plot(sfa.jd,rhb_u,'ro',sfa.jd,sfa_u,'k-'); grid;
%     axis([sfa.jd(1),sfa.jd(end),Umin,Umax]); xlabel('DOY, 2015'); ylabel('U, m/s');
%     title('Wind Speed');
% orient tall;
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_met_timeseries',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% % stress timeseries
% figure();
% plot(sfa.jd,rhb_tauc,'bo',sfa.jd,rhb_tauib,'gd',sfa.jd,rhb_taub,'r*',sfa.jd,sfa_tau,'k-'); grid;
% axis([sfa.jd(1),sfa.jd(end),taumin,taumax]); xlabel('DOY, 2015'); ylabel('\tau, N/m^2');
% legend('tauc','tauib','taub','SURFA','location','northwest'); title('Wind Stress');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_stress_timeseries',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% % latent heat timeseries
% figure();
% plot(sfa.jd,rhb_hlc,'bo',sfa.jd,rhb_hlib,'gd',sfa.jd,rhb_hlb,'r*',sfa.jd,sfa_hl,'k-'); grid;
% axis([sfa.jd(1),sfa.jd(end),hlmin,hlmax]); xlabel('DOY, 2015'); ylabel('hl, W/m^2');
% legend('hlc','hlib','hlb','SURFA','location','northwest'); title('Latent Heat Flux');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_hl_timeseries',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% % sensible heat timeseries
% figure();
% plot(sfa.jd,rhb_hsc,'bo',sfa.jd,rhb_hsib,'gd',sfa.jd,rhb_hsb,'r*',sfa.jd,sfa_hs,'k-'); grid;
% axis([sfa.jd(1),sfa.jd(end),hsmin,hsmax]); xlabel('DOY, 2015'); ylabel('hs, W/m^2');
% legend('hsc','hsib','hsb','SURFA','location','northwest'); title('Sensible Heat Flux');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_hs_timeseries',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% % downwelling solar timeseries
% figure();
% plot(sfa.jd,rhb_rs,'bo',sfa.jd,sfa_rs,'k-'); grid;
% axis([sfa.jd(1),sfa.jd(end),rsmin,rsmax]); xlabel('DOY, 2015'); ylabel('rs, W/m^2');
% legend('RHB','SURFA','location','northwest'); title('Downwelling Solar Heat Flux');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_rs_timeseries',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% % downwelling IR timeseries
% figure();
% plot(sfa.jd,rhb_rl,'bo',sfa.jd,sfa_rl,'k-'); grid;
% axis([sfa.jd(1),sfa.jd(end),rlmin,rlmax]); xlabel('DOY, 2015'); ylabel('rl, W/m^2');
% legend('RHB','SURFA','location','northwest'); title('Downwelling IR Heat Flux');
% if save_plots
%     ppath = fullfile(surfa_fig_path,['CalWater_2015_SURFA_RHB_rl_timeseries',graphformat]);
%     print(graphdevice,ppath);
% end;
% 
% close all

%% movie plots
if run_movies
    x_lim = [sfa.lon(1),sfa.lon(end)];
    y_lim = [sfa.lat(1),sfa.lat(end)];
    N = length(sfa.jd);
    
    % evaporation
    mfl = fullfile(surfa_fig_path,'SURFA_evap.mp4');
    v_ev = VideoWriter(mfl,'MPEG-4');
    v_ev.FrameRate = 2;
    open(v_ev);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3); % date/time string
        imagesc(x_lim,y_lim,sfa.evp(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Evaporation, g/m^2/s';
        caxis([-3e-4,1e-4]); c.Ticks = -3e-4:0.5e-4:1e-4;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Evaporation ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_ev,frame);
    end
    close(v_ev);
    
    % low cloud percent
    mfl = fullfile(surfa_fig_path,'SURFA_lowcld.mp4');
    v_lcld = VideoWriter(mfl,'MPEG-4');
    v_lcld.FrameRate = 2;
    open(v_lcld);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.lowcld(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Low Cloud Fraction, %';
        caxis([0,100]); c.Ticks = 0:20:100;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Low Cloud ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_lcld,frame);
    end
    close(v_lcld);
    
    % mid cloud percent
    mfl = fullfile(surfa_fig_path,'SURFA_midcld.mp4');
    v_mcld = VideoWriter(mfl,'MPEG-4');
    v_mcld.FrameRate = 2;
    open(v_mcld);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.midcld(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Mid Cloud Fraction, %';
        caxis([0,100]); c.Ticks = 0:20:100;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Mid Cloud ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_mcld,frame);
    end
    close(v_mcld);
    
    % cloud water
    mfl = fullfile(surfa_fig_path,'SURFA_wcloud.mp4');
    v_wcld = VideoWriter(mfl,'MPEG-4');
    v_wcld.FrameRate = 2;
    open(v_wcld);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.wcloud(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Cloud Water, g/m^2';
        caxis([0,3000]); c.Ticks = 0:500:3000;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Cloud Water ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%        hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_wcld,frame);
    end
    close(v_wcld);
    
    % pressure
    mfl = fullfile(surfa_fig_path,'SURFA_sfcP.mp4');
    v_P = VideoWriter(mfl,'MPEG-4');
    v_P.FrameRate = 2;
    open(v_P);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.P(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Pressure, mb';
        caxis([960,1020]); c.Ticks = 960:20:1020;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Sfc Pressure ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_P,frame);
    end
    close(v_P);
    
    % surface temperature
    mfl = fullfile(surfa_fig_path,'SURFA_ta.mp4');
    v_ta = VideoWriter(mfl,'MPEG-4');
    v_ta.FrameRate = 2;
    open(v_ta);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.ta(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Sfc Air Temp, C';
        caxis([0,30]); c.Ticks = 0:10:30;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF 2m Air Temp ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_ta,frame);
    end
    close(v_ta);
    
    % specific humidity
    mfl = fullfile(surfa_fig_path,'SURFA_qa.mp4');
    v_qa = VideoWriter(mfl,'MPEG-4');
    v_qa.FrameRate = 2;
    open(v_qa);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.qa(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Specific Humidity, g/kg';
        caxis([5,20]); c.Ticks = 5:5:20;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Specific Humidity ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_qa,frame);
    end
    close(v_qa);
    
    % wind speed
    mfl = fullfile(surfa_fig_path,'SURFA_wspd.mp4');
    v_wspd = VideoWriter(mfl,'MPEG-4');
    v_wspd.FrameRate = 2;
    open(v_wspd);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.u(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Wind Speed, m/s';
        caxis([0,20]); c.Ticks = 0:5:20;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Wind ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_wspd,frame);
    end
    close(v_wspd);
    
    % wind stress
    mfl = fullfile(surfa_fig_path,'SURFA_tau.mp4');
    v_tau = VideoWriter(mfl,'MPEG-4');
    v_tau.FrameRate = 2;
    open(v_tau);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.tau(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Stress, N/m^2';
        caxis([0,1]); c.Ticks = 0:0.2:1;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Wind Stress ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_tau,frame);
    end
    close(v_tau);
    
    % sensible heat flux
    mfl = fullfile(surfa_fig_path,'SURFA_hs.mp4');
    v_hs = VideoWriter(mfl,'MPEG-4');
    v_hs.FrameRate = 2;
    open(v_hs);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.hs(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Sens Heat, W/m^2';
        caxis([-100,200]); c.Ticks = -100:50:200;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Sensible Heat Flux ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_hs,frame);
    end
    close(v_hs);
    
    % latent heat flux
    mfl = fullfile(surfa_fig_path,'SURFA_hl.mp4');
    v_hl = VideoWriter(mfl,'MPEG-4');
    v_hl.FrameRate = 2;
    open(v_hl);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.hl(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Latent Heat, W/m^2';
        caxis([-100,600]); c.Ticks = -100:100:600;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Latent Heat Flux ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_hl,frame);
    end
    close(v_hl);
    
    % downwelling solar flux
    mfl = fullfile(surfa_fig_path,'SURFA_rs.mp4');
    v_rs = VideoWriter(mfl,'MPEG-4');
    v_rs.FrameRate = 2;
    open(v_rs);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.rs(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Sfc Solar Flux, W/m^2';
        caxis([0,1000]); c.Ticks = 0:200:1000;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Sfc Solar ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_rs,frame);
    end
    close(v_rs);
    
    % downwelling IR flux
    mfl = fullfile(surfa_fig_path,'SURFA_rl.mp4');
    v_rl = VideoWriter(mfl,'MPEG-4');
    v_rl.FrameRate = 2;
    open(v_rl);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.rl(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Sfc IR Flux, W/m^2';
        caxis([0,500]); c.Ticks = 0:100:500;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF Sfc IR ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_rl,frame);
    end
    close(v_rl);
    
    % SST / Soil Temperature
    mfl = fullfile(surfa_fig_path,'SURFA_ts.mp4');
    v_ts = VideoWriter(mfl,'MPEG-4');
    v_ts.FrameRate = 2;
    open(v_ts);
    figure(1);
    for ii = 1:N
        dstr = sfa.dfl{ii}(end-15:end-3);
        imagesc(x_lim,y_lim,sfa.ts(:,:,ii)'); axis xy; if ii==1; grid; end
        c = colorbar; c.Label.String = 'Sfc Temp, C';
        caxis([0,30]); c.Ticks = 0:10:30;
        ylabel('latitiude'); xlabel('longitude');
        title(['DOY ',sprintf('%07.3f',sfa.jd(ii)),'  SURFA-ECMWF SST ',dstr(1:4),' ',dstr(5:6),...
            ' ',dstr(7:8),' ',dstr(9:10),'+',dstr(12:13)]);
%         hold on; plot(rhb_lon(ii),rhb_lat(ii),'mo');
        drawnow
        frame = getframe(gcf);
        writeVideo(v_ts,frame);
    end
    close(v_ts);
    close (1);
end
