POP_heat_transport_basin_plots.m

% plot MHT as a function of latitude in each basin

path(path,'~/POP/')
path(path,'~/plotting_scripts/')
cd('/indopac/adelman/Global_mesoscale/POP/')

lat_transects_vec_South = ((-70):1:(-35))';
lat_transects_vec = [-34.7; ((-34):1:70)'];
depth_range = [0 6000];

year_offset = 1976;
spatial_scheme_file_id = 'errfcn_altbuffer_10_20';    % indicator of spatial scheme in file name
time_scale_separation = (3/12)*365;   % time scale to separate low frequency from high frequency variability

time_range_start = [1978 01 01];
time_range_end = [2010 01 01];

grid_res = 0.1;

load('Multitransect_arrays_POP_Atl_1deg.mat','delta_t','datenum_start','datenum_end','time_datenum_in_range')

datestr_title_start = [num2str(floor(datenum_start/365) + year_offset),'-',datestr([1990 01 (mod(datenum_start,365) + 1) 0 0 0],'mm-dd')];
datestr_title_end_minus1 = [num2str(floor((datenum_end - 1)/365) + year_offset),'-',datestr([1990 01 (mod(datenum_end - 1,365) + 1) 0 0 0],'mm-dd')];
datestr_start = [num2str(floor(datenum_start/365) + year_offset),datestr([1990 01 (mod(datenum_start,365) + 1) 0 0 0],'mmdd')];
datestr_end_minus1 = [num2str(floor((datenum_end - 1)/365) + year_offset),datestr([1990 01 (mod(datenum_end - 1,365) + 1) 0 0 0],'mmdd')];

n_bins = round(365/delta_t);


% plot heat transport time mean and standard deviation from each component

basin_names = {'Atl'; 'IndoPac'; 'South'};

for curr_basin_ind = 1:length(basin_names)
    curr_basin_name = basin_names{curr_basin_ind};
    
    load(['Multitransect_arrays_POP_',curr_basin_name,'_1deg.mat'],'*tseries_array')

    curr_array = heat_transport_zonmean_tmean_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_zonmean_tmean_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array = heat_transport_zonmean_lowfreq_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_zonmean_lowfreq_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_zonmean_lowfreq_tmean;
    heat_transport_zonmean_lowfreq_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_zonmean_lowfreq_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array = heat_transport_zonmean_highfreq_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_zonmean_highfreq_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_zonmean_highfreq_tmean;
    heat_transport_zonmean_highfreq_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_zonmean_highfreq_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array = heat_transport_zonmean_lowfreq_tseries_array + heat_transport_zonmean_highfreq_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_zonmean_lowphigh_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_zonmean_lowphigh_tmean;
    heat_transport_zonmean_lowphigh_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_zonmean_lowphigh_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array = heat_transport_largescale_tmean_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_largescale_tmean_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array = heat_transport_largescale_lowfreq_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_largescale_lowfreq_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_largescale_lowfreq_tmean;
    heat_transport_largescale_lowfreq_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_largescale_lowfreq_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array = heat_transport_largescale_highfreq_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_largescale_highfreq_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_largescale_highfreq_tmean;
    heat_transport_largescale_highfreq_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_largescale_highfreq_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array = heat_transport_largescale_lowfreq_tseries_array + heat_transport_largescale_highfreq_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_largescale_lowphigh_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_largescale_lowphigh_tmean;
    heat_transport_largescale_lowphigh_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_largescale_lowphigh_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array = heat_transport_mesoscale_tmean_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_mesoscale_tmean_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array = heat_transport_mesoscale_lowfreq_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_mesoscale_lowfreq_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_mesoscale_lowfreq_tmean;
    heat_transport_mesoscale_lowfreq_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_mesoscale_lowfreq_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array = heat_transport_mesoscale_highfreq_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_mesoscale_highfreq_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_mesoscale_highfreq_tmean;
    heat_transport_mesoscale_highfreq_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_mesoscale_highfreq_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array = heat_transport_mesoscale_lowfreq_tseries_array + heat_transport_mesoscale_highfreq_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_mesoscale_lowphigh_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_mesoscale_lowphigh_tmean;
    heat_transport_mesoscale_lowphigh_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_mesoscale_lowphigh_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);

    curr_array = heat_transport_all_tseries_array - (heat_transport_zonmean_tmean_tseries_array + heat_transport_zonmean_lowfreq_tseries_array + heat_transport_zonmean_highfreq_tseries_array + heat_transport_largescale_tmean_tseries_array + heat_transport_largescale_lowfreq_tseries_array + heat_transport_largescale_highfreq_tseries_array + heat_transport_mesoscale_tmean_tseries_array + heat_transport_mesoscale_lowfreq_tseries_array + heat_transport_mesoscale_highfreq_tseries_array);
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_residual_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_residual_tmean;
    heat_transport_residual_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_residual_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);

    curr_array = heat_transport_vT_timevar_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_vT_timevar_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_vT_timevar_tmean;
    heat_transport_vT_timevar_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_vT_timevar_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);

    curr_array = heat_transport_zonmean_highfreq_tseries_array + heat_transport_largescale_highfreq_tseries_array + heat_transport_mesoscale_highfreq_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_vT_highfreq_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_vT_highfreq_tmean;
    heat_transport_vT_highfreq_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_vT_highfreq_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    

    curr_array = heat_transport_all_tseries_array;
    curr_nan_mask = ones(size(curr_array));
    curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
    curr_array(isnan(curr_array) == 1) = 0;
    heat_transport_all_tmean = sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2));
    curr_array_mean = heat_transport_all_tmean;
    heat_transport_all_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    curr_array_noseason = NaN(size(curr_array));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_array_noseason(:,curr_in_bin_ind) = curr_array(:,curr_in_bin_ind) - repmat(sum(curr_nan_mask(:,curr_in_bin_ind).*curr_array(:,curr_in_bin_ind),2)./(sum(curr_nan_mask(:,curr_in_bin_ind),2)),[1 length(curr_in_bin_ind)]);
    end
    curr_array = curr_array_noseason + repmat(curr_array_mean,[1 size(curr_array,2)]);
    curr_array_ID = bandpass_err_fcn(curr_array - repmat(curr_array_mean,[1 size(curr_array,2)]),2,delta_t,1/(3*sum(diff(time_datenum_in_range))),1/426,5,1,1,1,0);
    curr_nan_mask((isnan(curr_array_ID) == 1) | (abs(curr_array_ID) < 1e-10)) = 0;
    curr_array_ID(isnan(curr_array_ID) == 1) = 0;
    curr_array_mean = sum(curr_nan_mask.*curr_array_ID,2)./(sum(curr_nan_mask,2));
    heat_transport_all_stddev_ID = (sum(curr_nan_mask.*((curr_array_ID - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2))).^(1/2);
    
    eval(['heat_transport_zonmean_tseries_array_',curr_basin_name,' = heat_transport_zonmean_tmean_tseries_array + heat_transport_zonmean_lowfreq_tseries_array + heat_transport_zonmean_highfreq_tseries_array;'])
    eval(['heat_transport_largescale_tseries_array_',curr_basin_name,' = heat_transport_largescale_tmean_tseries_array + heat_transport_largescale_lowfreq_tseries_array + heat_transport_largescale_highfreq_tseries_array;'])
    eval(['heat_transport_mesoscale_tseries_array_',curr_basin_name,' = heat_transport_mesoscale_tmean_tseries_array + heat_transport_mesoscale_lowfreq_tseries_array + heat_transport_mesoscale_highfreq_tseries_array;'])
    eval(['heat_transport_residual_tseries_array_',curr_basin_name,' = heat_transport_all_tseries_array - (heat_transport_zonmean_tmean_tseries_array + heat_transport_zonmean_lowfreq_tseries_array + heat_transport_zonmean_highfreq_tseries_array + heat_transport_largescale_tmean_tseries_array + heat_transport_largescale_lowfreq_tseries_array + heat_transport_largescale_highfreq_tseries_array + heat_transport_mesoscale_tmean_tseries_array + heat_transport_mesoscale_lowfreq_tseries_array + heat_transport_mesoscale_highfreq_tseries_array);'])
    eval(['heat_transport_vT_timevar_tseries_array_',curr_basin_name,' = heat_transport_vT_timevar_tseries_array;'])
    eval(['heat_transport_vT_highfreq_tseries_array_',curr_basin_name,' = heat_transport_zonmean_highfreq_tseries_array + heat_transport_largescale_highfreq_tseries_array + heat_transport_mesoscale_highfreq_tseries_array;'])
    eval(['heat_transport_all_tseries_array_',curr_basin_name,' = heat_transport_all_tseries_array;'])

    eval(['heat_transport_zonmean_tmean_tmean_',curr_basin_name,' = heat_transport_zonmean_tmean_tmean;'])
    eval(['heat_transport_zonmean_lowfreq_tmean_',curr_basin_name,' = heat_transport_zonmean_lowfreq_tmean;'])
    eval(['heat_transport_zonmean_highfreq_tmean_',curr_basin_name,' = heat_transport_zonmean_highfreq_tmean;'])
    eval(['heat_transport_zonmean_lowphigh_stddev_ID_',curr_basin_name,' = heat_transport_zonmean_lowphigh_stddev_ID;'])
    eval(['heat_transport_largescale_tmean_tmean_',curr_basin_name,' = heat_transport_largescale_tmean_tmean;'])
    eval(['heat_transport_largescale_lowfreq_tmean_',curr_basin_name,' = heat_transport_largescale_lowfreq_tmean;'])
    eval(['heat_transport_largescale_highfreq_tmean_',curr_basin_name,' = heat_transport_largescale_highfreq_tmean;'])
    eval(['heat_transport_largescale_lowphigh_stddev_ID_',curr_basin_name,' = heat_transport_largescale_lowphigh_stddev_ID;'])
    eval(['heat_transport_mesoscale_tmean_tmean_',curr_basin_name,' = heat_transport_mesoscale_tmean_tmean;'])
    eval(['heat_transport_mesoscale_lowfreq_tmean_',curr_basin_name,' = heat_transport_mesoscale_lowfreq_tmean;'])
    eval(['heat_transport_mesoscale_highfreq_tmean_',curr_basin_name,' = heat_transport_mesoscale_highfreq_tmean;'])
    eval(['heat_transport_mesoscale_lowphigh_stddev_ID_',curr_basin_name,' = heat_transport_mesoscale_lowphigh_stddev_ID;'])
    eval(['heat_transport_residual_tmean_',curr_basin_name,' = heat_transport_residual_tmean;'])
    eval(['heat_transport_residual_stddev_ID_',curr_basin_name,' = heat_transport_residual_stddev_ID;'])
    eval(['heat_transport_vT_timevar_tmean_',curr_basin_name,' = heat_transport_vT_timevar_tmean;'])
    eval(['heat_transport_vT_timevar_stddev_ID_',curr_basin_name,' = heat_transport_vT_timevar_stddev_ID;'])
    eval(['heat_transport_vT_highfreq_tmean_',curr_basin_name,' = heat_transport_vT_highfreq_tmean;'])
    eval(['heat_transport_vT_highfreq_stddev_ID_',curr_basin_name,' = heat_transport_vT_highfreq_stddev_ID;'])
    eval(['heat_transport_all_tmean_',curr_basin_name,' = heat_transport_all_tmean;'])
    eval(['heat_transport_all_stddev_ID_',curr_basin_name,' = heat_transport_all_stddev_ID;'])
end


aspect_ratio_baseline = 10;
xtick_spacing = 10;


fig17 = figure(17);
close(figure(1))
h = plot(lat_transects_vec,(1e-15)*(heat_transport_zonmean_tmean_tmean_Atl + heat_transport_zonmean_lowfreq_tmean_Atl + heat_transport_zonmean_highfreq_tmean_Atl),lat_transects_vec,(1e-15)*(heat_transport_largescale_tmean_tmean_Atl + heat_transport_largescale_lowfreq_tmean_Atl + heat_transport_largescale_highfreq_tmean_Atl),lat_transects_vec,(1e-15)*(heat_transport_mesoscale_tmean_tmean_Atl + heat_transport_mesoscale_lowfreq_tmean_Atl + heat_transport_mesoscale_highfreq_tmean_Atl),lat_transects_vec,(1e-15)*heat_transport_residual_tmean_Atl,lat_transects_vec,(1e-15)*heat_transport_all_tmean_Atl);
leg = legend('Overturning','Large scale','Mesoscale','Residual','Total','location','southeast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec) max(lat_transects_vec)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[-1 1],'ytick',(-1):0.2:1)
daspect([(((15/7)*aspect_ratio_baseline) + ((1.5/7)*(max(lat_transects_vec) - min(lat_transects_vec)))) 1 1])
set(h(1),'Color',[0.55 0 0.55],'LineWidth',2)
set(h(2),'Color',[0 0 0.8],'LineWidth',2)
set(h(3),'Color',[0.8 0 0],'LineWidth',2)
set(h(4),'Color',[0.8 0.6 0.4],'LineWidth',2)
set(h(5),'Color',[0 0 0],'LineWidth',2)
hold on
line([min(lat_transects_vec) max(lat_transects_vec)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('Heat transport (PW)')
title({'POP total time-averaged heat transport from spatial components,'; ['across Atlantic basin, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig17,['Heat_transport_spatial_components_POP_Atl_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_tavg.pdf'])
close(fig17)


fig17 = figure(17);
close(figure(1))
h = plot(lat_transects_vec,(1e-15)*heat_transport_zonmean_lowphigh_stddev_ID_Atl,lat_transects_vec,(1e-15)*heat_transport_largescale_lowphigh_stddev_ID_Atl,lat_transects_vec,(1e-15)*heat_transport_mesoscale_lowphigh_stddev_ID_Atl,lat_transects_vec,(1e-15)*heat_transport_residual_stddev_ID_Atl,lat_transects_vec,(1e-15)*heat_transport_all_stddev_ID_Atl);
leg = legend('Overturning','Large scale','Mesoscale','Residual','Total','location','NorthEast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec) max(lat_transects_vec)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[0 0.2],'ytick',0:0.05:0.3)
daspect([(((120/7)*aspect_ratio_baseline) + ((12/7)*(max(lat_transects_vec) - min(lat_transects_vec)))) 1 1])
set(h(1),'Color',[0.55 0 0.55],'LineWidth',2)
set(h(2),'Color',[0 0 0.8],'LineWidth',2)
set(h(3),'Color',[0.8 0 0],'LineWidth',2)
set(h(4),'Color',[0.8 0.6 0.4],'LineWidth',2)
set(h(5),'Color',[0 0 0],'LineWidth',2)
hold on
line([min(lat_transects_vec) max(lat_transects_vec)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('ID heat transport std. dev. (PW)')
title({'POP ID standard deviation of heat transport from spatial components,'; ['across Atlantic basin, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig17,['Heat_transport_spatial_components_POP_Atl_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_stddev_ID.pdf'])
close(fig17)


fig18 = figure(18);
close(figure(1))
h = plot(lat_transects_vec,(1e-15)*(heat_transport_mesoscale_tmean_tmean_Atl + heat_transport_mesoscale_lowfreq_tmean_Atl + heat_transport_mesoscale_highfreq_tmean_Atl),lat_transects_vec,(1e-15)*heat_transport_vT_timevar_tmean_Atl,lat_transects_vec,(1e-15)*heat_transport_vT_highfreq_tmean_Atl);
leg = legend('Mesoscale','Time-varying','High freq.','location','SouthEast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec) max(lat_transects_vec)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[-0.4 0.4],'ytick',(-0.4):0.1:0.4)
daspect([(((45/7)*aspect_ratio_baseline) + ((4.5/7)*(max(lat_transects_vec) - min(lat_transects_vec)))) 1 1])
set(h(1),'Color',[0.8 0 0],'LineWidth',2)
set(h(2),'Color',[0.4 0.3 0.7],'LineWidth',2)
set(h(3),'Color',[0.8 0 0.8],'LineWidth',2,'LineStyle','-.')
hold on
line([min(lat_transects_vec) max(lat_transects_vec)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('Heat transport (PW)')
title({'POP total time-averaged heat transport from various "eddy" formulations,'; ['across Atlantic basin, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig18,['Heat_transport_eddy_formulations_POP_Atl_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_tavg.pdf'])
close(fig18)


fig18 = figure(18);
close(figure(1))
h = plot(lat_transects_vec,(1e-15)*heat_transport_mesoscale_lowphigh_stddev_ID_Atl,lat_transects_vec,(1e-15)*heat_transport_vT_timevar_stddev_ID_Atl,lat_transects_vec,(1e-15)*heat_transport_vT_highfreq_stddev_ID_Atl);
leg = legend('Mesoscale','Time-varying','High freq.','location','NorthEast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec) max(lat_transects_vec)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[0 0.2],'ytick',0:0.05:0.3)
daspect([(((180/7)*aspect_ratio_baseline) + ((18/7)*(max(lat_transects_vec) - min(lat_transects_vec)))) 1 1])
set(h(1),'Color',[0.8 0 0],'LineWidth',2)
set(h(2),'Color',[0.4 0.3 0.7],'LineWidth',2)
set(h(3),'Color',[0.8 0 0.8],'LineWidth',2,'LineStyle','-.')
hold on
line([min(lat_transects_vec) max(lat_transects_vec)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('ID heat transport std. dev. (PW)')
title({'POP ID std. dev. of heat transport from various "eddy" formulations,'; ['across Atlantic basin, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)
keyboard
saveas(fig18,['Heat_transport_eddy_formulations_POP_Atl_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_stddev_ID.pdf'])
close(fig18)



fig17 = figure(17);
close(figure(1))
h = plot(lat_transects_vec,(1e-15)*(heat_transport_zonmean_tmean_tmean_IndoPac + heat_transport_zonmean_lowfreq_tmean_IndoPac + heat_transport_zonmean_highfreq_tmean_IndoPac),lat_transects_vec,(1e-15)*(heat_transport_largescale_tmean_tmean_IndoPac + heat_transport_largescale_lowfreq_tmean_IndoPac + heat_transport_largescale_highfreq_tmean_IndoPac),lat_transects_vec,(1e-15)*(heat_transport_mesoscale_tmean_tmean_IndoPac + heat_transport_mesoscale_lowfreq_tmean_IndoPac + heat_transport_mesoscale_highfreq_tmean_IndoPac),lat_transects_vec,(1e-15)*heat_transport_residual_tmean_IndoPac,lat_transects_vec,(1e-15)*heat_transport_all_tmean_IndoPac);
leg = legend('Overturning','Large scale','Mesoscale','Residual','Total','location','southeast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec) max(lat_transects_vec)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[-2 1.4],'ytick',(-2):0.2:2)
daspect([(((10/7)*aspect_ratio_baseline) + ((1/7)*(max(lat_transects_vec) - min(lat_transects_vec)))) 1 1])
set(h(1),'Color',[0.55 0 0.55],'LineWidth',2)
set(h(2),'Color',[0 0 0.8],'LineWidth',2)
set(h(3),'Color',[0.8 0 0],'LineWidth',2)
set(h(4),'Color',[0.8 0.6 0.4],'LineWidth',2)
set(h(5),'Color',[0 0 0],'LineWidth',2)
hold on
line([min(lat_transects_vec) max(lat_transects_vec)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('Heat transport (PW)')
title({'POP total time-averaged heat transport from various "eddy" formulations,'; ['across Indo-Pacific basin, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig17,['Heat_transport_spatial_components_POP_IndoPac_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_tavg.pdf'])
close(fig17)


fig17 = figure(17);
close(figure(1))
h = plot(lat_transects_vec,(1e-15)*heat_transport_zonmean_lowphigh_stddev_ID_IndoPac,lat_transects_vec,(1e-15)*heat_transport_largescale_lowphigh_stddev_ID_IndoPac,lat_transects_vec,(1e-15)*heat_transport_mesoscale_lowphigh_stddev_ID_IndoPac,lat_transects_vec,(1e-15)*heat_transport_residual_stddev_ID_IndoPac,lat_transects_vec,(1e-15)*heat_transport_all_stddev_ID_IndoPac);
leg = legend('Overturning','Large scale','Mesoscale','Residual','Total','location','NorthEast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec) max(lat_transects_vec)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[0 0.9],'ytick',0:0.1:0.9)
daspect([(((40/7)*aspect_ratio_baseline) + ((4/7)*(max(lat_transects_vec) - min(lat_transects_vec)))) 1 1])
set(h(1),'Color',[0.55 0 0.55],'LineWidth',2)
set(h(2),'Color',[0 0 0.8],'LineWidth',2)
set(h(3),'Color',[0.8 0 0],'LineWidth',2)
set(h(4),'Color',[0.8 0.6 0.4],'LineWidth',2)
set(h(5),'Color',[0 0 0],'LineWidth',2)
hold on
line([min(lat_transects_vec) max(lat_transects_vec)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('ID heat transport std. dev. (PW)')
title({'POP ID std. dev. of heat transport from various "eddy" formulations,'; ['across Indo-Pacific basin, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig17,['Heat_transport_spatial_components_POP_IndoPac_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_stddev_ID.pdf'])
close(fig17)


fig18 = figure(18);
close(figure(1))
h = plot(lat_transects_vec,(1e-15)*(heat_transport_mesoscale_tmean_tmean_IndoPac + heat_transport_mesoscale_lowfreq_tmean_IndoPac + heat_transport_mesoscale_highfreq_tmean_IndoPac),lat_transects_vec,(1e-15)*heat_transport_vT_timevar_tmean_IndoPac,lat_transects_vec,(1e-15)*heat_transport_vT_highfreq_tmean_IndoPac);
leg = legend('Mesoscale','Time-varying','High freq.','location','SouthEast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec) max(lat_transects_vec)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[-0.8 0.8],'ytick',(-0.8):0.2:0.8)
daspect([(((22.5/7)*aspect_ratio_baseline) + ((2.25/7)*(max(lat_transects_vec) - min(lat_transects_vec)))) 1 1])
set(h(1),'Color',[0.8 0 0],'LineWidth',2)
set(h(2),'Color',[0.4 0.3 0.7],'LineWidth',2)
set(h(3),'Color',[0.8 0 0.8],'LineWidth',2,'LineStyle','-.')
hold on
line([min(lat_transects_vec) max(lat_transects_vec)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('Heat transport (PW)')
title({'POP total time-averaged heat transport from spatial components,'; ['across Indo-Pacific basin, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig18,['Heat_transport_eddy_formulations_POP_IndoPac_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_tavg.pdf'])
close(fig18)


fig18 = figure(18);
close(figure(1))
h = plot(lat_transects_vec,(1e-15)*heat_transport_mesoscale_lowphigh_stddev_ID_IndoPac,lat_transects_vec,(1e-15)*heat_transport_vT_timevar_stddev_ID_IndoPac,lat_transects_vec,(1e-15)*heat_transport_vT_highfreq_stddev_ID_IndoPac);
leg = legend('Mesoscale','Time-varying','High freq.','location','NorthEast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec) max(lat_transects_vec)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[0 0.2],'ytick',0:0.05:0.3)
daspect([(((180/7)*aspect_ratio_baseline) + ((18/7)*(max(lat_transects_vec) - min(lat_transects_vec)))) 1 1])
set(h(1),'Color',[0.8 0 0],'LineWidth',2)
set(h(2),'Color',[0.4 0.3 0.7],'LineWidth',2)
set(h(3),'Color',[0.8 0 0.8],'LineWidth',2,'LineStyle','-.')
hold on
line([min(lat_transects_vec) max(lat_transects_vec)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('ID heat transport std. dev. (PW)')
title({'POP ID standard deviation of heat transport from spatial components,'; ['across Indo-Pacific basin, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig18,['Heat_transport_eddy_formulations_POP_IndoPac_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_stddev_ID.pdf'])
close(fig18)




% produce global plots as a function of latitude

lat_transects_vec_global = [lat_transects_vec_South; lat_transects_vec];
heat_transport_zonmean_tmean_tmean_global = [heat_transport_zonmean_tmean_tmean_South; (heat_transport_zonmean_tmean_tmean_Atl + heat_transport_zonmean_tmean_tmean_IndoPac)];
heat_transport_zonmean_lowfreq_tmean_global = [heat_transport_zonmean_lowfreq_tmean_South; (heat_transport_zonmean_lowfreq_tmean_Atl + heat_transport_zonmean_lowfreq_tmean_IndoPac)];
heat_transport_zonmean_highfreq_tmean_global = [heat_transport_zonmean_highfreq_tmean_South; (heat_transport_zonmean_highfreq_tmean_Atl + heat_transport_zonmean_highfreq_tmean_IndoPac)];
heat_transport_largescale_tmean_tmean_global = [heat_transport_largescale_tmean_tmean_South; (heat_transport_largescale_tmean_tmean_Atl + heat_transport_largescale_tmean_tmean_IndoPac)];
heat_transport_largescale_lowfreq_tmean_global = [heat_transport_largescale_lowfreq_tmean_South; (heat_transport_largescale_lowfreq_tmean_Atl + heat_transport_largescale_lowfreq_tmean_IndoPac)];
heat_transport_largescale_highfreq_tmean_global = [heat_transport_largescale_highfreq_tmean_South; (heat_transport_largescale_highfreq_tmean_Atl + heat_transport_largescale_highfreq_tmean_IndoPac)];
heat_transport_mesoscale_tmean_tmean_global = [heat_transport_mesoscale_tmean_tmean_South; (heat_transport_mesoscale_tmean_tmean_Atl + heat_transport_mesoscale_tmean_tmean_IndoPac)];
heat_transport_mesoscale_lowfreq_tmean_global = [heat_transport_mesoscale_lowfreq_tmean_South; (heat_transport_mesoscale_lowfreq_tmean_Atl + heat_transport_mesoscale_lowfreq_tmean_IndoPac)];
heat_transport_mesoscale_highfreq_tmean_global = [heat_transport_mesoscale_highfreq_tmean_South; (heat_transport_mesoscale_highfreq_tmean_Atl + heat_transport_mesoscale_highfreq_tmean_IndoPac)];
heat_transport_residual_tmean_global = [heat_transport_residual_tmean_South; (heat_transport_residual_tmean_Atl + heat_transport_residual_tmean_IndoPac)];
heat_transport_vT_timevar_tmean_global = [heat_transport_vT_timevar_tmean_South; (heat_transport_vT_timevar_tmean_Atl + heat_transport_vT_timevar_tmean_IndoPac)];
heat_transport_vT_highfreq_tmean_global = [heat_transport_vT_highfreq_tmean_South; (heat_transport_vT_highfreq_tmean_Atl + heat_transport_vT_highfreq_tmean_IndoPac)];
heat_transport_all_tmean_global = [heat_transport_all_tmean_South; (heat_transport_all_tmean_Atl + heat_transport_all_tmean_IndoPac)];

heat_transport_zonmean_tseries_array_global = [heat_transport_zonmean_tseries_array_South; (heat_transport_zonmean_tseries_array_Atl + heat_transport_zonmean_tseries_array_IndoPac)];
heat_transport_largescale_tseries_array_global = [heat_transport_largescale_tseries_array_South; (heat_transport_largescale_tseries_array_Atl + heat_transport_largescale_tseries_array_IndoPac)];
heat_transport_mesoscale_tseries_array_global = [heat_transport_mesoscale_tseries_array_South; (heat_transport_mesoscale_tseries_array_Atl + heat_transport_mesoscale_tseries_array_IndoPac)];
heat_transport_residual_tseries_array_global = [heat_transport_residual_tseries_array_South; (heat_transport_residual_tseries_array_Atl + heat_transport_residual_tseries_array_IndoPac)];
heat_transport_vT_timevar_tseries_array_global = [heat_transport_vT_timevar_tseries_array_South; (heat_transport_vT_timevar_tseries_array_Atl + heat_transport_vT_timevar_tseries_array_IndoPac)];
heat_transport_vT_highfreq_tseries_array_global = [heat_transport_vT_highfreq_tseries_array_South; (heat_transport_vT_highfreq_tseries_array_Atl + heat_transport_vT_highfreq_tseries_array_IndoPac)];
heat_transport_all_tseries_array_global = [heat_transport_all_tseries_array_South; (heat_transport_all_tseries_array_Atl + heat_transport_all_tseries_array_IndoPac)];


% compute ID standard deviations globally

n_bins = round(365/delta_t);
heat_transport_zonmean_tseries_array_global_noseason = NaN(size(heat_transport_zonmean_tseries_array_global));
heat_transport_largescale_tseries_array_global_noseason = NaN(size(heat_transport_largescale_tseries_array_global));
heat_transport_mesoscale_tseries_array_global_noseason = NaN(size(heat_transport_mesoscale_tseries_array_global));
heat_transport_residual_tseries_array_global_noseason = NaN(size(heat_transport_residual_tseries_array_global));
heat_transport_vT_timevar_tseries_array_global_noseason = NaN(size(heat_transport_vT_timevar_tseries_array_global));
heat_transport_vT_highfreq_tseries_array_global_noseason = NaN(size(heat_transport_vT_highfreq_tseries_array_global));
heat_transport_all_tseries_array_global_noseason = NaN(size(heat_transport_all_tseries_array_global));
for bin_ind = 1:n_bins
    curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
    curr_nan_mask_1 = ones(size(heat_transport_largescale_tseries_array_global(:,curr_in_bin_ind)));
    curr_nan_mask_2 = ones(size(heat_transport_mesoscale_tseries_array_global(:,curr_in_bin_ind)));
    curr_nan_mask_1((isnan(heat_transport_largescale_tseries_array_global(:,curr_in_bin_ind)) == 1) | (abs(heat_transport_largescale_tseries_array_global(:,curr_in_bin_ind)) < 1e-5)) = 0;
    curr_nan_mask_2((isnan(heat_transport_mesoscale_tseries_array_global(:,curr_in_bin_ind)) == 1) | (abs(heat_transport_mesoscale_tseries_array_global(:,curr_in_bin_ind)) < 1e-5)) = 0;
    curr_nan_mask = curr_nan_mask_1 & curr_nan_mask_2;

    curr_array = heat_transport_zonmean_tseries_array_global(:,curr_in_bin_ind);
    curr_array(abs(curr_nan_mask) < 1e-5) = 0;
    heat_transport_zonmean_tseries_array_global_noseason(:,curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
    curr_array = heat_transport_largescale_tseries_array_global(:,curr_in_bin_ind);
    curr_array(abs(curr_nan_mask) < 1e-5) = 0;
    heat_transport_largescale_tseries_array_global_noseason(:,curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
    curr_array = heat_transport_mesoscale_tseries_array_global(:,curr_in_bin_ind);
    curr_array(abs(curr_nan_mask) < 1e-5) = 0;
    heat_transport_mesoscale_tseries_array_global_noseason(:,curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
    curr_array = heat_transport_residual_tseries_array_global(:,curr_in_bin_ind);
    curr_array(abs(curr_nan_mask) < 1e-5) = 0;
    heat_transport_residual_tseries_array_global_noseason(:,curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
    curr_array = heat_transport_vT_timevar_tseries_array_global(:,curr_in_bin_ind);
    curr_array(abs(curr_nan_mask) < 1e-5) = 0;
    heat_transport_vT_timevar_tseries_array_global_noseason(:,curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
    curr_array = heat_transport_vT_highfreq_tseries_array_global(:,curr_in_bin_ind);
    curr_array(abs(curr_nan_mask) < 1e-5) = 0;
    heat_transport_vT_highfreq_tseries_array_global_noseason(:,curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
    curr_array = heat_transport_all_tseries_array_global(:,curr_in_bin_ind);
    curr_array(abs(curr_nan_mask) < 1e-5) = 0;
    heat_transport_all_tseries_array_global_noseason(:,curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
end


% filter time series for interannual & decadal frequencies

half_power_adj = exp(erfinv((2^(1/2)) - 1)/5);   % adjustment factor to set bounds at half-power (rather than half-amplitude)
[heat_transport_zonmean_tseries_array_global_ID,~,~] = bandpass_err_fcn(heat_transport_zonmean_tseries_array_global_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
[heat_transport_largescale_tseries_array_global_ID,~,~] = bandpass_err_fcn(heat_transport_largescale_tseries_array_global_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
[heat_transport_mesoscale_tseries_array_global_ID,~,~] = bandpass_err_fcn(heat_transport_mesoscale_tseries_array_global_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
[heat_transport_residual_tseries_array_global_ID,~,~] = bandpass_err_fcn(heat_transport_residual_tseries_array_global_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
[heat_transport_vT_timevar_tseries_array_global_ID,~,~] = bandpass_err_fcn(heat_transport_vT_timevar_tseries_array_global_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
[heat_transport_vT_highfreq_tseries_array_global_ID,~,~] = bandpass_err_fcn(heat_transport_vT_highfreq_tseries_array_global_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
[heat_transport_all_tseries_array_global_ID,~,~] = bandpass_err_fcn(heat_transport_all_tseries_array_global_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
heat_transport_zonmean_tseries_array_global_ID = heat_transport_zonmean_tseries_array_global_ID + repmat(mean(heat_transport_zonmean_tseries_array_global,2),[1 length(time_datenum_in_range)]);
heat_transport_largescale_tseries_array_global_ID = heat_transport_largescale_tseries_array_global_ID + repmat(mean(heat_transport_largescale_tseries_array_global,2),[1 length(time_datenum_in_range)]);
heat_transport_mesoscale_tseries_array_global_ID = heat_transport_mesoscale_tseries_array_global_ID + repmat(mean(heat_transport_mesoscale_tseries_array_global,2),[1 length(time_datenum_in_range)]);
heat_transport_residual_tseries_array_global_ID = heat_transport_residual_tseries_array_global_ID + repmat(mean(heat_transport_residual_tseries_array_global,2),[1 length(time_datenum_in_range)]);
heat_transport_vT_timevar_tseries_array_global_ID = heat_transport_vT_timevar_tseries_array_global_ID + repmat(mean(heat_transport_vT_timevar_tseries_array_global,2),[1 length(time_datenum_in_range)]);
heat_transport_vT_highfreq_tseries_array_global_ID = heat_transport_vT_highfreq_tseries_array_global_ID + repmat(mean(heat_transport_vT_highfreq_tseries_array_global,2),[1 length(time_datenum_in_range)]);
heat_transport_all_tseries_array_global_ID = heat_transport_all_tseries_array_global_ID + repmat(mean(heat_transport_all_tseries_array_global,2),[1 length(time_datenum_in_range)]);


curr_array = heat_transport_zonmean_tseries_array_global_ID;
curr_nan_mask = ones(size(curr_array));
curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
curr_array(abs(curr_nan_mask) < 1e-5) = 0;
curr_array_mean = sum(curr_nan_mask.*curr_array,2)./sum(curr_nan_mask,2);
curr_array_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2) - 1)).^(1/2);
heat_transport_zonmean_global_stddev_ID = curr_array_stddev;
curr_array = heat_transport_largescale_tseries_array_global_ID;
curr_nan_mask = ones(size(curr_array));
curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
curr_array(abs(curr_nan_mask) < 1e-5) = 0;
curr_array_mean = sum(curr_nan_mask.*curr_array,2)./sum(curr_nan_mask,2);
curr_array_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2) - 1)).^(1/2);
heat_transport_largescale_global_stddev_ID = curr_array_stddev;
curr_array = heat_transport_mesoscale_tseries_array_global_ID;
curr_nan_mask = ones(size(curr_array));
curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
curr_array(abs(curr_nan_mask) < 1e-5) = 0;
curr_array_mean = sum(curr_nan_mask.*curr_array,2)./sum(curr_nan_mask,2);
curr_array_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2) - 1)).^(1/2);
heat_transport_mesoscale_global_stddev_ID = curr_array_stddev;
curr_array = heat_transport_residual_tseries_array_global_ID;
curr_nan_mask = ones(size(curr_array));
curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
curr_array(abs(curr_nan_mask) < 1e-5) = 0;
curr_array_mean = sum(curr_nan_mask.*curr_array,2)./sum(curr_nan_mask,2);
curr_array_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2) - 1)).^(1/2);
heat_transport_residual_global_stddev_ID = curr_array_stddev;
curr_array = heat_transport_vT_timevar_tseries_array_global_ID;
curr_nan_mask = ones(size(curr_array));
curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
curr_array(abs(curr_nan_mask) < 1e-5) = 0;
curr_array_mean = sum(curr_nan_mask.*curr_array,2)./sum(curr_nan_mask,2);
curr_array_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2) - 1)).^(1/2);
heat_transport_vT_timevar_global_tmean = curr_array_mean;
heat_transport_vT_timevar_global_stddev_ID = curr_array_stddev;
curr_array = heat_transport_vT_highfreq_tseries_array_global_ID;
curr_nan_mask = ones(size(curr_array));
curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
curr_array(abs(curr_nan_mask) < 1e-5) = 0;
curr_array_mean = sum(curr_nan_mask.*curr_array,2)./sum(curr_nan_mask,2);
curr_array_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2) - 1)).^(1/2);
heat_transport_vT_highfreq_global_tmean = curr_array_mean;
heat_transport_vT_highfreq_global_stddev_ID = curr_array_stddev;
curr_array = heat_transport_all_tseries_array_global_ID;
curr_nan_mask = ones(size(curr_array));
curr_nan_mask((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
curr_array(abs(curr_nan_mask) < 1e-5) = 0;
curr_array_mean = sum(curr_nan_mask.*curr_array,2)./sum(curr_nan_mask,2);
curr_array_stddev = (sum(curr_nan_mask.*((curr_array - repmat(curr_array_mean,[1 size(curr_array,2)])).^2),2)./(sum(curr_nan_mask,2) - 1)).^(1/2);
heat_transport_all_global_stddev_ID = curr_array_stddev;



fig17 = figure(17);
close(figure(1))
h = plot(lat_transects_vec_global,(1e-15)*(heat_transport_zonmean_tmean_tmean_global + heat_transport_zonmean_lowfreq_tmean_global + heat_transport_zonmean_highfreq_tmean_global),lat_transects_vec_global,(1e-15)*(heat_transport_largescale_tmean_tmean_global + heat_transport_largescale_lowfreq_tmean_global + heat_transport_largescale_highfreq_tmean_global),lat_transects_vec_global,(1e-15)*(heat_transport_mesoscale_tmean_tmean_global + heat_transport_mesoscale_lowfreq_tmean_global + heat_transport_mesoscale_highfreq_tmean_global),lat_transects_vec_global,(1e-15)*heat_transport_residual_tmean_global,lat_transects_vec_global,(1e-15)*heat_transport_all_tmean_global);
leg = legend('Overturning','Large scale','Mesoscale','Residual','Total','location','southeast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec_global) max(lat_transects_vec_global)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[-1.75 2.25],'ytick',(-2):0.25:2.25)
set(gca,'xticklabels',{''; ''; '-60'; ''; '-40'; ''; '-20'; ''; '0'; ''; '20'; ''; '40'; ''; '60'; ''})
set(gca,'yticklabels',{''; ''; '-1.5'; ''; '-1'; ''; '-0.5'; ''; '0'; ''; '0.5'; ''; '1'; ''; '1.5'; ''; '2'; ''})
daspect([(((10/7)*aspect_ratio_baseline) + ((1/7)*(max(lat_transects_vec_global) - min(lat_transects_vec_global)))) 1 1])
set(h(1),'Color',[0.55 0 0.55],'LineWidth',2)
set(h(2),'Color',[0 0 0.8],'LineWidth',2)
set(h(3),'Color',[0.8 0 0],'LineWidth',2)
set(h(4),'Color',[0.8 0.6 0.4],'LineWidth',2)
set(h(5),'Color',[0 0 0],'LineWidth',2)
hold on
line([min(lat_transects_vec_global) max(lat_transects_vec_global)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('Heat transport (PW)')
title({'POP total time-averaged heat transport from spatial components,'; ['globally integrated, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig17,['Heat_transport_spatial_components_POP_global_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_tavg.pdf'])
close(fig17)


fig17 = figure(17);
close(figure(1))
h = plot(lat_transects_vec_global,(1e-15)*heat_transport_zonmean_global_stddev_ID,lat_transects_vec_global,(1e-15)*heat_transport_largescale_global_stddev_ID,lat_transects_vec_global,(1e-15)*heat_transport_mesoscale_global_stddev_ID,lat_transects_vec_global,(1e-15)*heat_transport_residual_global_stddev_ID,lat_transects_vec_global,(1e-15)*heat_transport_all_global_stddev_ID);
leg = legend('Overturning','Large scale','Mesoscale','Residual','Total','location','NorthEast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec_global) max(lat_transects_vec_global)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[0 0.8],'ytick',0:0.05:0.8)
set(gca,'xticklabels',{''; ''; '-60'; ''; '-40'; ''; '-20'; ''; '0'; ''; '20'; ''; '40'; ''; '60'; ''})
set(gca,'yticklabels',{'0'; ''; '0.1'; ''; '0.2'; ''; '0.3'; ''; '0.4'; ''; '0.5'; ''; '0.6'; ''; '0.7'; ''; '0.8'})
daspect([(((40/7)*aspect_ratio_baseline) + ((4/7)*(max(lat_transects_vec_global) - min(lat_transects_vec_global)))) 1 1])
set(h(1),'Color',[0.55 0 0.55],'LineWidth',2)
set(h(2),'Color',[0 0 0.8],'LineWidth',2)
set(h(3),'Color',[0.8 0 0],'LineWidth',2)
set(h(4),'Color',[0.8 0.6 0.4],'LineWidth',2)
set(h(5),'Color',[0 0 0],'LineWidth',2)
hold on
line([min(lat_transects_vec_global) max(lat_transects_vec_global)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('ID heat transport std. dev. (PW)')
title({'POP ID standard deviation of heat transport from spatial components,'; ['globally integrated, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig17,['Heat_transport_spatial_components_POP_global_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_stddev_ID.pdf'])
close(fig17)


fig18 = figure(18);
close(figure(1))
h = plot(lat_transects_vec_global,(1e-15)*(heat_transport_mesoscale_tmean_tmean_global + heat_transport_mesoscale_lowfreq_tmean_global + heat_transport_mesoscale_highfreq_tmean_global),lat_transects_vec_global,(1e-15)*heat_transport_vT_timevar_global_tmean,lat_transects_vec_global,(1e-15)*heat_transport_vT_highfreq_global_tmean);
leg = legend('Mesoscale','Time-varying','High freq.','location','NorthEast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec_global) max(lat_transects_vec_global)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[-1 1],'ytick',(-1):0.25:1)
set(gca,'xticklabels',{''; ''; '-60'; ''; '-40'; ''; '-20'; ''; '0'; ''; '20'; ''; '40'; ''; '60'; ''})
set(gca,'yticklabels',{''; ''; '-0.5'; ''; '-0'; ''; '0.5'; ''; '1'})
daspect([(((20/7)*aspect_ratio_baseline) + ((2/7)*(max(lat_transects_vec_global) - min(lat_transects_vec_global)))) 1 1])
set(h(1),'Color',[0.8 0 0],'LineWidth',2)
set(h(2),'Color',[0.4 0.3 0.7],'LineWidth',2)
set(h(3),'Color',[0.8 0 0.8],'LineWidth',2,'LineStyle','-.')
hold on
line([min(lat_transects_vec_global) max(lat_transects_vec_global)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('Heat transport (PW)')
title({'POP total time-averaged heat transport from various "eddy" formulations,'; ['globally integrated, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig18,['Heat_transport_eddy_formulations_POP_global_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_tavg.pdf'])
close(fig18)


fig18 = figure(18);
close(figure(1))
h = plot(lat_transects_vec_global,(1e-15)*heat_transport_mesoscale_global_stddev_ID,lat_transects_vec_global,(1e-15)*heat_transport_vT_timevar_global_stddev_ID,lat_transects_vec_global,(1e-15)*heat_transport_vT_highfreq_global_stddev_ID);
leg = legend('Mesoscale','Time-varying','High freq.','location','NorthEast');
set(leg,'FontSize',12)
set(gca,'FontSize',12,'xlim',[min(lat_transects_vec_global) max(lat_transects_vec_global)],'xtick',(-80):xtick_spacing:80,'xgrid','on','ylim',[0 0.2],'ytick',0:0.05:0.2)
set(gca,'xticklabels',{''; ''; '-60'; ''; '-40'; ''; '-20'; ''; '0'; ''; '20'; ''; '40'; ''; '60'; ''})
set(gca,'yticklabels',{'0'; '0.05'; '0.1'; '0.15'; '0.2'})
daspect([(((180/7)*aspect_ratio_baseline) + ((18/7)*(max(lat_transects_vec_global) - min(lat_transects_vec_global)))) 1 1])
set(h(1),'Color',[0.8 0 0],'LineWidth',2)
set(h(2),'Color',[0.4 0.3 0.7],'LineWidth',2)
set(h(3),'Color',[0.8 0 0.8],'LineWidth',2,'LineStyle','-.')
hold on
line([min(lat_transects_vec_global) max(lat_transects_vec_global)],[0 0],[1 1],'Color','k','LineStyle','-')
hold off
xlabel('Latitude')
ylabel('ID heat transport std. dev. (PW)')
title({'POP ID std. dev. of heat transport from various "eddy" formulations,'; ['globally integrated, error fcn filter scheme with ',num2str(round(time_scale_separation)),' day time separation']; ''},'FontSize',8)

saveas(fig18,['Heat_transport_eddy_formulations_POP_global_',spatial_scheme_file_id,'_spatsep_',datestr_start,'_',datestr_end_minus1,'_stddev_ID.pdf'])
close(fig18)

% % % % 


% plot time series of spatial components at a given latitude

basin_id = 'Atl';
lat_to_plot = 63;

if strcmp(basin_id,'global') == 1
    curr_lat_ind = find(abs(lat_transects_vec_global - lat_to_plot) < 0.5);
    curr_zonmean_tseries_ID = heat_transport_zonmean_tseries_array_global_ID(curr_lat_ind,:);
    curr_largescale_tseries_ID = heat_transport_largescale_tseries_array_global_ID(curr_lat_ind,:);        
    curr_mesoscale_tseries_ID = heat_transport_mesoscale_tseries_array_global_ID(curr_lat_ind,:);        
    curr_residual_tseries_ID = heat_transport_residual_tseries_array_global_ID(curr_lat_ind,:);        
    curr_all_tseries_ID = heat_transport_all_tseries_array_global_ID(curr_lat_ind,:);
    curr_zonmean_tseries = heat_transport_zonmean_tseries_array_global(curr_lat_ind,:);
    curr_largescale_tseries = heat_transport_largescale_tseries_array_global(curr_lat_ind,:);
    curr_mesoscale_tseries = heat_transport_mesoscale_tseries_array_global(curr_lat_ind,:);
    curr_residual_tseries = heat_transport_residual_tseries_array_global(curr_lat_ind,:);
    curr_all_tseries = heat_transport_all_tseries_array_global(curr_lat_ind,:);
else
    curr_lat_ind = find(abs(lat_transects_vec - lat_to_plot) < 0.5);
    n_bins = round(365/delta_t);
    curr_zonmean_tseries = eval(['heat_transport_zonmean_tseries_array_',basin_id,'(curr_lat_ind,:)']);
    curr_largescale_tseries = eval(['heat_transport_largescale_tseries_array_',basin_id,'(curr_lat_ind,:)']);
    curr_mesoscale_tseries = eval(['heat_transport_mesoscale_tseries_array_',basin_id,'(curr_lat_ind,:)']);
    curr_residual_tseries = eval(['heat_transport_residual_tseries_array_',basin_id,'(curr_lat_ind,:)']);
    curr_all_tseries = eval(['heat_transport_all_tseries_array_',basin_id,'(curr_lat_ind,:)']);
    curr_zonmean_tseries_noseason = NaN(size(curr_zonmean_tseries));
    curr_largescale_tseries_noseason = NaN(size(curr_largescale_tseries));
    curr_mesoscale_tseries_noseason = NaN(size(curr_mesoscale_tseries));
    curr_residual_tseries_noseason = NaN(size(curr_residual_tseries));
    curr_all_tseries_noseason = NaN(size(curr_all_tseries));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        curr_nan_mask_1 = ones(size(curr_largescale_tseries(curr_in_bin_ind)));
        curr_nan_mask_2 = ones(size(curr_mesoscale_tseries(curr_in_bin_ind)));
        curr_nan_mask_1((isnan(curr_largescale_tseries(curr_in_bin_ind)) == 1) | (abs(curr_largescale_tseries(curr_in_bin_ind)) < 1e-5)) = 0;
        curr_nan_mask_2((isnan(curr_mesoscale_tseries(curr_in_bin_ind)) == 1) | (abs(curr_mesoscale_tseries(curr_in_bin_ind)) < 1e-5)) = 0;
        curr_nan_mask = curr_nan_mask_1 & curr_nan_mask_2;

        curr_array = curr_zonmean_tseries(curr_in_bin_ind);
        curr_array(abs(curr_nan_mask) < 1e-5) = 0;
        curr_zonmean_tseries_noseason(curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
        curr_array = curr_largescale_tseries(curr_in_bin_ind);
        curr_array(abs(curr_nan_mask) < 1e-5) = 0;
        curr_largescale_tseries_noseason(curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
        curr_array = curr_mesoscale_tseries(curr_in_bin_ind);
        curr_array(abs(curr_nan_mask) < 1e-5) = 0;
        curr_mesoscale_tseries_noseason(curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
        curr_array = curr_residual_tseries(curr_in_bin_ind);
        curr_array(abs(curr_nan_mask) < 1e-5) = 0;
        curr_residual_tseries_noseason(curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
        curr_array = curr_all_tseries(curr_in_bin_ind);
        curr_array(abs(curr_nan_mask) < 1e-5) = 0;
        curr_all_tseries_noseason(curr_in_bin_ind) = curr_array - repmat(sum(curr_nan_mask.*curr_array,2)./(sum(curr_nan_mask,2)),[1 size(curr_array,2)]);
    end


    % filter time series for interannual & decadal frequencies

    half_power_adj = exp(erfinv((2^(1/2)) - 1)/5);   % adjustment factor to set bounds at half-power (rather than half-amplitude)
    [curr_zonmean_tseries_ID,~,~] = bandpass_err_fcn(curr_zonmean_tseries_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
    [curr_largescale_tseries_ID,~,~] = bandpass_err_fcn(curr_largescale_tseries_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
    [curr_mesoscale_tseries_ID,~,~] = bandpass_err_fcn(curr_mesoscale_tseries_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
    [curr_residual_tseries_ID,~,~] = bandpass_err_fcn(curr_residual_tseries_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
    [curr_all_tseries_ID,~,~] = bandpass_err_fcn(curr_all_tseries_noseason,2,delta_t,(1/(3*sum(diff(time_datenum_in_range))))/half_power_adj,(1/426)*half_power_adj,5,1,1,1,0);
    curr_zonmean_tseries_ID = curr_zonmean_tseries_ID + repmat(mean(curr_zonmean_tseries,2),[1 length(time_datenum_in_range)]);
    curr_largescale_tseries_ID = curr_largescale_tseries_ID + repmat(mean(curr_largescale_tseries,2),[1 length(time_datenum_in_range)]);
    curr_mesoscale_tseries_ID = curr_mesoscale_tseries_ID + repmat(mean(curr_mesoscale_tseries,2),[1 length(time_datenum_in_range)]);
    curr_residual_tseries_ID = curr_residual_tseries_ID + repmat(mean(curr_residual_tseries,2),[1 length(time_datenum_in_range)]);
    curr_all_tseries_ID = curr_all_tseries_ID + repmat(mean(curr_all_tseries,2),[1 length(time_datenum_in_range)]);

end

curr_zonmean_tseries_ID_minus_tmean = curr_zonmean_tseries_ID - repmat(mean(curr_zonmean_tseries,2),[1 length(time_datenum_in_range)]);
curr_largescale_tseries_ID_minus_tmean = curr_largescale_tseries_ID - repmat(mean(curr_largescale_tseries,2),[1 length(time_datenum_in_range)]);
curr_mesoscale_tseries_ID_minus_tmean = curr_mesoscale_tseries_ID - repmat(mean(curr_mesoscale_tseries,2),[1 length(time_datenum_in_range)]);
curr_residual_tseries_ID_minus_tmean = curr_residual_tseries_ID - repmat(mean(curr_residual_tseries,2),[1 length(time_datenum_in_range)]);
curr_all_tseries_ID_minus_tmean = curr_all_tseries_ID - repmat(mean(curr_all_tseries,2),[1 length(time_datenum_in_range)]);


field_plot_scale_factor = 1e15;

fig4 = figure(4);
% h = plot(time_datenum_in_range,(1/field_plot_scale_factor)*curr_zonmean_tseries_ID,time_datenum_in_range,(1/field_plot_scale_factor)*curr_largescale_tseries_ID,time_datenum_in_range,(1/field_plot_scale_factor)*curr_mesoscale_tseries_ID,time_datenum_in_range,(1/field_plot_scale_factor)*curr_residual_tseries_ID,time_datenum_in_range,(1/field_plot_scale_factor)*curr_all_tseries_ID);
h = plot(time_datenum_in_range,(1/field_plot_scale_factor)*curr_zonmean_tseries_ID_minus_tmean,time_datenum_in_range,(1/field_plot_scale_factor)*curr_largescale_tseries_ID_minus_tmean,time_datenum_in_range,(1/field_plot_scale_factor)*curr_mesoscale_tseries_ID_minus_tmean,time_datenum_in_range,(1/field_plot_scale_factor)*curr_residual_tseries_ID_minus_tmean,time_datenum_in_range,(1/field_plot_scale_factor)*curr_all_tseries_ID_minus_tmean);
set(h(1),'Color',[0.55 0 0.55],'LineWidth',2)
set(h(2),'Color',[0 0 0.8],'LineWidth',2)
set(h(3),'Color',[0.8 0 0],'LineWidth',2)
set(h(4),'Color',[0.8 0.6 0.4],'LineWidth',2)
set(h(5),'Color',[0 0 0],'LineWidth',2)
years_to_plot_ticks = ((ceil(time_range_start(1)/5)*5):5:(ceil(time_range_end(1)/5)*5))';
xtick_datenums_plot = 365*(years_to_plot_ticks - year_offset);
xtick_labels_plot = cell(length(years_to_plot_ticks),1);
for xtick_ind = 1:length(years_to_plot_ticks)
%     xtick_labels_plot{xtick_ind} = [num2str(years_to_plot_ticks(xtick_ind)),'-01-01'];
    xtick_labels_plot{xtick_ind} = num2str(years_to_plot_ticks(xtick_ind));
end
% n_ticks = (3*(time_range_end(1) - time_range_start(1) + 1)) + 1;
% xtick_datenums_plot = (365*time_range_start(1)) + (datenum([(1990*ones([n_ticks 1])) (1:4:(((n_ticks - 1)*4) + 1))' ones([n_ticks 1]) zeros([n_ticks 3])]) - datenum([1990 1 1 0 0 0])) - ceil((datenum([(1990*ones([n_ticks 1])) (1:4:(((n_ticks - 1)*4) + 1))' ones([n_ticks 1]) zeros([n_ticks 3])]) - datenum([1992 2 29 0 0 0]))/(4*365.25)); 
% xtick_labels_plot = cell(length(xtick_datenums_plot),1);
% for xtick_ind = 1:length(xtick_datenums_plot)
%     xtick_labels_plot{xtick_ind} = ['00',num2str(floor(xtick_datenums_plot(xtick_ind)/365)),datestr([1990 01 mod(xtick_datenums_plot(xtick_ind),365) 0 0 0],'-mm-dd')];
% end
set(gca,'FontSize',12,'xlim',[datenum_start datenum_end],'xtick',xtick_datenums_plot,'xticklabel',xtick_labels_plot,'xgrid','on','ylim',[-0.5 0.5],'ytick',(-1):0.2:1,'DataAspectRatio',[(12*365) 1 1])
% set(get(gca,'ylabel'),'String','Heat transport (PW)','FontSize',14)
set(get(gca,'ylabel'),'String','MHT anomaly (PW)','FontSize',14)
hold on
line([datenum_start datenum_end],[0 0],[0 0],'Color','k','LineStyle','-','LineWidth',1)
hold off
title({['POP time series of MHT spatial components across ',num2str(lat_to_plot),' ^{\circ} latitude, ',basin_id]; ''},'FontSize',10)
leg = legend('Overturning','Large scale','Mesoscale','Residual','Total','location','southoutside');
keyboard
% print(fig4,['Transect_tseries_POP_',basin_id,'_',num2str(lat_to_plot),'_lat_',datestr(datenum(time_range_start),'yyyymmdd'),'_',datestr(datenum(time_range_end) - 1,'yyyymmdd'),'.pdf'],'-dpdf','-r300')
print(fig4,['Transect_tseries_POP_nomean_',basin_id,'_',num2str(lat_to_plot),'_lat_',datestr(datenum(time_range_start),'yyyymmdd'),'_',datestr(datenum(time_range_end) - 1,'yyyymmdd'),'.pdf'],'-dpdf','-r300')
close(fig4)


fig4 = figure(4);
% h = plot(time_datenum_in_range,(1/field_plot_scale_factor)*curr_largescale_tseries_ID,time_datenum_in_range,(1/field_plot_scale_factor)*curr_mesoscale_tseries_ID,time_datenum_in_range,(1/field_plot_scale_factor)*curr_residual_tseries_ID,time_datenum_in_range,(1/field_plot_scale_factor)*(curr_all_tseries_ID - curr_zonmean_tseries_ID));
h = plot(time_datenum_in_range,(1/field_plot_scale_factor)*curr_largescale_tseries_ID_minus_tmean,time_datenum_in_range,(1/field_plot_scale_factor)*curr_mesoscale_tseries_ID_minus_tmean,time_datenum_in_range,(1/field_plot_scale_factor)*curr_residual_tseries_ID_minus_tmean,time_datenum_in_range,(1/field_plot_scale_factor)*(curr_all_tseries_ID_minus_tmean - curr_zonmean_tseries_ID_minus_tmean));
set(h(1),'Color',[0 0 0.8],'LineWidth',2)
set(h(2),'Color',[0.8 0 0],'LineWidth',2)
set(h(3),'Color',[0.8 0.6 0.4],'LineWidth',2)
set(h(4),'Color',[0 0 0],'LineWidth',2)
years_to_plot_ticks = ((ceil(time_range_start(1)/5)*5):5:(ceil(time_range_end(1)/5)*5))';
xtick_datenums_plot = 365*(years_to_plot_ticks - year_offset);
xtick_labels_plot = cell(length(years_to_plot_ticks),1);
for xtick_ind = 1:length(years_to_plot_ticks)
%     xtick_labels_plot{xtick_ind} = [num2str(years_to_plot_ticks(xtick_ind)),'-01-01'];
    xtick_labels_plot{xtick_ind} = num2str(years_to_plot_ticks(xtick_ind));
end
% n_ticks = (3*(time_range_end(1) - time_range_start(1) + 1)) + 1;
% xtick_datenums_plot = (365*time_range_start(1)) + (datenum([(1990*ones([n_ticks 1])) (1:4:(((n_ticks - 1)*4) + 1))' ones([n_ticks 1]) zeros([n_ticks 3])]) - datenum([1990 1 1 0 0 0])) - ceil((datenum([(1990*ones([n_ticks 1])) (1:4:(((n_ticks - 1)*4) + 1))' ones([n_ticks 1]) zeros([n_ticks 3])]) - datenum([1992 2 29 0 0 0]))/(4*365.25)); 
% xtick_labels_plot = cell(length(xtick_datenums_plot),1);
% for xtick_ind = 1:length(xtick_datenums_plot)
%     xtick_labels_plot{xtick_ind} = ['00',num2str(floor(xtick_datenums_plot(xtick_ind)/365)),datestr([1990 01 mod(xtick_datenums_plot(xtick_ind),365) 0 0 0],'-mm-dd')];
% end
set(gca,'FontSize',12,'xlim',[datenum_start datenum_end],'xtick',xtick_datenums_plot,'xticklabel',xtick_labels_plot,'xgrid','on','ylim',[-0.3 0.3],'ytick',(-1):0.1:1,'DataAspectRatio',[(20*365) 1 1])
% set(get(gca,'ylabel'),'String','Heat transport (PW)','FontSize',14)
set(get(gca,'ylabel'),'String','MHT anomaly (PW)','FontSize',14)
hold on
line([datenum_start datenum_end],[0 0],[0 0],'Color','k','LineStyle','-','LineWidth',1)
hold off
title({['POP time series of MHT spatial components across ',num2str(lat_to_plot),' ^{\circ} latitude, ',basin_id,', no overturn']; ''},'FontSize',10)
leg = legend('Large scale','Mesoscale','Residual','Total - overturning','location','southoutside');
keyboard
% print(fig4,['Transect_tseries_POP_minusoverturn_',basin_id,'_',num2str(lat_to_plot),'_lat_',datestr(datenum(time_range_start),'yyyymmdd'),'_',datestr(datenum(time_range_end) - 1,'yyyymmdd'),'.pdf'],'-dpdf','-r300')
print(fig4,['Transect_tseries_POP_minusoverturn_nomean_',basin_id,'_',num2str(lat_to_plot),'_lat_',datestr(datenum(time_range_start),'yyyymmdd'),'_',datestr(datenum(time_range_end) - 1,'yyyymmdd'),'.pdf'],'-dpdf','-r300')
close(fig4)