POP_heat_transport_multitransects_EKE_mask.m

% compute meridional temperature fluxes and MHT, masking out some times so that the time-mean model EKE approximates satellite estimates

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


% transect .mat file to load fluxes from
transect_file = '/indopac/adelman/Global_mesoscale/POP/Multitransect_arrays_POP_IndoPac_1deg.mat';
plot_name = 'IndoPac_1deg';

% time period to include
time_range_start = [1978 01 01];
time_range_end = [2010 01 01];

% minimum number of days required in unmasked time (in adj_alt fields)
min_days_alt = 3*365;


% nominal grid spacing (in degrees)
grid_res = 0.1;

% model year offset
year_offset = 1976;


gridded_SSH_fileform_1 = '/indopac/adelman/POP/pop.h.nday5.JC.SSH_global.';
gridded_SSH_fileform_2 = '.nc';
delta_t = 5;


% load meridional fluxes

load(transect_file,'lat_transects_vec','lon_bounds_array','lon_in_range_vel_all','lon_in_range_all','time_datenum_in_range','spatial_scale_separation_vec','cumint_heat_flux_mesoscale*across_transect_array','integrated_heat_flux_mesoscale*tmean_array')

cumint_heat_flux_mesoscale_all_across_transect_array = cumint_heat_flux_mesoscale_tmean_across_transect_array + cumint_heat_flux_mesoscale_lowfreq_across_transect_array + cumint_heat_flux_mesoscale_highfreq_across_transect_array;

local_heat_flux_mesoscale_all_array = diff(cumint_heat_flux_mesoscale_all_across_transect_array,1,2)./(111100*repmat(cosd(lat_transects_vec),[1 length(lon_in_range_all) length(time_datenum_in_range)]).*repmat(diff(lon_in_range_vel_all'),[length(lat_transects_vec) 1 length(time_datenum_in_range)]));
% land masks
land_mask_local_plot = ones([length(lat_transects_vec) length(lon_in_range_all)]);
for transect_ind = 1:length(lat_transects_vec)
    land_mask_local_plot(transect_ind,(isnan(cumint_heat_flux_mesoscale_all_across_transect_array(transect_ind,1:length(lon_in_range_all),ceil(length(time_datenum_in_range)/2))) == 1) & (isnan(cumint_heat_flux_mesoscale_all_across_transect_array(transect_ind,2:(length(lon_in_range_all) + 1),ceil(length(time_datenum_in_range)/2))) == 1)) = 0;
    diff2_mesoscale_hflux = diff(cumint_heat_flux_mesoscale_all_across_transect_array(transect_ind,:,ceil(length(time_datenum_in_range)/2)),2,2);
    diff2_0_ind = find(abs(diff2_mesoscale_hflux) < 1);
    gap_ind = find(diff(diff2_0_ind) > 1.5);
    if isempty(diff2_0_ind) == 1
        continue
    end
    seg_start_ind = diff2_0_ind([1 (gap_ind + 1)]);
    seg_end_ind = 1 + diff2_0_ind([gap_ind length(diff2_0_ind)]);
    for curr_seg_ind = 1:length(seg_start_ind)
        land_mask_local_plot(transect_ind,seg_start_ind(curr_seg_ind):seg_end_ind(curr_seg_ind)) = 0;
    end
end


n_files = time_range_end(1) - time_range_start(1);

time_bound = NaN([0 2]);
file_num_vec = [];
in_file_ind_vec = [];
for file_num_ind = 1:n_files
    curr_model_yr = time_range_start(1) + (file_num_ind - 1) - year_offset;
    if curr_model_yr < 10
        curr_yearstr = ['000',num2str(curr_model_yr)];
    else
        curr_yearstr = ['00',num2str(curr_model_yr)];
    end
    curr_filename = [gridded_SSH_fileform_1,curr_yearstr,gridded_SSH_fileform_2];
    
    time_bound_curr_source_file = (ncread(curr_filename,'time_bound'))';
    
    time_bound = [time_bound; time_bound_curr_source_file];
    
    file_num_vec = [file_num_vec; (file_num_ind*ones([size(time_bound_curr_source_file,1) 1]))];
    in_file_ind_vec = [in_file_ind_vec; ((1:1:size(time_bound_curr_source_file,1))')];
end

time_datenum = mean(time_bound,2);

datenum_start = (365*(time_range_start(1) - year_offset)) + (datenum([1990 time_range_start(2:3)]) - datenum([1990 01 01]));
datenum_end = (365*(time_range_end(1) - year_offset)) + (datenum([1990 time_range_end(2:3)]) - datenum([1990 01 01]));


% % % % 

% extract POP output

curr_yearstr = '0002';
curr_filename = [gridded_SSH_fileform_1,curr_yearstr,gridded_SSH_fileform_2];
longitude = ncread(curr_filename,'TLONG');
land_ind = find(longitude < -1e-5);
longitude(land_ind) = NaN;
latitude = ncread(curr_filename,'TLAT');
latitude(land_ind) = NaN;

time_datenum_range_start = datenum(time_range_start);
time_datenum_range_end = datenum(time_range_end);


in_time_range_ind = find((time_datenum - datenum_start >= -1e-5) & (time_datenum - datenum_end <= 1e-5));

% remove possible duplication of times across files
[unique_times,unique_time_in_range_ind,~] = unique(time_datenum(in_time_range_ind));
[~,sorted_unique_ind] = sort(unique_times,'ascend');



nc_file_name_satobs = '/indopac/adelman/Global_mesoscale/SSH/SSH_lp_elong_eddy_SSALTO_DUACS_5day_global.nc';
grid_res_satobs = 1/4;
lon_satobs = ncread(nc_file_name_satobs,'longitude');
lat_satobs = ncread(nc_file_name_satobs,'latitude');
time_datenum_satobs = ncread(nc_file_name_satobs,'time');


EKE_along_transect_smoothed_array = NaN([length(lat_transects_vec) length(lon_in_range_all) length(time_datenum_in_range)]);
EKE_mean_along_transect_smoothed_array = NaN([length(lat_transects_vec) length(lon_in_range_all)]);
EKE_mean_satobs_model_grid_smoothed_array = NaN([length(lat_transects_vec) length(lon_in_range_all)]);
EKE_low_threshold_adj = NaN([length(lat_transects_vec) length(lon_in_range_all)]);
EKE_high_threshold_adj = NaN([length(lat_transects_vec) length(lon_in_range_all)]);
EKE_low_threshold_adj_alt = NaN([length(lat_transects_vec) length(lon_in_range_all)]);
EKE_high_threshold_adj_alt = NaN([length(lat_transects_vec) length(lon_in_range_all)]);
n_days_unmasked_along_transect = NaN([length(lat_transects_vec) length(lon_in_range_all)]);
for transect_ind = 1:length(lat_transects_vec)
    lat_transect = lat_transects_vec(transect_ind);
    
    if abs(lat_transect) < 5.1
        continue
    end
    
    lon_min = lon_in_range_all(find(abs(land_mask_local_plot(transect_ind,:)) > 1e-5,1,'first'));
    lon_max = lon_in_range_all(find(abs(land_mask_local_plot(transect_ind,:)) > 1e-5,1,'last'));
    longwest = lon_min - 1;
    longeast = lon_max + 1;
    latsouth = lat_transect - 1;
    latnorth = lat_transect + 1;
    
    if lat_transect < 28.5
        
        near_lat_ind = find((latitude - latsouth >= 0) & (latitude - latnorth <= 0));
        [~,nearest_west_ind] = min(abs((mod(longitude(near_lat_ind) - longwest + 180,360) - 180) + (1i*(latitude(near_lat_ind) - lat_transect))));
        nearest_west_ind = near_lat_ind(nearest_west_ind);
        [~,nearest_east_ind] = min(abs((mod(longitude(near_lat_ind) - longeast + 180,360) - 180) + (1i*(latitude(near_lat_ind) - lat_transect))));
        nearest_east_ind = near_lat_ind(nearest_east_ind);
        
        nearest_west_i_ind = mod(nearest_west_ind - 1,size(longitude,1)) + 1;
        nearest_east_i_ind = mod(nearest_east_ind - 1,size(longitude,1)) + 1;
        near_lat_j_ind = ceil(near_lat_ind/size(latitude,1));
        min_j_ind = min(near_lat_j_ind);
        max_j_ind = max(near_lat_j_ind);
        
        if nearest_east_i_ind - nearest_west_i_ind < (3/grid_res)
            in_lon_range_ind = [(nearest_west_i_ind:1:size(longitude,1))'; (1:1:nearest_east_i_ind)'];
        else
            in_lon_range_ind = (nearest_west_i_ind:1:nearest_east_i_ind)';
        end
        in_lat_range_ind = (min_j_ind:1:max_j_ind)';
        
        lon_in_range = longitude(in_lon_range_ind,in_lat_range_ind);
        lat_in_range = latitude(in_lon_range_ind,in_lat_range_ind);
        time_datenum_in_range_all = time_datenum(in_time_range_ind);
        time_bound_datenum_in_range_all = time_bound(in_time_range_ind,:);
        file_num_in_range_vec_all = file_num_vec(in_time_range_ind);
        in_file_ind_in_range_vec_all = in_file_ind_vec(in_time_range_ind);
        
        
        size_array = [length(in_lon_range_ind) length(in_lat_range_ind) 1 length(in_time_range_ind)];
        
        
        lon_in_range(isnan(lon_in_range) == 1) = 1e4;
        diff_lon_in_range = diff([lon_in_range; lon_in_range(1,:)],1,1); 
        diff_lon_in_range((abs(diff_lon_in_range) > 180) & (abs(diff_lon_in_range) < 5e3)) = mod(diff_lon_in_range((abs(diff_lon_in_range) > 180) & (abs(diff_lon_in_range) < 5e3)) + 180,360) - 180;
        lon_in_range_aug_x = repmat(lon_in_range(1,:),[(size_array(1) + 1) 1]) + [zeros([1 size_array(2)]); cumsum(diff_lon_in_range,1,'forward')];
        lon_in_range_aug_x(lon_in_range_aug_x > 5e3) = NaN;
        lon_in_range_aug_x(lon_in_range_aug_x < min(lon_in_range_aug_x(1,:))) = lon_in_range_aug_x(lon_in_range_aug_x < min(lon_in_range_aug_x(1,:))) + (360*ceil((min(lon_in_range_aug_x(1,:)) - lon_in_range_aug_x(lon_in_range_aug_x < min(lon_in_range_aug_x(1,:))))/360));
        lon_in_range_aug_x(size_array(1) + 1,lon_in_range_aug_x(size_array(1) + 1,:) < min(lon_in_range_aug_x(1,:)) + 180) = lon_in_range_aug_x(size_array(1) + 1,lon_in_range_aug_x(size_array(1) + 1,:) < min(lon_in_range_aug_x(1,:)) + 180) + 360;
        lon_in_range_aug_x(lon_in_range_aug_x > max(lon_in_range_aug_x(size_array(1) + 1,:))) = lon_in_range_aug_x(lon_in_range_aug_x > max(lon_in_range_aug_x(size_array(1) + 1,:))) + (360*floor((max(lon_in_range_aug_x(size_array(1) + 1,:)) - lon_in_range_aug_x(lon_in_range_aug_x > max(lon_in_range_aug_x(size_array(1) + 1,:))))/360));
        lon_in_range = lon_in_range_aug_x(1:size_array(1),:);
        lon_in_range(abs(lon_in_range - 1e4) < 1e-5) = NaN;
        lon_in_range_midpts = lon_in_range_aug_x(1:size_array(1),:) + (diff(lon_in_range_aug_x,1,1)/2);
        lon_in_range_midpts = [(lon_in_range_midpts(:,2:size_array(2)) - (diff(lon_in_range_midpts,1,2)/2)) NaN([size_array(1) 1])];
        lat_in_range_aug_x = [lat_in_range; lat_in_range(1,:)];
        lat_in_range_midpts = lat_in_range_aug_x(1:size_array(1),:) + (diff(lat_in_range_aug_x,1,1)/2);
        lat_in_range_midpts = [(lat_in_range_midpts(:,2:size_array(2)) - (diff(lat_in_range_midpts,1,2)/2)) NaN([size_array(1) 1])];
        dx = 111100*(abs((cosd(lat_in_range_midpts).*diff([(lon_in_range_aug_x(:,2:size_array(2)) - (diff(lon_in_range_aug_x,1,2)/2)) NaN([(size_array(1) + 1) 1])],1,1)) + (1i*diff([(lat_in_range_aug_x(:,2:size_array(2)) - (diff(lat_in_range_aug_x,1,2)/2)) NaN([(size_array(1) + 1) 1])],1,1))));
        dy = 111100*(abs((cosd(lat_in_range_midpts).*diff([(lon_in_range_aug_x(1:size_array(1),:) + (diff(lon_in_range_aug_x,1,1)/2)) NaN([size_array(1) 1])],1,2)) + (1i*[diff(lat_in_range_aug_x(1:size_array(1),:) + (diff(lat_in_range_aug_x,1,1)/2),1,2) NaN([size_array(1) 1])])));
        
%         in_time_range_ind = in_time_range_ind_all;
%         time_datenum_in_range = time_datenum_in_range_all;
        file_num_in_range_vec = file_num_in_range_vec_all;
        in_file_ind_in_range_vec = in_file_ind_in_range_vec_all;
        
        
        size_array = [length(in_lon_range_ind) length(in_lat_range_ind) 1 length(in_time_range_ind)];
%         delta_t = mean(diff(time_datenum_in_range));
        
        
        SSH_in_range = NaN(size_array([1 2 4]));
        for file_num_ind = 1:(max(file_num_in_range_vec) - min(file_num_in_range_vec) + 1)
            curr_model_yr = time_range_start(1) + (file_num_ind - 1) - year_offset;
            if curr_model_yr < 10
                curr_yearstr = ['000',num2str(curr_model_yr)];
            else
                curr_yearstr = ['00',num2str(curr_model_yr)];
            end
            curr_filename = [gridded_SSH_fileform_1,curr_yearstr,gridded_SSH_fileform_2];
    
            curr_time_subset_ind = find(file_num_in_range_vec == file_num_ind);
            if isempty(curr_time_subset_ind) == 1
                continue
            end
            curr_in_time_range_ind = in_file_ind_in_range_vec(curr_time_subset_ind);
            
            if min(diff(in_lon_range_ind)) < 0
                gap_ind = find(diff(in_lon_range_ind) < 0);
                
                start_vec = [in_lon_range_ind(1) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
                count_vec = [(diff(in_lon_range_ind([1 gap_ind])) + 1) (max(in_lat_range_ind) - min(in_lat_range_ind) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                curr_SSH_in_range = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
                start_vec_2 = [in_lon_range_ind(gap_ind + 1) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
                count_vec_2 = [(diff(in_lon_range_ind([(gap_ind + 1) length(in_lon_range_ind)])) + 1) (max(in_lat_range_ind) - min(in_lat_range_ind) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                curr_SSH_in_range_2 = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec_2([1 2 4]),count_vec_2([1 2 4]))));
                curr_SSH_in_range = [curr_SSH_in_range; curr_SSH_in_range_2];
            else
                lon_range_ind_span = max(in_lon_range_ind) - min(in_lon_range_ind) + 1;
                lat_range_ind_span = max(in_lat_range_ind) - min(in_lat_range_ind) + 1;
                time_range_ind_span = max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1;
        
                start_vec = [min(in_lon_range_ind) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
                count_vec = [lon_range_ind_span lat_range_ind_span 1 time_range_ind_span];
        
                curr_SSH_in_range = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
            end
            curr_SSH_in_range(abs(curr_SSH_in_range) > 1e15) = 0;
            SSH_in_range(:,:,curr_time_subset_ind) = curr_SSH_in_range(:,:,curr_in_time_range_ind - min(curr_in_time_range_ind) + 1);
    
        end
        
        
        size_array = size(SSH_in_range);
        
        nan_mask = ones(size(SSH_in_range));
        nan_mask(isnan(SSH_in_range) == 1) = 0;
        SSH_in_range_nans_zeroed = SSH_in_range;
        SSH_in_range_nans_zeroed(isnan(SSH_in_range) == 1) = 0;
        sum_nan_mask = sum(nan_mask,3);
        nan_mask_not_enough = ones(size(sum_nan_mask));
        nan_mask_not_enough(sum_nan_mask < (0.9*max(max(sum_nan_mask)))) = 0;
        nan_mask_not_enough_3D = repmat(nan_mask_not_enough,[1 1 size(nan_mask,3)]);
        nan_mask_SSH = nan_mask & nan_mask_not_enough_3D;
        
        clear nan_mask nan_mask_not_enough*
        
        SSH_mean_in_range = sum(nan_mask_SSH.*SSH_in_range_nans_zeroed,3)./sum(nan_mask_SSH,3);
        
        SSH_nomean_in_range = SSH_in_range - repmat(SSH_mean_in_range,[1 1 size(SSH_in_range,3)]);
        
        
        diff_x_SSH_nomean_in_r = diff(SSH_nomean_in_range([(1:1:size_array(1))'; 1],:,:),1,1);
        diff_x_SSH_nomean_in_range = [(diff_x_SSH_nomean_in_r(:,2:size(diff_x_SSH_nomean_in_r,2),:) - (diff(diff_x_SSH_nomean_in_r,1,2)/2)) NaN([size_array(1) 1 size_array(3)])];
        diff_y_SSH_nomean_in_r = diff(SSH_nomean_in_range,1,2);
        diff_y_SSH_nomean_in_range = [(diff_y_SSH_nomean_in_r([(2:1:size(diff_y_SSH_nomean_in_r,1))'; 1],:,:) - (diff(diff_y_SSH_nomean_in_r([(1:1:size(diff_y_SSH_nomean_in_r,1))'; 1],:,:),1,1)/2)) NaN([size_array(1) 1 size_array(3)])];
    
    
        f = 2*((2*pi)/86164)*sind(lat_in_range_midpts);
        
        uvel_geostr_nomean_in_range = -9.81*repmat((f.^(-1)),[1 1 size(SSH_nomean_in_range,3)]).*(diff_y_SSH_nomean_in_range./repmat(dy,[1 1 size(SSH_nomean_in_range,3)]));
        vvel_geostr_nomean_in_range = 9.81*repmat((f.^(-1)),[1 1 size(SSH_nomean_in_range,3)]).*(diff_x_SSH_nomean_in_range./repmat(dx,[1 1 size(SSH_nomean_in_range,3)]));
    
        clear diff_*
        
        EKE_in_range = (0.5*((uvel_geostr_nomean_in_range.^2) + (vvel_geostr_nomean_in_range.^2)));
    
        clear uvel_geostr* vvel_geostr*
    
        % mask out low-latitude points
        lat_mask = ones(size(lat_in_range_midpts));
        lat_mask(abs(lat_in_range_midpts) < 5) = 0;
        lat_mask_3d = repmat(lat_mask,[1 1 size(EKE_in_range,3)]);
        lat_mask_ind = find(abs(lat_mask_3d) < 1e-5);
    
        EKE_in_range(lat_mask_ind) = NaN;
        
        % interpolate to transect points
        
%         EKE_in_range((isnan(EKE_in_range) == 1) | (isinf(EKE_in_range) == 1)) = 0;
        
        lon_pts_in_range_ind = find((lon_in_range_all >= lon_min) & (lon_in_range_all <= lon_max));
        EKE_along_transect = NaN([1 length(lon_in_range_all) length(time_datenum_in_range)]);
%         for t = 1:length(time_datenum_in_range)
%             EKE_along_transect(1,lon_pts_in_range_ind,t) = (interp2_fft(mod(lon_in_range_midpts(:,ceil(size(lon_in_range_midpts,2)/2)) - min(lon_in_range_all(lon_pts_in_range_ind)),360) + min(lon_in_range_all(lon_pts_in_range_ind)),lat_in_range_midpts(ceil(size(lon_in_range_midpts,1)/2),1:(size(lat_in_range_midpts,2) - 1)),EKE_in_range(:,1:(size(lat_in_range_midpts,2) - 1),t),lon_in_range_all(lon_pts_in_range_ind),lat_transect*ones([length(lon_pts_in_range_ind) 1])))';
%         end
        south_ind = find(lat_in_range_midpts(ceil(size(lon_in_range_midpts,1)/2),:) <= lat_transect);
        north_ind = find(lat_in_range_midpts(ceil(size(lon_in_range_midpts,1)/2),:) > lat_transect);
        if ((isempty(south_ind) == 1) || (isempty(north_ind) == 1))
            south_ind = find(lat_in_range_midpts(round(2/grid_res),:) <= lat_transect);
            [~,nearest_south_ind] = min(abs(lat_in_range_midpts(round(2/grid_res),south_ind) - lat_transect));
            nearest_south_ind = south_ind(nearest_south_ind);
            north_ind = find(lat_in_range_midpts(round(2/grid_res),:) > lat_transect);
            [~,nearest_north_ind] = min(abs(lat_in_range_midpts(round(2/grid_res),north_ind) - lat_transect));
            nearest_north_ind = north_ind(nearest_north_ind);
            weighted_south = (lat_in_range_midpts(round(2/grid_res),nearest_north_ind) - lat_transect)/diff(lat_in_range_midpts(round(2/grid_res),[nearest_south_ind nearest_north_ind]));
        else
            [~,nearest_south_ind] = min(abs(lat_in_range_midpts(ceil(size(lon_in_range_midpts,1)/2),south_ind) - lat_transect));
            nearest_south_ind = south_ind(nearest_south_ind);
            [~,nearest_north_ind] = min(abs(lat_in_range_midpts(ceil(size(lon_in_range_midpts,1)/2),north_ind) - lat_transect));
            nearest_north_ind = north_ind(nearest_north_ind);
            weighted_south = (lat_in_range_midpts(ceil(size(lon_in_range_midpts,1)/2),nearest_north_ind) - lat_transect)/diff(lat_in_range_midpts(ceil(size(lon_in_range_midpts,1)/2),[nearest_south_ind nearest_north_ind]));
        end            
        weighted_north = 1 - weighted_south;
        EKE_at_lat = squeeze((weighted_south*EKE_in_range(:,nearest_south_ind,:)) + (weighted_north*EKE_in_range(:,nearest_north_ind,:)));
        good_coords_ind = find(isnan(lon_in_range_midpts(:,ceil(size(lon_in_range_midpts,2)/2))) == 0);
        EKE_along_transect(1,lon_pts_in_range_ind,:) = reshape(interp1(lon_in_range_midpts(good_coords_ind,ceil(size(lon_in_range_midpts,2)/2)) + (round((min(lon_in_range_all(lon_pts_in_range_ind)) - min(lon_in_range_midpts(good_coords_ind,ceil(size(lon_in_range_midpts,2)/2))))/360)*360),EKE_at_lat(good_coords_ind,:),lon_in_range_all(lon_pts_in_range_ind)),[1 length(lon_pts_in_range_ind) length(time_datenum_in_range)]);
    else
        near_lat_ind = find((latitude - latsouth >= 0) & (latitude - latnorth <= 0));
        
        lon_width_subsets = 10;
        lon_min_subsets = (lon_min:lon_width_subsets:lon_max)';
        lon_max_subsets = lon_min_subsets + lon_width_subsets;
        lon_max_subsets(lon_max_subsets > lon_max) = lon_max;
        longwest_subsets = lon_min_subsets - 1;
        longeast_subsets = lon_max_subsets + 1;
        
        EKE_along_transect = NaN([1 length(lon_in_range_all) length(time_datenum_in_range)]);
        for subset_ind = 1:length(lon_min_subsets)
            curr_lon_min = lon_min_subsets(subset_ind);
            curr_lon_max = lon_max_subsets(subset_ind);
            curr_longwest = longwest_subsets(subset_ind);
            curr_longeast = longeast_subsets(subset_ind);
            
            curr_lon_range_ind = find((mod(longitude(near_lat_ind) - curr_longwest + 180,360) - 180 >= 0) & (mod(longitude(near_lat_ind) - curr_longeast + 180,360) - 180 <= 0));
            curr_region_ind = near_lat_ind(curr_lon_range_ind);
            if isempty(curr_region_ind) == 1
                continue
            end
            curr_region_i_ind = unique(mod(curr_region_ind - 1,size(longitude,1)) + 1);
            curr_region_j_ind = unique(ceil(curr_region_ind/size(latitude,1)));
            
            gap_ind = find(diff(curr_region_i_ind) > 5);
            gap_ind_i_connect = find(diff([curr_region_i_ind; (curr_region_i_ind + size(longitude,1))]) > 5);
            curr_region_i_ind_halfdomain = min([curr_region_i_ind (size(longitude,1) + 1 - curr_region_i_ind)],[],2);
            gap_ind_j_connect = find(diff(sort(curr_region_i_ind_halfdomain,'ascend')) > 5);
            if ((lat_transect > 60) && (((length(gap_ind_j_connect) < length(gap_ind)) && (max(curr_region_i_ind) - min(curr_region_i_ind) < size(longitude,1) - 10)) || (max(curr_region_j_ind) >= size(latitude,2) - 1)))
                if max(curr_region_i_ind) - min(curr_region_i_ind) > size(longitude,1) - 10
                    in_lon_range_ind = [((size(longitude,1) + 1 - max(curr_region_i_ind_halfdomain)):1:size(longitude,1))'; (1:1:max(curr_region_i_ind_halfdomain))'];
                elseif max(curr_region_i_ind_halfdomain) > (size(longitude,1)/2) - 1
                    in_lon_range_ind = (min(curr_region_i_ind_halfdomain):1:(size(longitude,1) + 1 - min(curr_region_i_ind_halfdomain)))';
                else
                    in_lon_range_ind = (min(curr_region_i_ind_halfdomain):1:max(curr_region_i_ind_halfdomain))';
                end
                in_lon_range_ind_2 = size(longitude,1) + 1 - in_lon_range_ind;
                in_lat_range_ind = (min(curr_region_j_ind):1:size(latitude,2))';
                in_lat_range_ind_2 = flip(in_lat_range_ind);
                
                lon_in_range = [longitude(in_lon_range_ind,in_lat_range_ind) longitude(in_lon_range_ind_2,in_lat_range_ind_2)];
                lat_in_range = [latitude(in_lon_range_ind,in_lat_range_ind) latitude(in_lon_range_ind_2,in_lat_range_ind_2)];
                time_datenum_in_range_all = time_datenum(in_time_range_ind);
                time_bound_datenum_in_range_all = time_bound(in_time_range_ind,:);
                file_num_in_range_vec_all = file_num_vec(in_time_range_ind);
                in_file_ind_in_range_vec_all = in_file_ind_vec(in_time_range_ind);
                
                
%                 in_time_range_ind = in_time_range_ind_all;
%                 time_datenum_in_range = time_datenum_in_range_all;
                file_num_in_range_vec = file_num_in_range_vec_all;
                in_file_ind_in_range_vec = in_file_ind_in_range_vec_all;
                
                
                size_array = [length(in_lon_range_ind) length([in_lat_range_ind; in_lat_range_ind_2]) 1 length(in_time_range_ind)];
%                 delta_t = mean(diff(time_datenum_in_range));
                
                
                SSH_in_range = NaN(size_array([1 2 4]));
                for file_num_ind = 1:(max(file_num_in_range_vec) - min(file_num_in_range_vec) + 1)
                    curr_model_yr = time_range_start(1) + (file_num_ind - 1) - year_offset;
                    if curr_model_yr < 10
                        curr_yearstr = ['000',num2str(curr_model_yr)];
                    else
                        curr_yearstr = ['00',num2str(curr_model_yr)];
                    end
                    curr_filename = [gridded_SSH_fileform_1,curr_yearstr,gridded_SSH_fileform_2];
            
                    curr_time_subset_ind = find(file_num_in_range_vec == file_num_ind);
                    if isempty(curr_time_subset_ind) == 1
                        continue
                    end
                    curr_in_time_range_ind = in_file_ind_in_range_vec(curr_time_subset_ind);
                    
                    if min(diff(in_lon_range_ind)) < 0
                        gap_ind = find(diff(in_lon_range_ind) < 0);
                        
                        start_vec = [in_lon_range_ind(1) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
                        count_vec = [(diff(in_lon_range_ind([1 gap_ind])) + 1) (max(in_lat_range_ind) - min(in_lon_range_ind) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                        curr_SSH_in_range = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
                        start_vec_2 = [in_lon_range_ind(gap_ind + 1) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
                        count_vec_2 = [(diff(in_lon_range_ind([(gap_ind + 1) length(in_lon_range_ind)])) + 1) (max(in_lat_range_ind) - min(in_lon_range_ind) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                        curr_SSH_in_range_2 = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec_2([1 2 4]),count_vec_2([1 2 4]))));
                        curr_SSH_in_range = [curr_SSH_in_range; curr_SSH_in_range_2];
                        
                        start_vec = [in_lon_range_ind_2(length(in_lon_range_ind_2)) min(in_lat_range_ind_2) 1 min(curr_in_time_range_ind)];
                        count_vec = [(diff(in_lon_range_ind_2([length(in_lon_range_ind_2) (gap_ind + 1)])) + 1) (max(in_lat_range_ind_2) - min(in_lat_range_ind_2) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                        curr_SSH_in_range_3 = 0.01*flip(flip(squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4])))),1),2);
                        start_vec_2 = [in_lon_range_ind_2(gap_ind) min(in_lat_range_ind_2) 1 min(curr_in_time_range_ind)];
                        count_vec_2 = [(diff(in_lon_range_ind([gap_ind 1])) + 1) (max(in_lat_range_ind_2) - min(in_lat_range_ind_2) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                        curr_SSH_in_range_4 = 0.01*flip(flip(squeeze(double(ncread(curr_filename,'SSH',start_vec_2([1 2 4]),count_vec_2([1 2 4])))),1),2);
                        curr_SSH_in_range = [curr_SSH_in_range [curr_SSH_in_range_3; curr_SSH_in_range_4]];
                    else
                        start_vec = [min(in_lon_range_ind) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
                        count_vec = [(max(in_lon_range_ind) - min(in_lon_range_ind) + 1) (max(in_lat_range_ind) - min(in_lat_range_ind) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                        curr_SSH_in_range = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
                        
                        start_vec = [min(in_lon_range_ind_2) min(in_lat_range_ind_2) 1 min(curr_in_time_range_ind)];
                        count_vec = [(max(in_lon_range_ind_2) - min(in_lon_range_ind_2) + 1) (max(in_lat_range_ind_2) - min(in_lat_range_ind_2) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                        curr_SSH_in_range_3 = 0.01*flip(flip(squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4])))),1),2);
                        curr_SSH_in_range = [curr_SSH_in_range curr_SSH_in_range_3];
                    end
                    curr_SSH_in_range(abs(curr_SSH_in_range) > 1e15) = 0;
                    SSH_in_range(:,:,curr_time_subset_ind) = curr_SSH_in_range(:,:,curr_in_time_range_ind - min(curr_in_time_range_ind) + 1);
                    
                end
                
            else
                if max(curr_region_i_ind) - min(curr_region_i_ind) > size(longitude,1) - 10
                    in_lon_range_ind = [((size(longitude,1) + 1 - max(curr_region_i_ind)):1:size(longitude,1))'; (1:1:max(curr_region_i_ind))'];
                else
                    in_lon_range_ind = (min(curr_region_i_ind):1:max(curr_region_i_ind))';
                end
                in_lat_range_ind = (min(curr_region_j_ind):1:max(curr_region_j_ind))';
                
                lon_in_range = longitude(in_lon_range_ind,in_lat_range_ind);
                lat_in_range = latitude(in_lon_range_ind,in_lat_range_ind);
                time_datenum_in_range_all = time_datenum(in_time_range_ind);
                time_bound_datenum_in_range_all = time_bound(in_time_range_ind,:);
                file_num_in_range_vec_all = file_num_vec(in_time_range_ind);
                in_file_ind_in_range_vec_all = in_file_ind_vec(in_time_range_ind);
                
                
%                 in_time_range_ind = in_time_range_ind_all;
%                 time_datenum_in_range = time_datenum_in_range_all;
                file_num_in_range_vec = file_num_in_range_vec_all;
                in_file_ind_in_range_vec = in_file_ind_in_range_vec_all;
                
                
                size_array = [length(in_lon_range_ind) length(in_lat_range_ind) 1 length(in_time_range_ind)];
%                 delta_t = mean(diff(time_datenum_in_range));
                
                
                SSH_in_range = NaN(size_array([1 2 4]));
                for file_num_ind = 1:(max(file_num_in_range_vec) - min(file_num_in_range_vec) + 1)
                    curr_model_yr = time_range_start(1) + (file_num_ind - 1) - year_offset;
                    if curr_model_yr < 10
                        curr_yearstr = ['000',num2str(curr_model_yr)];
                    else
                        curr_yearstr = ['00',num2str(curr_model_yr)];
                    end
                    curr_filename = [gridded_SSH_fileform_1,curr_yearstr,gridded_SSH_fileform_2];
            
                    curr_time_subset_ind = find(file_num_in_range_vec == file_num_ind);
                    if isempty(curr_time_subset_ind) == 1
                        continue
                    end
                    curr_in_time_range_ind = in_file_ind_in_range_vec(curr_time_subset_ind);
                    
                    if min(diff(in_lon_range_ind)) < 0
                        gap_ind = find(diff(in_lon_range_ind) < 0);
                        
                        start_vec = [in_lon_range_ind(1) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
                        count_vec = [(diff(in_lon_range_ind([1 gap_ind])) + 1) (max(in_lat_range_ind) - min(in_lon_range_ind) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                        curr_SSH_in_range = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
                        start_vec_2 = [in_lon_range_ind(gap_ind + 1) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
                        count_vec_2 = [(diff(in_lon_range_ind([(gap_ind + 1) length(in_lon_range_ind)])) + 1) (max(in_lat_range_ind) - min(in_lon_range_ind) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                        curr_SSH_in_range_2 = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec_2([1 2 4]),count_vec_2([1 2 4]))));
                        curr_SSH_in_range = [curr_SSH_in_range; curr_SSH_in_range_2];
                    else
                        start_vec = [min(in_lon_range_ind) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
                        count_vec = [(max(in_lon_range_ind) - min(in_lon_range_ind) + 1) (max(in_lat_range_ind) - min(in_lat_range_ind) + 1) 1 (max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1)];
                        curr_SSH_in_range = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
                    end
                    
                    curr_SSH_in_range(abs(curr_SSH_in_range) > 1e15) = 0;
                    SSH_in_range(:,:,curr_time_subset_ind) = curr_SSH_in_range(:,:,curr_in_time_range_ind - min(curr_in_time_range_ind) + 1);
                    
                end
                
            end
            
            lon_in_range(isnan(lon_in_range) == 1) = 1e4;
            diff_lon_in_range = diff([lon_in_range; lon_in_range(1,:)],1,1); 
            diff_lon_in_range((abs(diff_lon_in_range) > 180) & (abs(diff_lon_in_range) < 5e3)) = mod(diff_lon_in_range((abs(diff_lon_in_range) > 180) & (abs(diff_lon_in_range) < 5e3)) + 180,360) - 180;
            lon_in_range_aug_x = repmat(lon_in_range(1,:),[(size_array(1) + 1) 1]) + [zeros([1 size_array(2)]); cumsum(diff_lon_in_range,1,'forward')];
            lon_in_range_aug_x(lon_in_range_aug_x > 5e3) = NaN;
            lon_in_range_aug_x(lon_in_range_aug_x < min(lon_in_range_aug_x(1,:))) = lon_in_range_aug_x(lon_in_range_aug_x < min(lon_in_range_aug_x(1,:))) + (360*ceil((min(lon_in_range_aug_x(1,:)) - lon_in_range_aug_x(lon_in_range_aug_x < min(lon_in_range_aug_x(1,:))))/360));
            lon_in_range_aug_x(size_array(1) + 1,lon_in_range_aug_x(size_array(1) + 1,:) < min(lon_in_range_aug_x(1,:)) + 180) = lon_in_range_aug_x(size_array(1) + 1,lon_in_range_aug_x(size_array(1) + 1,:) < min(lon_in_range_aug_x(1,:)) + 180) + 360;
            lon_in_range_aug_x(lon_in_range_aug_x > max(lon_in_range_aug_x(size_array(1) + 1,:))) = lon_in_range_aug_x(lon_in_range_aug_x > max(lon_in_range_aug_x(size_array(1) + 1,:))) + (360*floor((max(lon_in_range_aug_x(size_array(1) + 1,:)) - lon_in_range_aug_x(lon_in_range_aug_x > max(lon_in_range_aug_x(size_array(1) + 1,:))))/360));
            lon_in_range = lon_in_range_aug_x(1:size_array(1),:);
            lon_in_range(abs(lon_in_range - 1e4) < 1e-5) = NaN;
            lon_in_range_midpts = lon_in_range_aug_x(1:size_array(1),:) + (diff(lon_in_range_aug_x,1,1)/2);
            lon_in_range_midpts = [(lon_in_range_midpts(:,2:size_array(2)) - (diff(lon_in_range_midpts,1,2)/2)) NaN([size_array(1) 1])];
            lat_in_range_aug_x = [lat_in_range; lat_in_range(1,:)];
            lat_in_range_midpts = lat_in_range_aug_x(1:size_array(1),:) + (diff(lat_in_range_aug_x,1,1)/2);
            lat_in_range_midpts = [(lat_in_range_midpts(:,2:size_array(2)) - (diff(lat_in_range_midpts,1,2)/2)) NaN([size_array(1) 1])];
            dx = 111100*(abs((cosd(lat_in_range_midpts).*diff([(lon_in_range_aug_x(:,2:size_array(2)) - (diff(lon_in_range_aug_x,1,2)/2)) NaN([(size_array(1) + 1) 1])],1,1)) + (1i*diff([(lat_in_range_aug_x(:,2:size_array(2)) - (diff(lat_in_range_aug_x,1,2)/2)) NaN([(size_array(1) + 1) 1])],1,1))));
            dy = 111100*(abs((cosd(lat_in_range_midpts).*diff([(lon_in_range_aug_x(1:size_array(1),:) + (diff(lon_in_range_aug_x,1,1)/2)) NaN([size_array(1) 1])],1,2)) + (1i*[diff(lat_in_range_aug_x(1:size_array(1),:) + (diff(lat_in_range_aug_x,1,1)/2),1,2) NaN([size_array(1) 1])])));
            
            
            size_array = size(SSH_in_range);
            
            nan_mask = ones(size(SSH_in_range));
            nan_mask(isnan(SSH_in_range) == 1) = 0;
            SSH_in_range_nans_zeroed = SSH_in_range;
            SSH_in_range_nans_zeroed(isnan(SSH_in_range) == 1) = 0;
            sum_nan_mask = sum(nan_mask,3);
            nan_mask_not_enough = ones(size(sum_nan_mask));
            nan_mask_not_enough(sum_nan_mask < (0.9*max(max(sum_nan_mask)))) = 0;
            nan_mask_not_enough_3D = repmat(nan_mask_not_enough,[1 1 size(nan_mask,3)]);
            nan_mask_SSH = nan_mask & nan_mask_not_enough_3D;
            
            clear nan_mask nan_mask_not_enough*
            
            SSH_mean_in_range = sum(nan_mask_SSH.*SSH_in_range_nans_zeroed,3)./sum(nan_mask_SSH,3);
            
            SSH_nomean_in_range = SSH_in_range - repmat(SSH_mean_in_range,[1 1 size(SSH_in_range,3)]);
            
            
            diff_x_SSH_nomean_in_r = diff(SSH_nomean_in_range([(1:1:size_array(1))'; 1],:,:),1,1);
            diff_x_SSH_nomean_in_range = [(diff_x_SSH_nomean_in_r(:,2:size(diff_x_SSH_nomean_in_r,2),:) - (diff(diff_x_SSH_nomean_in_r,1,2)/2)) NaN([size_array(1) 1 size_array(3)])];
            diff_y_SSH_nomean_in_r = diff(SSH_nomean_in_range,1,2);
            diff_y_SSH_nomean_in_range = [(diff_y_SSH_nomean_in_r([(2:1:size(diff_y_SSH_nomean_in_r,1))'; 1],:,:) - (diff(diff_y_SSH_nomean_in_r([(1:1:size(diff_y_SSH_nomean_in_r,1))'; 1],:,:),1,1)/2)) NaN([size_array(1) 1 size_array(3)])];
        
        
            f = 2*((2*pi)/86164)*sind(lat_in_range_midpts);
            
            uvel_geostr_nomean_in_range = -9.81*repmat((f.^(-1)),[1 1 size(SSH_nomean_in_range,3)]).*(diff_y_SSH_nomean_in_range./repmat(dy,[1 1 size(SSH_nomean_in_range,3)]));
            vvel_geostr_nomean_in_range = 9.81*repmat((f.^(-1)),[1 1 size(SSH_nomean_in_range,3)]).*(diff_x_SSH_nomean_in_range./repmat(dx,[1 1 size(SSH_nomean_in_range,3)]));
        
            clear diff_*        
        
            EKE_in_range = (0.5*((uvel_geostr_nomean_in_range.^2) + (vvel_geostr_nomean_in_range.^2)));
        
            clear uvel_geostr* vvel_geostr*
        
            % mask out low-latitude points
            lat_mask = ones(size(lat_in_range_midpts));
            lat_mask(abs(lat_in_range_midpts) < 5) = 0;
            lat_mask_3d = repmat(lat_mask,[1 1 size(EKE_in_range,3)]);
            lat_mask_ind = find(abs(lat_mask_3d) < 1e-5);
        
            EKE_in_range(lat_mask_ind) = NaN;
            
            % interpolate to transect points
            
            EKE_in_range(isnan(EKE_in_range) == 1) = 0;
            
            lon_pts_in_range_ind = find(abs(mod(lon_in_range_all - mean([curr_lon_min curr_lon_max]) + 180,360) - 180) <= ((curr_lon_max - curr_lon_min)/2));
            EKE_snapshot = EKE_in_range(:,:,ceil(size(EKE_in_range,3)/2));
            good_ind = find(isnan(EKE_snapshot) == 0);
            for curr_lon_ind = 1:length(lon_pts_in_range_ind)
                [sorted_dist,sorted_dist_ind] = sort(abs((cosd(lat_transect)*(mod(lon_in_range_midpts(good_ind) - lon_in_range_all(lon_pts_in_range_ind(curr_lon_ind)) + 180,360) - 180)) + (1i*(lat_in_range_midpts(good_ind) - lat_transect))),'ascend');
                dist_radius = 2*grid_res;
                sorted_dist = sorted_dist(1:4);
                sorted_dist_ind = good_ind(sorted_dist_ind(1:4));
                sorted_dist_ind = sorted_dist_ind(sorted_dist < dist_radius);
                sorted_dist = sorted_dist(sorted_dist < dist_radius);
                weighted_dist = dist_radius - sorted_dist;
                weighted_dist = weighted_dist/sum(weighted_dist);
                for t = 1:length(time_datenum_in_range)
                    curr_EKE_in_range = EKE_in_range(:,:,t);
                    EKE_along_transect(1,lon_pts_in_range_ind(curr_lon_ind),t) = sum(weighted_dist.*curr_EKE_in_range(sorted_dist_ind));
                end
            end
            
        end        
    end
    
    nan_mask_EKE = ones(size(EKE_along_transect));
    nan_mask_EKE(isnan(EKE_along_transect) == 1) = 0;
    EKE_along_transect_nans_zeroed = EKE_along_transect;
    EKE_along_transect_nans_zeroed(isnan(EKE_along_transect) == 1) = 0;
    
    EKE_mean_along_transect = sum(nan_mask_EKE.*EKE_along_transect_nans_zeroed,3)./sum(nan_mask_EKE,3);
    
    
    % load satellite SSH and compute EKE
    longwest_satobs = mod(longwest - (min(lon_satobs) - 1e-10),360) + (min(lon_satobs) - 1e-10);
    longeast_satobs = mod(longeast - (min(lon_satobs) - 1e-10),360) + (min(lon_satobs) - 1e-10);
%     if ((longeast_satobs - longwest_satobs < 0) || (mod(longeast_satobs - longwest_satobs,360) < 2 + 1e-10))
%         in_lon_range_sat_ind = find((lon_satobs >= longwest_satobs) || (lon_satobs <= longeast_satobs));
%     else
%         in_lon_range_sat_ind = find((lon_satobs >= longwest_satobs) && (lon_satobs <= longeast_satobs));
%     end
    in_lat_range_sat_ind = find((lat_satobs >= latsouth) & (lat_satobs <= latnorth));
    lat_in_range_satobs = lat_satobs(in_lat_range_sat_ind);
    [~,start_lon_sat_ind] = min(abs(mod(lon_satobs - (longwest_satobs + (0.4999*grid_res_satobs)) + 180,360) - 180));
    [~,end_lon_sat_ind] = min(abs(mod(lon_satobs - (longeast_satobs - (0.4999*grid_res_satobs)) + 180,360) - 180));
    if end_lon_sat_ind - start_lon_sat_ind <= (2.1/grid_res_satobs)
        in_lon_range_sat_ind = [(start_lon_sat_ind:1:length(lon_satobs))'; (1:1:end_lon_sat_ind)'];
        lon_in_range_satobs = lon_satobs(in_lon_range_sat_ind);
        start_vec = [in_lon_range_sat_ind(1) in_lat_range_sat_ind(1) 1];
        count_vec = [(length(lon_satobs) - in_lon_range_sat_ind(1) + 1) length(in_lat_range_sat_ind) length(time_datenum_satobs)];
        SSH_satobs = ncread(nc_file_name_satobs,'SSH',start_vec,count_vec);
        start_vec_2 = [1 in_lat_range_sat_ind(1) 1];
        count_vec_2 = [in_lon_range_sat_ind(length(in_lon_range_sat_ind)) length(in_lat_range_sat_ind) length(time_datenum_satobs)];
        SSH_satobs_2 = ncread(nc_file_name_satobs,'SSH',start_vec_2,count_vec_2);
        SSH_satobs = [SSH_satobs; SSH_satobs_2];
    else
        in_lon_range_sat_ind = (start_lon_sat_ind:1:end_lon_sat_ind)';
        lon_in_range_satobs = lon_satobs(in_lon_range_sat_ind);
        start_vec = [in_lon_range_sat_ind(1) in_lat_range_sat_ind(1) 1];
        count_vec = [length(in_lon_range_sat_ind) length(in_lat_range_sat_ind) length(time_datenum_satobs)];
        SSH_satobs = ncread(nc_file_name_satobs,'SSH',start_vec,count_vec);
    end
    
    diff_lon_in_range_satobs = diff(lon_in_range_satobs);
    diff_lon_in_range_satobs((diff_lon_in_range_satobs < 0) | (diff_lon_in_range_satobs >= 360)) = mod(diff_lon_in_range_satobs((diff_lon_in_range_satobs < 0) | (diff_lon_in_range_satobs >= 360)),360);
    lon_in_range_satobs = repmat(lon_in_range_satobs(1),[length(lon_in_range_satobs) 1]) + [0; cumsum(diff_lon_in_range_satobs)];
    

    size_array_satobs = size(SSH_satobs);
    
    nan_mask = ones(size(SSH_satobs));
    nan_mask(isnan(SSH_satobs) == 1) = 0;
    sum_nan_mask = sum(nan_mask,3);
    not_enough_ind = sum_nan_mask < (0.9*max(max(sum_nan_mask)));
    
%     nan_mask_array_2d = (1e-10)*ones([size_array(1) size_array(2)]);
%     nan_mask_array_2d(not_enough_ind) = -1;
    
    
    
    % compute EKE associated with each SSH component
    
    nan_mask = ones(size(SSH_satobs));
    nan_mask(isnan(SSH_satobs) == 1) = 0;
    SSH_satobs_nans_zeroed = SSH_satobs;
    SSH_satobs_nans_zeroed(isnan(SSH_satobs) == 1) = 0;
    sum_nan_mask = sum(nan_mask,3);
    nan_mask_not_enough = ones(size(sum_nan_mask));
    nan_mask_not_enough(sum_nan_mask < (0.9*max(max(sum_nan_mask)))) = 0;
    nan_mask_not_enough_3D = repmat(nan_mask_not_enough,[1 1 size(nan_mask,3)]);
    nan_mask_SSH = nan_mask & nan_mask_not_enough_3D;
    
    clear nan_mask nan_mask_not_enough*
    
    SSH_mean_satobs = sum(nan_mask_SSH.*SSH_satobs_nans_zeroed,3)./sum(nan_mask_SSH,3);
    
    SSH_nomean_satobs = SSH_satobs - repmat(SSH_mean_satobs,[1 1 size(SSH_satobs,3)]);
    
    clear SSH_*_nans_zeroed
    
    
    % compute EKE from geostrophic gradients
    
    lon_satobs_midpts = lon_in_range_satobs(2:length(lon_in_range_satobs)) - (diff(lon_in_range_satobs)/2);
    lat_satobs_midpts = lat_in_range_satobs(2:length(lat_in_range_satobs)) - (diff(lat_in_range_satobs)/2);
    dx = 111100*(repmat(cosd(lat_satobs_midpts'),[(size_array_satobs(1) - 1) 1]).*repmat(diff(lon_in_range_satobs),[1 (size_array_satobs(2) - 1)]));
    dy = 111100*repmat(diff(lat_in_range_satobs'),[(size_array_satobs(1) - 1) 1]);
    
    diff_x_SSH_nomean_in_r = diff(SSH_nomean_satobs,1,1);
    diff_x_SSH_nomean_satobs = diff_x_SSH_nomean_in_r(:,2:size(diff_x_SSH_nomean_in_r,2),:) - (diff(diff_x_SSH_nomean_in_r,1,2)/2);
    diff_y_SSH_nomean_in_r = diff(SSH_nomean_satobs,1,2);
    diff_y_SSH_nomean_satobs = diff_y_SSH_nomean_in_r(2:size(diff_y_SSH_nomean_in_r,1),:,:) - (diff(diff_y_SSH_nomean_in_r,1,1)/2);
    
    
    f = repmat(2*((2*pi)/86164)*sind(lat_satobs_midpts'),[length(lon_satobs_midpts) 1]);
        
    uvel_geostr_nomean_satobs = -9.81*repmat((f.^(-1)),[1 1 size(SSH_nomean_satobs,3)]).*(diff_y_SSH_nomean_satobs./repmat(dy,[1 1 size(SSH_nomean_satobs,3)]));
    vvel_geostr_nomean_satobs = 9.81*repmat((f.^(-1)),[1 1 size(SSH_nomean_satobs,3)]).*(diff_x_SSH_nomean_satobs./repmat(dx,[1 1 size(SSH_nomean_satobs,3)]));
    
    clear diff_*
    
    
    EKE_satobs = (0.5*((uvel_geostr_nomean_satobs.^2) + (vvel_geostr_nomean_satobs.^2)));
    
    clear uvel_geostr* vvel_geostr*
    
    % mask out low-latitude points
    lat_mask = ones(size(lat_satobs_midpts));
    lat_mask(abs(lat_satobs_midpts) < 5) = 0;
    lat_mask_3d = repmat(lat_mask',[size(EKE_satobs,1) 1 size(EKE_satobs,3)]);
    lat_mask_ind = find(abs(lat_mask_3d) < 1e-5);
    
    EKE_satobs(lat_mask_ind) = NaN;
    
    nan_mask_EKE = ones(size(EKE_satobs));
    nan_mask_EKE(isnan(EKE_satobs) == 1) = 0;
    EKE_satobs_nans_zeroed = EKE_satobs;
    EKE_satobs_nans_zeroed(isnan(EKE_satobs) == 1) = 0;
    
    EKE_mean_satobs = sum(nan_mask_EKE.*EKE_satobs_nans_zeroed,3)./sum(nan_mask_EKE,3);
    EKE_mean_satobs(isnan(EKE_mean_satobs) == 1) = 0;
    
    diff_lon_satobs_midpts = diff(lon_satobs_midpts);
    diff_lon_satobs_midpts(diff_lon_satobs_midpts < 0) = mod(diff_lon_satobs_midpts(diff_lon_satobs_midpts < 0),360);
    lon_satobs_midpts_adj = lon_satobs_midpts(1) + [0; cumsum(diff_lon_satobs_midpts)] + (360*round((lon_in_range_all(1) - lon_satobs_midpts(1))/360));
    EKE_mean_satobs_model_grid = interp2_fft(lon_satobs_midpts_adj,lat_satobs_midpts,EKE_mean_satobs,lon_in_range_all',lat_transect*ones([1 length(lon_in_range_vel_all)]));
    EKE_mean_satobs_model_grid(isnan(EKE_mean_along_transect) == 1) = NaN;
    
    
    % apply spatial smoothing filter
    half_power_adj = 1;
%     EKE_along_transect_smoothed_array = NaN([1 length(lon_in_range_all) length(time_datenum_in_range)]);
%     EKE_mean_along_transect_smoothed_array = NaN([1 length(lon_in_range_all)]);
%     EKE_mean_satobs_model_grid_smoothed_array = NaN([1 length(lon_in_range_all)]);
    spatial_scale_separation = spatial_scale_separation_vec(transect_ind);
    if sum(diff(lon_in_range_all)) >= 359
        curr_size_array = size(EKE_along_transect);
        buffer_size = round(spatial_scale_separation/(mean(diff(lon_in_range_all))));
        buffered_ind = [(curr_size_array(2) + ((-buffer_size + 1):1:0)'); (1:1:curr_size_array(2))'; (1:1:buffer_size)'];
        unbuffered_ind = buffer_size + (1:1:curr_size_array(2))';
        [curr_array,~,~] = bandpass_err_fcn(EKE_along_transect(1,buffered_ind,:),2,grid_res,1/(4*grid_res*length(lon_in_range_all))/half_power_adj,(1/(2*spatial_scale_separation))*half_power_adj,2,1,1,0,1);
        EKE_along_transect_smoothed_array(transect_ind,:,:) = curr_array(:,unbuffered_ind,:);
        [curr_array,~,~] = bandpass_err_fcn(EKE_mean_along_transect(1,buffered_ind),2,grid_res,1/(4*grid_res*length(lon_in_range_all))/half_power_adj,(1/(2*spatial_scale_separation))*half_power_adj,2,1,1,0,1);
        EKE_mean_along_transect_smoothed_array(transect_ind,:) = curr_array(:,unbuffered_ind);
        [curr_array,~,~] = bandpass_err_fcn(EKE_mean_satobs_model_grid(1,buffered_ind),2,grid_res,1/(4*grid_res*length(lon_in_range_all))/half_power_adj,(1/(2*spatial_scale_separation))*half_power_adj,2,1,1,0,1);
        EKE_mean_satobs_model_grid_smoothed_array(transect_ind,:) = curr_array(:,unbuffered_ind);
    else
        [curr_array,~,~] = bandpass_err_fcn(EKE_along_transect(1,:,:),2,grid_res,1/(4*grid_res*length(lon_in_range_all))/half_power_adj,(1/(2*spatial_scale_separation))*half_power_adj,2,1,1,0,1);
        EKE_along_transect_smoothed_array(transect_ind,:,:) = curr_array;
        [curr_array,~,~] = bandpass_err_fcn(EKE_mean_along_transect(1,:),2,grid_res,1/(4*grid_res*length(lon_in_range_all))/half_power_adj,(1/(2*spatial_scale_separation))*half_power_adj,2,1,1,0,1);
        EKE_mean_along_transect_smoothed_array(transect_ind,:) = curr_array;
        [curr_array,~,~] = bandpass_err_fcn(EKE_mean_satobs_model_grid(1,:),2,grid_res,1/(4*grid_res*length(lon_in_range_all))/half_power_adj,(1/(2*spatial_scale_separation))*half_power_adj,2,1,1,0,1);
        EKE_mean_satobs_model_grid_smoothed_array(transect_ind,:) = curr_array;
    end
    
%     EKE_along_transect_smoothed_array(EKE_along_transect_smoothed_array < 0) = 0;
    EKE_mean_along_transect_smoothed_array(transect_ind,EKE_mean_along_transect_smoothed_array(transect_ind,:) < 0) = 0;
    EKE_mean_satobs_model_grid_smoothed_array(transect_ind,EKE_mean_satobs_model_grid_smoothed_array(transect_ind,:) < 0) = 0;
    
    
    % "clip" distributions of spatially-smoothed EKE
    
    EKE_along_transect_smoothed_sorted = sort(EKE_along_transect_smoothed_array(transect_ind,:,:),3,'ascend');
    curr_nan_mask = ones(size(EKE_along_transect_smoothed_sorted));
    curr_nan_mask((isnan(EKE_along_transect_smoothed_sorted) == 1) | (abs(EKE_along_transect_smoothed_sorted) < 1e-10)) = 0;
    EKE_along_transect_smoothed_sorted(abs(curr_nan_mask) < 1e-5) = 0;
    n_times_at_loc = sum(curr_nan_mask,3);
    
    EKE_lower_avg = NaN(size(EKE_along_transect_smoothed_sorted));
    EKE_higher_avg = NaN(size(EKE_along_transect_smoothed_sorted));
    for t = 1:max(n_times_at_loc)
        EKE_lower_avg(1,:,t) = sum(curr_nan_mask(1,:,1:t).*EKE_along_transect_smoothed_sorted(1,:,1:t),3)./(sum(curr_nan_mask(1,:,1:t),3));
        EKE_higher_avg(1,:,t) = sum(curr_nan_mask(1,:,t:max(n_times_at_loc)).*EKE_along_transect_smoothed_sorted(1,:,t:max(n_times_at_loc)),3)./(sum(curr_nan_mask(1,:,t:max(n_times_at_loc)),3));
    end

    for loc_ind = 1:size(EKE_along_transect_smoothed_sorted,2)
        if n_times_at_loc(loc_ind) < 0.5*max(n_times_at_loc)
            continue
        end
        if abs(EKE_mean_satobs_model_grid_smoothed_array(transect_ind,loc_ind) - EKE_mean_along_transect_smoothed_array(transect_ind,loc_ind)) < 1e-10
            EKE_low_threshold_adj(transect_ind,loc_ind) = 0;
            EKE_high_threshold_adj(transect_ind,loc_ind) = 1e2;
        elseif EKE_mean_satobs_model_grid_smoothed_array(transect_ind,loc_ind) - EKE_mean_along_transect_smoothed_array(transect_ind,loc_ind) > 0
            [~,unique_ind] = unique(squeeze(EKE_higher_avg(1,loc_ind,:)));
            EKE_low_threshold_adj(transect_ind,loc_ind) = interp1([(EKE_higher_avg(1,loc_ind,unique_ind(1)) - 1e-3); squeeze(EKE_higher_avg(1,loc_ind,unique_ind))],[(EKE_along_transect_smoothed_sorted(1,loc_ind,unique_ind(1)) - 1e-10); squeeze(EKE_along_transect_smoothed_sorted(1,loc_ind,unique_ind))],EKE_mean_satobs_model_grid_smoothed_array(transect_ind,loc_ind));
            if isnan(EKE_low_threshold_adj(transect_ind,loc_ind)) == 1
                EKE_low_threshold_adj(transect_ind,loc_ind) = 1e2;
            end
            EKE_high_threshold_adj(transect_ind,loc_ind) = 1e2;
        else
            EKE_low_threshold_adj(transect_ind,loc_ind) = 0;
            [~,unique_ind] = unique(squeeze(EKE_lower_avg(1,loc_ind,:)));
            EKE_high_threshold_adj(transect_ind,loc_ind) = interp1([squeeze(EKE_lower_avg(1,loc_ind,unique_ind)); (EKE_lower_avg(1,loc_ind,unique_ind(length(unique_ind))) + 1e-3)],[squeeze(EKE_along_transect_smoothed_sorted(1,loc_ind,unique_ind)); (EKE_along_transect_smoothed_sorted(1,loc_ind,unique_ind(length(unique_ind))) + 1e-10)],EKE_mean_satobs_model_grid_smoothed_array(transect_ind,loc_ind));
            if isnan(EKE_high_threshold_adj(transect_ind,loc_ind)) == 1
                EKE_high_threshold_adj(transect_ind,loc_ind) = 0;
            end
        end
        n_days_unmasked_along_transect(transect_ind,loc_ind) = delta_t*length(find((EKE_along_transect_smoothed_sorted(1,loc_ind,:) >= EKE_low_threshold_adj(transect_ind,loc_ind)) & (EKE_along_transect_smoothed_sorted(1,loc_ind,:) <= EKE_high_threshold_adj(transect_ind,loc_ind))));
        
        EKE_low_threshold_adj_alt(transect_ind,loc_ind) = EKE_low_threshold_adj(transect_ind,loc_ind);
        EKE_high_threshold_adj_alt(transect_ind,loc_ind) = EKE_high_threshold_adj(transect_ind,loc_ind);
        
        % provide alternate threshold if unmasked length of time is below min_days_alt
        if n_days_unmasked_along_transect(transect_ind,loc_ind) < min_days_alt
            if EKE_mean_satobs_model_grid_smoothed_array(transect_ind,loc_ind) - EKE_mean_along_transect_smoothed_array(transect_ind,loc_ind) > 0
                EKE_low_threshold_adj_alt(transect_ind,loc_ind) = interp1((1:1:n_times_at_loc(loc_ind))',squeeze(EKE_along_transect_smoothed_sorted(1,loc_ind,abs(EKE_along_transect_smoothed_sorted(1,loc_ind,:)) > 1e-10)),n_times_at_loc(loc_ind) + 0.5 - (min_days_alt/delta_t));
            else
                EKE_high_threshold_adj_alt(transect_ind,loc_ind) = interp1((1:1:n_times_at_loc(loc_ind))',squeeze(EKE_along_transect_smoothed_sorted(1,loc_ind,abs(EKE_along_transect_smoothed_sorted(1,loc_ind,:)) > 1e-10)),0.5 + (min_days_alt/delta_t));
            end
        end
        
    end
    
    disp(['transect_ind = ',num2str(transect_ind),' (out of ',num2str(length(lat_transects_vec)),')'])
end

save(['Mesoscale_EKE_composite_masking_',plot_name,'.mat'],'lat_transects_vec','lon_in_range_all','land_mask_local_plot','EKE_along_transect_smoothed_array','EKE_mean*smoothed_array','EKE_low_threshold_ad*','EKE_high_threshold_ad*','min_days_alt','-v7.3')
% % % % 
% load(['Mesoscale_EKE_composite_masking_',plot_name,'.mat'])

in_adj_range_mask = zeros(size(EKE_along_transect_smoothed_array));
in_adj_range_mask((EKE_along_transect_smoothed_array - repmat(EKE_low_threshold_adj_alt,[1 1 length(time_datenum_in_range)]) > 0) & (EKE_along_transect_smoothed_array - repmat(EKE_high_threshold_adj_alt,[1 1 length(time_datenum_in_range)]) <= 0)) = 1;
EKE_along_transect_smoothed_array_zeronans = EKE_along_transect_smoothed_array;
EKE_along_transect_smoothed_array_zeronans(isnan(EKE_along_transect_smoothed_array) == 1) = 0;

EKE_mean_adj_range_smoothed_array = sum(in_adj_range_mask.*EKE_along_transect_smoothed_array_zeronans,3)./(sum(in_adj_range_mask,3));


lon_bounds_all = [min(lon_bounds_array(:,1)) max(lon_bounds_array(:,2))];
lat_transects_vec_lowbound = [((1.5*lat_transects_vec(1)) + ((-0.5)*lat_transects_vec(2))); (lat_transects_vec(2:length(lat_transects_vec)) - (diff(lat_transects_vec)/2)); ((1.5*lat_transects_vec(length(lat_transects_vec))) + ((-0.5)*lat_transects_vec(length(lat_transects_vec) - 1)))];
lon_in_range_vel_pcolor_plot = [((1.5*lon_in_range_vel_all(1)) + ((-0.5)*lon_in_range_vel_all(2))); (lon_in_range_vel_all(2:length(lon_in_range_vel_all)) - (diff(lon_in_range_vel_all)/2)); ((1.5*lon_in_range_vel_all(length(lon_in_range_vel_all))) + ((-0.5)*lon_in_range_vel_all(length(lon_in_range_vel_all) - 1)))] - (floor((mean(lon_in_range_vel_all) - lon_bounds_all(1))/360)*360);
lon_in_range_pcolor_plot = lon_in_range_vel_pcolor_plot(2:length(lon_in_range_vel_pcolor_plot)) - (diff(lon_in_range_vel_pcolor_plot)/2);

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


cmap = colormap(0.85*bcyr(19));    % colormap for contours
cmap = [[0 0 0.7]; cmap; [0.7 0 0]];

c_levels = [0; (0.05:0.1:1.95)'; 5];
[irregular_clevels_array,~] = irregular_clevels_plot(EKE_mean_adj_range_smoothed_array./EKE_mean_satobs_model_grid_smoothed_array,c_levels);

fig3 = figure(3);
close(figure(1))
colormap(cmap)
[~,cbar] = pcolor_plot_with_landmask(fig3,lon_in_range_pcolor_plot',lat_transects_vec_lowbound,[zeros([length(lat_transects_vec_lowbound) 1]) [irregular_clevels_array; zeros([1 (length(lon_in_range_pcolor_plot) - 1)])]],0:1:(length(c_levels) - 1),cmap,'southoutside',lon_bounds_all(1),lon_bounds_all(2),min(lat_transects_vec_lowbound),max(lat_transects_vec) + 0.5);
shading flat
set(gca,'FontSize',12,'xlim',lon_bounds_all,'xtick',(-360):(10*ceil((diff(lon_bounds_all)/10)/10)):360,'ylim',[min(lat_transects_vec_lowbound) (max(lat_transects_vec) + 0.5)],'ydir','normal','DataAspectRatio',[1 mean(cosd(lat_transects_vec)) 1],'clim',[0 size(cmap,1)])
hold on
line(lon_bounds_all,[0 0],[0 0],'Color','k','LineStyle','-','LineWidth',1)
hold off
xlabel('Longitude')
ylabel('Latitude')
title({'Ratio of model "composite" to satellite EKE,'; ['starting from ',num2str(lon_bounds_all(1)),'^{\circ} longitude, averaged ',datestr_title_start,' to ',datestr_title_end_minus1]; ''},'FontSize',8)
set(cbar,'ticks',0.5:5:20.5,'FontSize',12,'xticklabel',{'0' '0.5' '1' '1.5' '2'})
set(get(cbar,'xlabel'),'String','Ratio of model "composite" to satellite EKE','FontSize',14)

print(fig3,['EKE_ratio_model_composite_satobs_POP_',plot_name,'_from_',num2str(lon_bounds_all(1)),'_to_',num2str(lon_bounds_all(2)),'_lon_',datestr_start,'_',datestr_end_minus1,'_tavg.bmp'],'-dbmp','-r300')
close(fig3)



% average EKE zonally

curr_array = EKE_mean_along_transect_smoothed_array + EKE_mean_satobs_model_grid_smoothed_array + EKE_mean_adj_range_smoothed_array;
EKE_mask_all = ones(size(curr_array));
EKE_mask_all((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;

curr_array = EKE_mean_along_transect_smoothed_array;
curr_array_zeronans = curr_array;
curr_array_zeronans((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
EKE_mean_model_transect_avg = sum(EKE_mask_all.*curr_array_zeronans,2)./(sum(EKE_mask_all,2));
curr_array = EKE_mean_satobs_model_grid_smoothed_array;
curr_array_zeronans = curr_array;
curr_array_zeronans((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
EKE_mean_satobs_transect_avg = sum(EKE_mask_all.*curr_array_zeronans,2)./(sum(EKE_mask_all,2));
curr_array = EKE_mean_adj_range_smoothed_array;
curr_array_zeronans = curr_array;
curr_array_zeronans((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
EKE_mean_adj_range_transect_avg = sum(EKE_mask_all.*curr_array_zeronans,2)./(sum(EKE_mask_all,2));


fig4 = figure(4);
h = plot((1e4)*EKE_mean_model_transect_avg,lat_transects_vec,'r-',(1e4)*EKE_mean_satobs_transect_avg,lat_transects_vec,'b-',(1e4)*EKE_mean_adj_range_transect_avg,lat_transects_vec,'r--');
set(h,'LineWidth',1)
set(h(1),'Color',[0.8 0 0])
set(h(2),'Color',[0.2 0.2 1])
set(h(3),'Color',[0.8 0 0])
set(gca,'xlim',[0 600],'ylim',[min(lat_transects_vec) max(lat_transects_vec)],'DataAspectRatio',[700 (max(lat_transects_vec) - min(lat_transects_vec)) 1],'FontSize',12)
hold on
line([0 600],[0 0],'Color',[0 0 0],'LineWidth',0.5)
hold off
xlabel('EKE (cm^2 s^{-2})')
ylabel('Latitude')
title('Zonally-averaged EKE','FontSize',10)
leg = legend('POP unadjusted','CMEMS obs.','POP composite','location','southeast');

saveas(fig4,['EKE_zonavg_model_satobs_composite_POP_',plot_name,'_',datestr_start,'_',datestr_end_minus1,'_tavg.pdf'])
close(fig4)


% compute adjusted composite temperature/heat fluxes

integrated_heat_flux_mesoscale_all_tmean_array = integrated_heat_flux_mesoscale_tmean_tmean_array + integrated_heat_flux_mesoscale_lowfreq_tmean_array + integrated_heat_flux_mesoscale_highfreq_tmean_array;
local_heat_flux_mesoscale_all_tmean_array = diff(integrated_heat_flux_mesoscale_all_tmean_array,1,2)./(111100*repmat(cosd(lat_transects_vec),[1 length(lon_in_range_all)]).*repmat(diff(lon_in_range_vel_all'),[length(lat_transects_vec) 1]));
local_heat_flux_mesoscale_all_tmean_array(abs(land_mask_local_plot) < 1e-10) = NaN;
local_heat_flux_mesoscale_all_array = diff(cumint_heat_flux_mesoscale_tmean_across_transect_array + cumint_heat_flux_mesoscale_lowfreq_across_transect_array + cumint_heat_flux_mesoscale_highfreq_across_transect_array,1,2)./(111100*repmat(repmat(cosd(lat_transects_vec),[1 length(lon_in_range_all)]).*repmat(diff(lon_in_range_vel_all'),[length(lat_transects_vec) 1]),[1 1 length(time_datenum_in_range)]));
local_heat_flux_mesoscale_all_array(abs(land_mask_local_plot) < 1e-10) = NaN;


half_power_adj = 1;
local_heat_flux_mesoscale_all_tmean_smoothed_array = NaN(size(local_heat_flux_mesoscale_all_tmean_array));
local_heat_flux_mesoscale_all_smoothed_array = NaN(size(local_heat_flux_mesoscale_all_array));
for transect_ind = 1:length(lat_transects_vec)
    spatial_scale_separation = spatial_scale_separation_vec(transect_ind);
    if sum(diff(lon_in_range_all)) >= 359
        curr_size_array = size(local_heat_flux_mesoscale_all_tmean_array);
        buffer_size = round(spatial_scale_separation/(mean(diff(lon_in_range_all))));
        buffered_ind = [(curr_size_array(2) + ((-buffer_size + 1):1:0)'); (1:1:curr_size_array(2))'; (1:1:buffer_size)'];
        unbuffered_ind = buffer_size + (1:1:curr_size_array(2))';
        [curr_array,~,~] = bandpass_err_fcn(local_heat_flux_mesoscale_all_tmean_array(transect_ind,buffered_ind),2,grid_res,1/(4*grid_res*length(lon_in_range_vel_all))/half_power_adj,(1/(2*spatial_scale_separation))*half_power_adj,2,1,1,0,1);
        local_heat_flux_mesoscale_all_tmean_smoothed_array(transect_ind,:) = curr_array(:,unbuffered_ind);
        [curr_array,~,~] = bandpass_err_fcn(local_heat_flux_mesoscale_all_array(transect_ind,buffered_ind,:),2,grid_res,1/(4*grid_res*length(lon_in_range_vel_all))/half_power_adj,(1/(2*spatial_scale_separation))*half_power_adj,2,1,1,0,1);
        local_heat_flux_mesoscale_all_smoothed_array(transect_ind,:,:) = curr_array(:,unbuffered_ind,:);
    else
        [curr_array,~,~] = bandpass_err_fcn(local_heat_flux_mesoscale_all_tmean_array(transect_ind,:),2,grid_res,1/(4*grid_res*length(lon_in_range_vel_all))/half_power_adj,(1/(2*spatial_scale_separation))*half_power_adj,2,1,1,0,1);
        local_heat_flux_mesoscale_all_tmean_smoothed_array(transect_ind,:) = curr_array;
        [curr_array,~,~] = bandpass_err_fcn(local_heat_flux_mesoscale_all_array(transect_ind,:,:),2,grid_res,1/(4*grid_res*length(lon_in_range_vel_all))/half_power_adj,(1/(2*spatial_scale_separation))*half_power_adj,2,1,1,0,1);
        local_heat_flux_mesoscale_all_smoothed_array(transect_ind,:,:) = curr_array;
    end
end

% compute composite temperature/heat fluxes
curr_array_zeronans = local_heat_flux_mesoscale_all_smoothed_array;
curr_array_zeronans((isnan(curr_array_zeronans) == 1) | (abs(curr_array_zeronans) < 1e-10)) = 0;
local_heat_flux_adj_tmean_meso_all_smoothed_array = sum(in_adj_range_mask.*curr_array_zeronans,3)./(sum(in_adj_range_mask,3));


% c_levels = (-1):0.02:1;     % specified levels for contours
cmap = colormap(0.85*bcyr(80));    % colormap for contours
cmap = [((2*cmap(1,:)) + ((-1)*cmap(3,:))); cmap; (((-1)*cmap(78,:)) + (2*cmap(80,:)))];

c_levels = [(-2); ((-0.2):0.005:0.2)'; 2];
[irregular_clevels_array,~] = irregular_clevels_plot((1e-9)*local_heat_flux_adj_tmean_meso_all_smoothed_array,c_levels);

fig5 = figure(5);
close(figure(1))
colormap(cmap)
[~,cbar] = pcolor_plot_with_landmask(fig5,lon_in_range_pcolor_plot',lat_transects_vec_lowbound,[zeros([length(lat_transects_vec_lowbound) 1]) [irregular_clevels_array; zeros([1 (length(lon_in_range_pcolor_plot) - 1)])]],0:1:(length(c_levels) - 1),cmap,'southoutside',lon_bounds_all(1),lon_bounds_all(2),min(lat_transects_vec_lowbound),max(lat_transects_vec) + 0.5);
shading flat
set(gca,'FontSize',12,'xlim',lon_bounds_all,'xtick',(-360):(10*ceil((diff(lon_bounds_all)/10)/10)):360,'ylim',[min(lat_transects_vec_lowbound) (max(lat_transects_vec) + 0.5)],'ydir','normal','DataAspectRatio',[1 mean(cosd(lat_transects_vec)) 1],'clim',[0 size(cmap,1)])
hold on
line(lon_bounds_all,[0 0],[0 0],'Color','k','LineStyle','-','LineWidth',1)
hold off
xlabel('Longitude')
ylabel('Latitude')
title({'POP depth-integrated temp. flux associated with composite mesoscale,'; ['starting from ',num2str(lon_bounds_all(1)),'^{\circ} longitude, averaged ',datestr_title_start,' to ',datestr_title_end_minus1]; ''},'FontSize',8)
cbar = colorbar('southoutside');
set(cbar,'ticks',1:10:81,'FontSize',12,'xticklabel',{'-0.2' '-0.15' '-0.1' '-0.05' '0' '0.05' '0.1' '0.15' '0.2'})
set(get(cbar,'xlabel'),'String','Temperature flux, TW km^{-1}','FontSize',14)

print(fig5,['Heat_flux_mesoscale_all_composite_satobs_POP_',plot_name,'_local_smoothed_from_',num2str(lon_bounds_all(1)),'_to_',num2str(lon_bounds_all(2)),'_lon_',datestr_start,'_',datestr_end_minus1,'_tavg.bmp'],'-dbmp','-r300')
close(fig5)



% integrate temperature/heat flux zonally

delta_lon_grid = 0.1;

curr_array = local_heat_flux_mesoscale_all_tmean_smoothed_array;
curr_mask_all = ones(size(curr_array));
curr_mask_all((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
curr_array_zeronans = curr_array;
curr_array_zeronans((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
local_heat_flux_mesoscale_transect_avg = sum(repmat(111100*cosd(lat_transects_vec)*delta_lon_grid,[1 size(curr_array_zeronans,2)]).*curr_mask_all.*curr_array_zeronans,2);
curr_array = local_heat_flux_adj_tmean_meso_all_smoothed_array;
curr_mask_all = ones(size(curr_array));
curr_mask_all((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
curr_array_zeronans = curr_array;
curr_array_zeronans((isnan(curr_array) == 1) | (abs(curr_array) < 1e-10)) = 0;
local_heat_flux_mesoscale_adj_transect_avg = sum(repmat(111100*cosd(lat_transects_vec)*delta_lon_grid,[1 size(curr_array_zeronans,2)]).*curr_mask_all.*curr_array_zeronans,2);

local_heat_flux_mesoscale_adj_transect_avg(abs(sum(curr_mask_all,2)) < 1e-5) = NaN;


fig6 = figure(6);
h = plot((1e-15)*local_heat_flux_mesoscale_transect_avg,lat_transects_vec,'r-',(1e-15)*local_heat_flux_mesoscale_adj_transect_avg,lat_transects_vec,'r--');
set(h,'LineWidth',1)
set(h(:),'Color',[0.8 0 0])
set(gca,'xlim',[-1 1],'xtick',(-1):0.5:1,'ylim',[min(lat_transects_vec) max(lat_transects_vec)],'DataAspectRatio',[3.5 (max(lat_transects_vec) - min(lat_transects_vec)) 1],'FontSize',12)
hold on
line([-1 1],[0 0],'Color',[0 0 0],'LineWidth',0.5)
line([0 0],[min(lat_transects_vec) max(lat_transects_vec)],'Color',[0 0 0],'LineWidth',0.5)
hold off
xlabel('Heat transport (PW)')
ylabel('Latitude')
title({'Mesoscale contribution to HT composite estimate'; ' '},'FontSize',10)
leg = legend('POP unadjusted','POP composite','location','northeast');

saveas(fig6,['Heat_transport_mesoscale_model_composite_POP_',plot_name,'_',datestr_start,'_',datestr_end_minus1,'_tavg.pdf'])
close(fig6)