EKE_variability_map_global_POP.m

% map time variability EKE (in specified frequency range) in POP simulation

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



% SSH file identifier
SSH_file_id = 'POP_5day_global';

% boundaries of region to be plotted

longwest = -340;
longeast = 20;
latsouth = -90;
latnorth = 90;


% frequency range to filter EKE variability
low_freq_bound = 1/(4*32*365);
high_freq_bound = 1/426;

% map projection
projection = 'robinson';


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


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


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_lon_range_ind = (1:1:size(longitude,1))';
in_lat_range_ind = (1:1:size(longitude,2))';
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');


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

east_end_ind = find(mod(lon_in_range(:,320) - mean([longwest longeast]) + 180,360) - 180 >= 180 - grid_res);
lon_in_range = mod(lon_in_range - mean([longwest longeast]) + 180,360) + mean([longwest longeast]) - 180;
lon_in_range = lon_in_range([((east_end_ind + 1):1:size(lon_in_range,1))'; (1:1:east_end_ind)'],:);
lat_in_range = lat_in_range([((east_end_ind + 1):1:size(lon_in_range,1))'; (1:1:east_end_ind)'],:);

% compute grid differences in global domain
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)) (lon_in_range_midpts(:,size_array(2)) + ((mod(diff([lon_in_range_midpts(:,size_array(2)) lon_in_range_midpts(mod(size_array(1) - (2*east_end_ind) - (1:1:size_array(1))' - 1,size_array(1)) + 1,size_array(2))],1,2) + 180,360) - 180)/2))];
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)) (lat_in_range_midpts(:,size_array(2)) + (diff([lat_in_range_midpts(:,size_array(2)) lat_in_range_midpts(mod(size_array(1) - (2*east_end_ind) - (1:1:size_array(1))' - 1,size_array(1)) + 1,size_array(2))],1,2)/2))];
dx = 111100*(abs((cosd(lat_in_range_midpts).*(mod(diff([(lon_in_range_aug_x(:,2:size_array(2)) - (diff(lon_in_range_aug_x,1,2)/2)) (lon_in_range_aug_x(:,size_array(2)) + ((mod(diff([lon_in_range_aug_x(:,size_array(2)) lon_in_range_aug_x(mod(size_array(1) + 1 - (2*east_end_ind) - (1:1:(size_array(1) + 1))' - 1,size_array(1)) + 1,size_array(2))],1,2) + 180,360) - 180)/2))],1,1) + 180,360) - 180)) + (1i*diff([(lat_in_range_aug_x(:,2:size_array(2)) - (diff(lat_in_range_aug_x,1,2)/2)) (lat_in_range_aug_x(:,size_array(2)) + (diff([lat_in_range_aug_x(:,size_array(2)) lat_in_range_aug_x(mod(size_array(1) + 1 - (2*east_end_ind) - (1:1:(size_array(1) + 1))' - 1,size_array(1)) + 1,size_array(2))],1,2)/2))],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])])));
dy = 111100*(abs((cosd(lat_in_range_midpts).*(mod(diff([(lon_in_range_aug_x(1:size_array(1),:) + (diff(lon_in_range_aug_x,1,1)/2)) (lon_in_range_aug_x(mod(size_array(1) + 1 - (2*east_end_ind) - (1:1:size_array(1))' - 1,size_array(1)) + 1,size_array(2)) + ((mod(diff(lon_in_range_aug_x(mod(size_array(1) + 1 - (2*east_end_ind) - (1:1:(size_array(1) + 1))' - 1,size_array(1)) + 1,size_array(2)),1,1) + 180,360) - 180)/2))],1,2) + 180,360) - 180)) + (1i*diff([(lat_in_range_aug_x(1:size_array(1),:) + (diff(lat_in_range_aug_x,1,1)/2)) (lat_in_range_aug_x(mod(size_array(1) + 1 - (2*east_end_ind) - (1:1:size_array(1))' - 1,size_array(1)) + 1,size_array(2)) + (diff(lat_in_range_aug_x(mod(size_array(1) + 1 - (2*east_end_ind) - (1:1:(size_array(1) + 1))' - 1,size_array(1)) + 1,size_array(2)),1,1)/2))],1,2))));


% fill in lon, lat over land areas that have NaN values
lon_in_range(1,:) = interp1(find(isnan(lon_in_range(1,:)) == 0),lon_in_range(1,isnan(lon_in_range(1,:)) == 0),1:1:size(lon_in_range,2));
lat_in_range(1,:) = interp1(find(isnan(lat_in_range(1,:)) == 0),lat_in_range(1,isnan(lat_in_range(1,:)) == 0),1:1:size(lat_in_range,2));
lon_in_range(size(lon_in_range,1),:) = interp1(find(isnan(lon_in_range(size(lon_in_range,1),:)) == 0),lon_in_range(size(lon_in_range,1),isnan(lon_in_range(size(lon_in_range,1),:)) == 0),1:1:size(lon_in_range,2));
lat_in_range(size(lat_in_range,1),:) = interp1(find(isnan(lat_in_range(size(lat_in_range,1),:)) == 0),lat_in_range(size(lat_in_range,1),isnan(lat_in_range(size(lat_in_range,1),:)) == 0),1:1:size(lat_in_range,2));
for j = 1:size(lon_in_range,2)
    curr_lon = lon_in_range(:,j);
    curr_lat = lat_in_range(:,j);
    good_coord_ind = find(isnan(curr_lon) == 0);
    lon_in_range(:,j) = interp1(good_coord_ind,curr_lon(good_coord_ind),(1:1:size(lon_in_range,1))');
    lat_in_range(:,j) = interp1(good_coord_ind,curr_lat(good_coord_ind),(1:1:size(lat_in_range,1))');
    if length(find(isnan(lon_in_range(:,j)) == 0)) < size(lon_in_range,1)
        first_good_pair_ind = find((isnan(lon_in_range(1:(size(lon_in_range,1) - 1),j)) == 0) & (isnan(lon_in_range(2:size(lon_in_range,1),j)) == 0),1,'first');
        curr_lon = lon_in_range(mod(first_good_pair_ind + (1:1:size(lon_in_range,1))' - 1,size(lon_in_range,1)) + 1,j);
        curr_lat = lat_in_range(mod(first_good_pair_ind + (1:1:size(lat_in_range,1))' - 1,size(lat_in_range,1)) + 1,j);
        good_coord_ind = find(isnan(curr_lon) == 0);
        lon_in_range(mod(first_good_pair_ind + (1:1:size(lon_in_range,1))' - 1,size(lon_in_range,1)) + 1,j) = mod(interp1(good_coord_ind,mod(curr_lon(good_coord_ind) - (curr_lon(1) - 1e-5),360) + (curr_lon(1) - 1e-5),(1:1:size(lon_in_range,1))') - longwest,360) + longwest;
        lat_in_range(mod(first_good_pair_ind + (1:1:size(lat_in_range,1))' - 1,size(lat_in_range,1)) + 1,j) = interp1(good_coord_ind,curr_lat(good_coord_ind),(1:1:size(lat_in_range,1))');
    end
end


in_time_range_ind_all = in_time_range_ind;

% n_good_in_subset_SSH = NaN([length(in_lon_range_ind) length(in_lat_range_ind) n_t_subsets]);
% n_good_in_subset_EKE = NaN([length(in_lon_range_ind) length(in_lat_range_ind) n_t_subsets]);
% SSH_mean_in_range_subset = NaN([length(in_lon_range_ind) length(in_lat_range_ind) n_t_subsets]);
% EKE_mean_in_range_subset = NaN([length(in_lon_range_ind) length(in_lat_range_ind) n_t_subsets]);
    
%     curr_file_nums = ((((curr_t_subset - 1)*n_years_in_time_subset) + 1):1:(curr_t_subset*n_years_in_time_subset))';
%     curr_t_in_subset_ind = find((file_num_in_range_vec_all >= min(curr_file_nums)) & (file_num_in_range_vec_all <= max(curr_file_nums)));
curr_t_in_subset_ind = (1:1:length(file_num_in_range_vec_all))';
in_time_range_ind = in_time_range_ind_all(curr_t_in_subset_ind);
time_datenum_in_range = time_datenum_in_range_all(curr_t_in_subset_ind);
%     time_bound_datenum_in_range = time_bound_datenum_in_range_all(curr_t_in_subset_ind,:);
file_num_in_range_vec = file_num_in_range_vec_all(curr_t_in_subset_ind);
in_file_ind_in_range_vec = in_file_ind_in_range_vec_all(curr_t_in_subset_ind);

load(['SSH_mean_',SSH_file_id,'_1978_2009.mat'],'SSH_mean_in_range')


n_i_subsets = 10;
n_j_subsets = 10;

i_size_subset = ceil(size(lon_in_range,1)/n_i_subsets);
j_size_subset = ceil(size(lon_in_range,2)/n_j_subsets);
n_i_subsets = ceil(size(lon_in_range,1)/i_size_subset);
n_j_subsets = ceil(size(lon_in_range,2)/j_size_subset);
base_i_ind_subset = mod(east_end_ind + ((i_size_subset*(0:1:(n_i_subsets - 1))') + 1) - 1,size(lon_in_range,1)) + 1;
base_j_ind_subset = (j_size_subset*(0:1:(n_j_subsets - 1))') + 1;
buffer_ind_size = 5;
start_i_ind_subset = base_i_ind_subset - buffer_ind_size;
end_i_ind_subset = base_i_ind_subset + (i_size_subset - 1) + buffer_ind_size;
start_j_ind_subset = base_j_ind_subset - buffer_ind_size;
end_j_ind_subset = base_j_ind_subset + (j_size_subset - 1) + buffer_ind_size;

EKE_stddev_in_freq_range = NaN(size(lon_in_range));
for i_subset = 1:n_i_subsets
    for j_subset = 1:n_j_subsets
        
        in_lon_range_ind = mod((start_i_ind_subset(i_subset):1:end_i_ind_subset(i_subset))' - 1,size(lon_in_range,1)) + 1;
        if j_subset == 1
            in_lat_range_ind = [(buffer_ind_size:(-1):1)'; (1:1:end_j_ind_subset(j_subset))'];
        elseif j_subset == n_j_subsets
            in_lat_range_ind = [(start_j_ind_subset(j_subset):1:size(lon_in_range,2))'; (size(lon_in_range,2) + (0:(-1):(-buffer_ind_size + 1)))'];
        else
            in_lat_range_ind = (start_j_ind_subset(j_subset):1:end_j_ind_subset(j_subset))';
        end
        
        
        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: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];
    
            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 in_lon_range_ind(length(in_lon_range_ind)) - in_lon_range_ind(1) < 0
                lon_range_ind_start = in_lon_range_ind(1);
                lon_range_ind_span = max(in_lon_range_ind) - in_lon_range_ind(1) + 1;
                lon_range_ind_start_2 = 1;
                lon_range_ind_span_2 = in_lon_range_ind(length(in_lon_range_ind)) - 1 + 1;
            else
                lon_range_ind_start = min(in_lon_range_ind);
                lon_range_ind_span = max(in_lon_range_ind) - min(in_lon_range_ind) + 1;
                lon_range_ind_start_2 = 1;
                lon_range_ind_span_2 = 0;
            end
            if isempty(find(diff(in_lat_range_ind) < 0,1,'first')) == 0
                if j_subset == 1
                    lat_range_ind_start = 1;
                    lat_range_ind_span = in_lat_range_ind(1) - 1 + 1;
                    lat_range_ind_start_2 = 1;
                    lat_range_ind_span_2 = in_lat_range_ind(length(in_lat_range_ind)) - 1 + 1;
                elseif j_subset == n_j_subsets
                    lat_range_ind_start = in_lat_range_ind(1);
                    lat_range_ind_span = max(in_lat_range_ind) - in_lat_range_ind(1) + 1;
                    lat_range_ind_span_2 = max(in_lat_range_ind) - in_lat_range_ind(length(in_lat_range_ind)) + 1;
                    lat_range_ind_start_2 = in_lat_range_ind(length(in_lat_range_ind));
                end
            else
                lat_range_ind_start = min(in_lat_range_ind);
                lat_range_ind_span = max(in_lat_range_ind) - min(in_lat_range_ind) + 1;
                lat_range_ind_start_2 = 1;
                lat_range_ind_span_2 = 0;
            end
%             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 = [lon_range_ind_start lat_range_ind_start 1 min(curr_in_time_range_ind)];
            count_vec = [lon_range_ind_span lat_range_ind_span 1 time_range_ind_span];
            if isempty(NaN(count_vec)) == 0
                curr_SSH_in_range = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
            else
                curr_SSH_in_range = NaN(count_vec([1 2 4]));
            end
            if ((lon_range_ind_span_2 > 1e-5) || (lat_range_ind_span_2 > 1e-5))
                start_vec = [lon_range_ind_start_2 lat_range_ind_start 1 min(curr_in_time_range_ind)];
                count_vec = [lon_range_ind_span_2 lat_range_ind_span 1 time_range_ind_span];
                if isempty(NaN(count_vec)) == 0
                    curr_SSH_in_range_2 = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
                else
                    curr_SSH_in_range_2 = NaN(count_vec([1 2 4]));
                end
                if j_subset == 1
                    curr_SSH_in_range = flip(curr_SSH_in_range,2);
                    curr_SSH_in_range_2 = flip(curr_SSH_in_range_2,2);
                    start_vec = [lon_range_ind_start lat_range_ind_start_2 1 min(curr_in_time_range_ind)];
                    count_vec = [lon_range_ind_span lat_range_ind_span_2 1 time_range_ind_span];
                    if isempty(NaN(count_vec)) == 0
                        curr_SSH_in_range_3 = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
                    else
                        curr_SSH_in_range_3 = NaN(count_vec([1 2 4]));
                    end
                    start_vec = [lon_range_ind_start_2 lat_range_ind_start_2 1 min(curr_in_time_range_ind)];
                    count_vec = [lon_range_ind_span_2 lat_range_ind_span_2 1 time_range_ind_span];
                    if isempty(NaN(count_vec)) == 0
                        curr_SSH_in_range_4 = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
                    else
                        curr_SSH_in_range_4 = NaN(count_vec([1 2 4]));
                    end                        
                elseif j_subset == n_j_subsets
                    start_vec = [(size(lon_in_range,1) - (lon_range_ind_start - 1 + lon_range_ind_span) + 1) lat_range_ind_start_2 1 min(curr_in_time_range_ind)];
                    count_vec = [lon_range_ind_span lat_range_ind_span_2 1 time_range_ind_span];
                    if isempty(NaN(count_vec)) == 0
                        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);
                    else
                        curr_SSH_in_range_3 = NaN(count_vec([1 2 4]));
                    end
                    start_vec = [(size(lon_in_range,1) - (lon_range_ind_start_2 - 1 + lon_range_ind_span_2) + 1) lat_range_ind_start_2 1 min(curr_in_time_range_ind)];
                    count_vec = [lon_range_ind_span_2 lat_range_ind_span_2 1 time_range_ind_span];
                    if isempty(NaN(count_vec)) == 0
                        curr_SSH_in_range_4 = 0.01*flip(flip(squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4])))),1),2);
                    else
                        curr_SSH_in_range_4 = NaN(count_vec([1 2 4]));
                    end
                else
                    start_vec = [lon_range_ind_start lat_range_ind_start_2 1 min(curr_in_time_range_ind)];
                    count_vec = [lon_range_ind_span lat_range_ind_span_2 1 time_range_ind_span];
                    if isempty(NaN(count_vec)) == 0
                        curr_SSH_in_range_3 = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
                    else
                        curr_SSH_in_range_3 = NaN(count_vec([1 2 4]));
                    end
                    start_vec = [lon_range_ind_start_2 lat_range_ind_start_2 1 min(curr_in_time_range_ind)];
                    count_vec = [lon_range_ind_span_2 lat_range_ind_span_2 1 time_range_ind_span];
                    if isempty(NaN(count_vec)) == 0
                        curr_SSH_in_range_4 = 0.01*squeeze(double(ncread(curr_filename,'SSH',start_vec([1 2 4]),count_vec([1 2 4]))));
                    else
                        curr_SSH_in_range_4 = NaN(count_vec([1 2 4]));
                    end
                end
                
                curr_SSH_in_range = [curr_SSH_in_range curr_SSH_in_range_3; curr_SSH_in_range_2 curr_SSH_in_range_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
        
        SSH_in_range(abs(SSH_in_range - (-0.01)) < 1e-10) = NaN;
        
        
        size_array = size(SSH_in_range);
        
        
        % spatial filter to remove wavenumbers higher than those possible in gridded altimetry
    
        delta_x = grid_res;
        delta_y = grid_res;
        altimetry_grid_res = 0.25;
    
    
        nan_mask = ones(size(SSH_in_range));
        nan_mask(isnan(SSH_in_range) == 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;    
    
    
    %     SSH_in_range_lp_expanded = NaN([1 2 1].*(size(SSH_in_range)));
        buffer_size = 100;
        SSH_in_range_lp_expanded = NaN(size(SSH_in_range) + [buffer_size buffer_size 0]);
        n_t_subset = 10;
        length_t_subset = ceil(size_array(3)/n_t_subset);
        subset_start_ind_t = (1 + (ceil(size_array(3)/n_t_subset)*(0:1:(n_t_subset - 1))))';
        subset_end_ind_t = min([(length_t_subset + (ceil(size_array(3)/n_t_subset)*(0:1:(n_t_subset - 1))))' ((size_array(3))*ones(n_t_subset,1))],[],2);
        for t_subset = 1:n_t_subset
            curr_t_min = subset_start_ind_t(t_subset);
            curr_t_max = subset_end_ind_t(t_subset);
            
%             curr_lon_expanded = [lon_in_range NaN([size_array(1) buffer_size])];
%             curr_lon_expanded = [NaN([ceil(buffer_size/2) (size_array(2) + buffer_size)]); curr_lon_expanded; NaN([ceil(buffer_size/2) (size_array(2) + buffer_size)])];
%             curr_lat_expanded = [lat_in_range NaN([size_array(1) buffer_size])];
%             curr_lat_expanded = [NaN([ceil(buffer_size/2) (size_array(2) + buffer_size)]); curr_lat_expanded; NaN([ceil(buffer_size/2) (size_array(2) + buffer_size)])];
            curr_SSH_expanded = [SSH_in_range(:,:,curr_t_min:curr_t_max) NaN([size_array(1) buffer_size (curr_t_max - curr_t_min + 1)])];
            curr_SSH_expanded = [NaN([ceil(buffer_size/2) (size_array(2) + buffer_size) (curr_t_max - curr_t_min + 1)]); curr_SSH_expanded; NaN([ceil(buffer_size/2) (size_array(2) + buffer_size) (curr_t_max - curr_t_min + 1)])];
            curr_SSH_in_range_lp_expanded = Gaussian_2D_filter(curr_SSH_expanded,mean([abs(delta_x) abs(delta_y)]),altimetry_grid_res,2*altimetry_grid_res,altimetry_grid_res,0,1);
            SSH_in_range_lp_expanded(:,:,curr_t_min:curr_t_max) = curr_SSH_in_range_lp_expanded;
            
            disp(['t_subset = ',num2str(t_subset),' (out of ',num2str(n_t_subset),')'])
        end
        clear curr_SSH*x_bp curr_SSH*lp_expanded
    
        curr_nan_2d_expanded = [nan_mask_array_2d ((-1)*ones([size_array(1) buffer_size]))];
        curr_nan_2d_expanded = [((-1)*ones([ceil(buffer_size/2) (size_array(2) + buffer_size)])); curr_nan_2d_expanded; ((-1)*ones([ceil(buffer_size/2) (size_array(2) + buffer_size)]))];
        nan_2d_spat_bp_expanded = Gaussian_2D_filter(curr_nan_2d_expanded,mean([abs(delta_x) abs(delta_y)]),altimetry_grid_res,2*altimetry_grid_res,altimetry_grid_res,0,1);
        SSH_in_range_lp = SSH_in_range_lp_expanded(ceil(buffer_size/2) + (1:1:size_array(1))',1:1:size_array(2),:);
        nan_2d_spat_bp = nan_2d_spat_bp_expanded(ceil(buffer_size/2) + (1:1:size_array(1))',1:1:size_array(2));
        
        clear *_expanded
    
        SSH_in_range_lp(repmat(abs(nan_2d_spat_bp - 1e-10),[1 1 size_array(3)]) > 0.1) = NaN;
        
        SSH_in_range = SSH_in_range_lp;
        clear SSH_in_range_lp
    
        
%         load(['SSH_mean_',SSH_file_id,'_1978_2009.mat'],'SSH_mean_in_range')
        if j_subset == n_j_subsets
            curr_SSH_mean_in_range = [SSH_mean_in_range(mod(in_lon_range_ind - east_end_ind - 1,size(lon_in_range,1)) + 1,in_lat_range_ind(1:(length(in_lat_range_ind) - buffer_ind_size))) SSH_mean_in_range(mod((size(lon_in_range,1) + 1 - in_lon_range_ind) - east_end_ind - 1,size(lon_in_range,1)) + 1,size(lon_in_range,2) + (0:(-1):(-buffer_ind_size + 1))')];
        else
            curr_SSH_mean_in_range = SSH_mean_in_range(mod(in_lon_range_ind - east_end_ind - 1,size(lon_in_range,1)) + 1,in_lat_range_ind);
        end
%         clear SSH_mean_in_range
        curr_lat_in_range_midpts = lat_in_range_midpts(mod(in_lon_range_ind - east_end_ind - 1,size(lon_in_range,1)) + 1,in_lat_range_ind);
        curr_dx = dx(mod(in_lon_range_ind - east_end_ind - 1,size(lon_in_range,1)) + 1,in_lat_range_ind);
        curr_dy = dy(mod(in_lon_range_ind - east_end_ind - 1,size(lon_in_range,1)) + 1,in_lat_range_ind);
        SSH_nomean_in_range = SSH_in_range - repmat(curr_SSH_mean_in_range,[1 1 size(SSH_in_range,3)]);
        
        
        size_array = size(SSH_nomean_in_range);
        
        curr_SSH_nomean_expanded = [SSH_nomean_in_range NaN([size_array(1) buffer_size size_array(3)])];
        curr_SSH_nomean_expanded = [NaN([ceil(buffer_size/2) (size_array(2) + buffer_size) size_array(3)]); curr_SSH_nomean_expanded; NaN([ceil(buffer_size/2) (size_array(2) + buffer_size) size_array(3)])];
        
        diff_x_SSH_nomean_in_r = diff(curr_SSH_nomean_expanded,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);
        diff_y_SSH_nomean_in_r = diff(curr_SSH_nomean_expanded,1,2);
        diff_y_SSH_nomean_in_range = 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);
        diff_x_SSH_nomean_in_range = diff_x_SSH_nomean_in_range(ceil(buffer_size/2) + (1:1:size_array(1))',1:1:size_array(2),:);
        diff_y_SSH_nomean_in_range = diff_y_SSH_nomean_in_range(ceil(buffer_size/2) + (1:1:size_array(1))',1:1:size_array(2),:);
        
        
        f = 2*((2*pi)/86164)*sind(curr_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(curr_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(curr_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(curr_lat_in_range_midpts));
        lat_mask(abs(curr_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;
        
        nan_mask_EKE = ones(size(EKE_in_range));
        nan_mask_EKE(isnan(EKE_in_range) == 1) = 0;
        EKE_in_range_nans_zeroed = EKE_in_range;
        EKE_in_range_nans_zeroed(isnan(EKE_in_range) == 1) = 0;
        
        
        % remove seasonal averages
        n_bins = round(365/delta_t);
        EKE_in_range_noseason = NaN(size(EKE_in_range));
        for curr_bin = 1:n_bins
            in_bin_ind = find(abs(mod(time_datenum_in_range - ((curr_bin - 0.5)*(365/n_bins)) + 180,360) - 180) <= (delta_t/2));
            curr_nan_mask = nan_mask_EKE(:,:,in_bin_ind);
            curr_EKE_nans_zeroed = EKE_in_range_nans_zeroed(:,:,in_bin_ind);
            EKE_in_range_noseason(:,:,in_bin_ind) = EKE_in_range(:,:,in_bin_ind) - repmat(sum(curr_nan_mask.*curr_EKE_nans_zeroed,3)./(sum(curr_nan_mask,3)),[1 1 length(in_bin_ind)]);
        end
        % add the EKE mean back
        EKE_mean_in_range = sum(nan_mask_EKE.*EKE_in_range_nans_zeroed,3)./(sum(nan_mask_EKE,3));
        EKE_in_range_noseason = EKE_in_range_noseason + repmat(EKE_mean_in_range,[1 1 size(EKE_in_range_noseason,3)]);
        
        % filter for specified frequency range
        steepness_factor = 5;
        half_power_adj = exp(erfinv((2^(1/2)) - 1)/steepness_factor);
        [EKE_in_range_filtered,~,~] = bandpass_err_fcn(EKE_in_range_noseason,3,delta_t,low_freq_bound/half_power_adj,high_freq_bound*half_power_adj,steepness_factor,1,1,1,0);
        
        curr_nan_mask = ones(size(EKE_in_range_filtered));
        curr_nan_mask((isnan(EKE_in_range_filtered) == 1) | (abs(EKE_in_range_filtered) < 1e-5)) = 0;
        EKE_in_range_filtered(abs(curr_nan_mask) < 1e-5) = 0;
        curr_array_mean = sum(curr_nan_mask.*EKE_in_range_filtered,3)./(sum(curr_nan_mask,3));
        curr_array_stddev = (sum(curr_nan_mask.*((EKE_in_range_filtered - repmat(curr_array_mean,[1 1 size(EKE_in_range_filtered,3)])).^2),3)./(sum(curr_nan_mask,3))).^(1/2);
        
        
        curr_i_ind = mod(base_i_ind_subset(i_subset) - east_end_ind + (0:1:(i_size_subset - 1))' - 1,size(lon_in_range,1)) + 1;
        curr_j_ind = base_j_ind_subset(j_subset) + (0:1:(j_size_subset - 1))';
        curr_i_ind = curr_i_ind(curr_i_ind <= size(lon_in_range,1));
        curr_j_ind = curr_j_ind(curr_j_ind <= size(lon_in_range,2));
        
        EKE_stddev_in_freq_range(curr_i_ind,curr_j_ind) = curr_array_stddev(buffer_ind_size + (1:1:length(curr_i_ind))',buffer_ind_size + (1:1:length(curr_j_ind))');

    end
    
    disp(['Completed subset ',num2str(i_subset),' (out of ',num2str(n_i_subsets),')'])
    
end


time_datenum_in_range = time_datenum_in_range_all;
time_bound_datenum_in_range = time_bound_datenum_in_range_all;

save(['EKE_stddev_',num2str(1/high_freq_bound),'_',num2str(1/low_freq_bound),'_periods_',SSH_file_id,'_',num2str(time_range_start(1)),'_',num2str(time_range_end(1) - 1),'.mat'],'EKE_stddev_in_freq_range')
% load(['EKE_stddev_',num2str(1/high_freq_bound),'_',num2str(1/low_freq_bound),'_periods_',SSH_file_id,'_',num2str(time_range_start(1)),'_',num2str(time_range_end(1) - 1),'.mat'])


% re-grid to lat/lon grid

delta_grid = 0.1;
lon_grid_vec = lon_in_range_midpts(:,450);
% lat_grid_vec = ((delta_grid*round(min(min(lat_in_range_midpts))/delta_grid)) + (delta_grid/2)):delta_grid:max(max(lat_in_range_midpts));
lat_grid_vec = (delta_grid*round(min(min(lat_in_range_midpts))/delta_grid)):delta_grid:max(max(lat_in_range_midpts));

[~,closest_ind] = min(abs(mod(lon_grid_vec - (-163.83) + 180,360) - 180));
EKE_stddev_in_freq_range_regridded = NaN([length(lon_grid_vec) length(lat_grid_vec)]);
if latsouth < 28.5
    lat_old_merc = lat_in_range_midpts(closest_ind,:);
    [~,max_merc_ind] = min(abs(lat_old_merc - 28.52));
    lat_old_merc = lat_old_merc(1:max_merc_ind);
    lat_new_merc_ind = find(lat_grid_vec <= lat_old_merc(max_merc_ind));
    EKE_stddev_in_freq_range_regridded(:,lat_new_merc_ind) = (interp1(lat_old_merc',EKE_stddev_in_freq_range(:,1:max_merc_ind)',lat_grid_vec(lat_new_merc_ind)'))';
    j_nonmerc_start_ind = max(lat_new_merc_ind);
else
    j_nonmerc_start_ind = 1;
end

% if j_start_ind <= 1480
%     n_j_no_regridding = 1480 + 1 - j_start_ind;
%     EKE_stddev_in_range_regridded(:,1:min([n_j_no_regridding size_array(2)])) = EKE_stddev_in_range(:,1:min([n_j_no_regridding size_array(2)]));
% end
% if j_start_ind - 1 + size_array(2) > 1480
if latnorth > 28.5
    for i_grid = 1:length(lon_grid_vec)
%         for j_grid = (1481 + 1 - j_start_ind):length(lat_grid_vec)
        j_grid = j_nonmerc_start_ind;
        while j_grid <= length(lat_grid_vec)
            curr_lon_center = lon_grid_vec(i_grid);
            curr_lat_center = lat_grid_vec(j_grid);
            
            j_region_spacing = 200;
            % limit region being searched (to save time)
            if ((j_grid == j_nonmerc_start_ind) || (mod(j_grid,j_region_spacing) == 1))
                region_start_i_ind = max([1 (i_grid - 5)]);
                region_end_i_ind = min([length(lon_grid_vec) (i_grid + 5)]);
                region_start_j_ind = max([1 (j_grid - 5)]);
                region_end_j_ind = min([length(lat_grid_vec) (j_grid - 1 + j_region_spacing + 5)]);
                lon_min_curr_region = min(mod(lon_grid_vec(region_start_i_ind:region_end_i_ind) - (lon_grid_vec(region_start_i_ind) - 180),360) + (lon_grid_vec(region_start_i_ind) - 180)) - delta_grid;
                lon_max_curr_region = max(mod(lon_grid_vec(region_start_i_ind:region_end_i_ind) - (lon_grid_vec(region_start_i_ind) - 180),360) + (lon_grid_vec(region_start_i_ind) - 180)) + delta_grid;
                lat_min_curr_region = min(lat_grid_vec(region_start_j_ind:region_end_j_ind)) - delta_grid;
                lat_max_curr_region = max(lat_grid_vec(region_start_j_ind:region_end_j_ind)) + delta_grid;
                in_region_old_range_ind = find((mod(lon_in_range_midpts - lon_min_curr_region + 180,360) - 180 >= 0) & (mod(lon_in_range_midpts - lon_max_curr_region + 180,360) - 180 <= 0) & (lat_in_range_midpts - lat_min_curr_region >= 0) & (lat_in_range_midpts - lat_max_curr_region <= 0));
                if isempty(in_region_old_range_ind) == 1
                    j_grid = (j_region_spacing*ceil(j_grid/j_region_spacing)) + 1;
                    continue
                end
                in_region_old_i_ind = unique(mod(in_region_old_range_ind - 1,size(lon_in_range_midpts,1)) + 1);
                in_region_old_i_ind = mod(sort(mod(in_region_old_i_ind - (in_region_old_i_ind(1) - (size(lon_in_range_midpts,1)/2)),size(lon_in_range_midpts,1)) + (in_region_old_i_ind(1) - (size(lon_in_range_midpts,1)/2)),'ascend') - 1,size(lon_in_range_midpts,1)) + 1;
                in_region_old_j_ind = unique(ceil(in_region_old_range_ind/size(lon_in_range_midpts,1)));
                lon_in_range_midpts_curr_region = lon_in_range_midpts(in_region_old_i_ind,in_region_old_j_ind);
                lat_in_range_midpts_curr_region = lat_in_range_midpts(in_region_old_i_ind,in_region_old_j_ind);
                EKE_stddev_in_freq_range_curr_region = EKE_stddev_in_freq_range(in_region_old_i_ind,in_region_old_j_ind);
            end

            in_range_ind = find((mod(lon_in_range_midpts_curr_region - curr_lon_center + 180,360) - 180 >= -(0.6*delta_grid/cosd(curr_lat_center))) & (mod(lon_in_range_midpts_curr_region - curr_lon_center + 180,360) - 180 < (0.6*delta_grid/cosd(curr_lat_center))) & (lat_in_range_midpts_curr_region - curr_lat_center >= -0.6*delta_grid) & (lat_in_range_midpts_curr_region - curr_lat_center < 0.6*delta_grid));
            in_range_ind = in_range_ind(isnan(EKE_stddev_in_freq_range_curr_region(in_range_ind)) == 0);

            EKE_stddev_in_freq_range_regridded(i_grid,j_grid) = mean(EKE_stddev_in_freq_range_curr_region(in_range_ind));
            j_grid = j_grid + 1;
        end
        if mod(i_grid,200) == 0
            disp(['i_grid = ',num2str(i_grid),' (out of ',num2str(length(lon_grid_vec)),')'])
        end
    end
end



% c_levels = [(0:10:40) (60:20:100) (100:40:300) (350:50:500) 1000];     % specified levels for contours
c_levels = [(0:12.5:50) (75:25:150) (200:50:350) (450:100:750) 1000];     % specified levels for contours
cmap = colormap(0.85*bcyr(31));
cmap = cmap([(1:4:17)'; (20:1:31)'],:);

[irregular_clevels_array,~] = irregular_clevels_plot((1e4)*real(EKE_stddev_in_freq_range_regridded),c_levels);

fig3 = figure(3);
fig_pos = get(fig3,'Position');
fig_pos(3:4) = 4*fig_pos(3:4);
% fig_paper_pos = get(fig3,'PaperPosition');
% fig_paper_pos(3:4) = [14 10];
set(fig3,'Position',fig_pos,'PaperSize',[11 8.5])
close(figure(1))
[~,axm,cbar] = contourm_plot_with_landmask(fig3,lon_grid_vec',lat_grid_vec',irregular_clevels_array',0:1:(length(c_levels) - 1),cmap,'southoutside',longwest,longeast,latsouth,latnorth,'robinson',[0 mean([longwest longeast]) 0]);
setm(axm,'frame','on','grid','on','flonlimit',[longwest longeast] - mean([longwest longeast]),'flatlimit',[latsouth latnorth],'mlinelocation',30*round(((longeast - longwest)/6)/30),'plinelocation',30,'glinestyle','-','meridianlabel','on','parallellabel','on','mlabellocation',60,'mlabelparallel',latnorth - 0.01,'plabellocation',30,'plabelmeridian',longwest - 10,'FontSize',12)
set(axm,'XColor','none','YColor','none')
children_handles = get(axm,'Children');
for curr_handle = 8:13
    try
        curr_pos = get(children_handles(curr_handle),'Position');
        set(children_handles(curr_handle),'Position',curr_pos + [0 0.07 0])
    catch
    end
end
hold on
linem([0 0],[longwest longeast],'Color','k','LineStyle','-','LineWidth',0.5)
hold off
title({['POP EKE std. dev. ',num2str(1/high_freq_bound),'-',num2str(1/low_freq_bound),' day periods, computed from data ',num2str(time_range_start(1)),'-',num2str(time_range_end(1) - 1),',']; ['and spatially low-pass filtered with a ',num2str(2*altimetry_grid_res),'^ {\circ} threshold']; ' '},'FontSize',12)
set(gca,'clim',[0 size(cmap,1)])
set(cbar,'ticks',[0 4 8 12 16],'xticklabel',{'0' '50' '150' '350' '750'})
set(get(cbar,'label'),'String','Eddy kinetic energy std. dev. (cm^2 s^{-2})','FontSize',14)

print(fig3,['EKE_std_dev_',num2str(1/high_freq_bound),'_',num2str(1/low_freq_bound),'_periods_',SSH_file_id,'_',num2str(2*altimetry_grid_res),'_spat_Gaussian_',num2str(time_range_start(1)),'_',num2str(time_range_end(1) - 1),'.pdf'],'-dpdf','-r300')
close(fig3)