EKE_zonavg_tseries_POP.m

% zonally average EKE in POP output and produce time series at given latitude

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

% grid spacing (in degrees)
grid_res = 0.1;


% transect to be zonally averaged

tseries_plot_id = 'Atlantic_40N_tseries';
tseries_title_id = 'Atlantic at 40N';
lat_transect = 40.0;
lon_zonavg_bounds = [-95; 42];

% time period to include

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

SSH_file_id = 'POP_5day';

% 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);
n_files = 32;

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;
    curr_model_yr = 1978 + (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;

longwest_min = lon_zonavg_bounds(1);
longeast_max = lon_zonavg_bounds(2);
longwest_zonavg = lon_zonavg_bounds(1) + (360*ceil(((min(min(longitude)) - (grid_res/2)) - longwest_min)/360));
longeast_zonavg = lon_zonavg_bounds(2) + (360*ceil(((min(min(longitude)) - (grid_res/2)) - longwest_min)/360));


if longeast_zonavg > max(max(longitude)) + (grid_res/2)
    in_region_ind_1 = find((longitude >= longwest_zonavg - (5*grid_res)) & (latitude >= lat_transect - (10*grid_res)) & (latitude <= lat_transect + (10*grid_res)));
    in_region_i_ind_1 = unique(mod(in_region_ind_1 - 1,size(longitude,1)) + 1);
    in_region_j_ind_1 = unique(ceil(in_region_ind_1/size(longitude,1)));
    in_region_ind_2 = find((longitude <= longeast_zonavg + (5*grid_res) - 360) & (latitude >= lat_transect - (10*grid_res)) & (latitude <= lat_transect + (10*grid_res)));
    in_region_i_ind_2 = unique(mod(in_region_ind_2 - 1,size(longitude,1)) + 1);
    in_region_j_ind_2 = unique(ceil(in_region_ind_2/size(longitude,1)));
    in_lat_range_ind = unique([in_region_j_ind_1; in_region_j_ind_2]);
    curr_unique = unique([in_region_i_ind_1; in_region_i_ind_2]);
    [~,gap_ind] = max(diff(curr_unique));
    if max(diff(curr_unique)) > 1.5
        in_lon_range_ind = [(curr_unique(gap_ind + 1):1:size_array(1))'; (1:1:curr_unique(gap_ind))'];
    else
        in_lon_range_ind = (curr_unique(1):1:curr_unique(length(curr_unique)))';
    end
else
    in_region_ind = find((longitude >= longwest_zonavg - (5*grid_res)) & (longitude <= longeast_zonavg + (5*grid_res)) & (latitude >= lat_transect - (10*grid_res)) & (latitude <= lat_transect + (10*grid_res)));
    in_lon_range_ind = unique(mod(in_region_ind - 1,size(longitude,1)) + 1);
    in_lat_range_ind = unique(ceil(in_region_ind/size(longitude,1)));
end


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



lon_in_range = longitude(in_lon_range_ind,in_lat_range_ind);
lat_in_range = latitude(in_lon_range_ind,in_lat_range_ind);
    
% adjust longitude values to account for branch cut
lon_in_range(isnan(lon_in_range) == 1) = 1e4;
diff_lon = diff(lon_in_range,1,2);
diff_lon((abs(diff_lon) > 180) & (abs(diff_lon) < 5e3)) = mod(diff_lon((abs(diff_lon) > 180) & (abs(diff_lon) < 5e3)) + 180,360) - 180;
start_lon = lon_in_range(:,1);
start_lon(abs(start_lon) < 5e3) = start_lon(abs(start_lon) < 5e3) + (360*round((min(longwest_zonavg) - start_lon(abs(start_lon) < 5e3))/360));
lon_in_range = repmat(start_lon,[1 size(lon_in_range,2)]) + [zeros([size(lon_in_range,1) 1]) cumsum(diff_lon,2,'forward')];
lon_in_range(lon_in_range > 5e3) = 1e4;
diff_lon = diff(lon_in_range,1,1);
diff_lon((abs(diff_lon) > 180) & (abs(diff_lon) < 5e3)) = mod(diff_lon((abs(diff_lon) > 180) & (abs(diff_lon) < 5e3)) + 180,360) - 180;
start_lon = lon_in_range(1,:);
start_lon(abs(start_lon) < 5e3) = start_lon(abs(start_lon) < 5e3) + (360*round((min(longwest_zonavg) - start_lon(abs(start_lon) < 5e3))/360));
lon_in_range = repmat(start_lon,[size(lon_in_range,1) 1]) + [zeros([1 size(lon_in_range,2)]); cumsum(diff_lon,1,'forward')];
lon_in_range(lon_in_range > 5e3) = NaN;


% compute SSH mean over all available times

size_array = [length(in_lon_range_ind) length(in_lat_range_ind) 1 length(time_datenum)];
% delta_t = mean(diff(time_datenum));

SSH_in_range = NaN(size_array([1 2 4]));
for file_num_ind = min(file_num_vec):max(file_num_vec)
%     curr_model_yr = time_range_start(1) + (file_num_ind - 1) - year_offset;
    curr_model_yr = 1978 + (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_vec == file_num_ind);
    if isempty(curr_time_subset_ind) == 1
        continue
    end
    curr_in_time_range_ind = in_file_ind_vec(curr_time_subset_ind);
    time_range_ind_span = max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1;
    
    lat_range_ind_span = max(in_lat_range_ind) - min(in_lat_range_ind) + 1;
    
    if min(diff(in_lon_range_ind)) < -0.5
        lon_range_ind_span_1 = max(in_lon_range_ind) - in_lon_range_ind(1) + 1;
        start_vec = [in_lon_range_ind(1) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
        count_vec = [lon_range_ind_span_1 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]))));
        curr_SSH_in_range(abs(curr_SSH_in_range) > 1e15) = 0;
        SSH_in_range(1:1:count_vec(1),:,curr_time_subset_ind) = curr_SSH_in_range(:,:,curr_in_time_range_ind - min(curr_in_time_range_ind) + 1);
        lon_range_ind_span_2 = in_lon_range_ind(length(in_lon_range_ind)) - min(in_lon_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_2 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]))));
        curr_SSH_in_range(abs(curr_SSH_in_range) > 1e15) = 0;
        SSH_in_range(lon_range_ind_span_1 + (1:1:count_vec(1)),:,curr_time_subset_ind) = curr_SSH_in_range(:,:,curr_in_time_range_ind - min(curr_in_time_range_ind) + 1);
    else
        lon_range_ind_span = max(in_lon_range_ind) - min(in_lon_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]))));
        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


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;
steepness_factor_spatfilt = 5;
waven_high_bound = (1/(2*altimetry_grid_res))*(exp(erfinv((2^(1/2)) - 1)/steepness_factor_spatfilt));
waven_x_low_bound = 1/(2*delta_x*size(SSH_in_range,1));
waven_y_low_bound = 1/(4*delta_y*size(SSH_in_range,2));


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 = find(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)));
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);
    if diff(lon_zonavg_bounds) > 358
        [curr_SSH_in_range_x_bp,~,~] = bandpass_err_fcn_no_pad([SSH_in_range(:,:,curr_t_min:curr_t_max) NaN([size_array(1:2) (curr_t_max - curr_t_min + 1)])],1,delta_x,waven_x_low_bound,waven_high_bound,steepness_factor_spatfilt,0,1,1,0);
    else
        [curr_SSH_in_range_x_bp,~,~] = bandpass_err_fcn_no_pad([SSH_in_range(:,:,curr_t_min:curr_t_max) NaN([size_array(1:2) (curr_t_max - curr_t_min + 1)])],1,delta_x,waven_x_low_bound,waven_high_bound,steepness_factor_spatfilt,1,1,1,0);
    end
    [curr_SSH_in_range_lp_expanded,~,~] = bandpass_err_fcn_no_pad(curr_SSH_in_range_x_bp,2,delta_y,waven_y_low_bound,waven_high_bound,steepness_factor_spatfilt,1,1,1,0);

    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

if diff(lon_zonavg_bounds) > 358
    [nan_2d_x_bp,~,~] = bandpass_err_fcn_no_pad([nan_mask_array_2d ((-1)*ones(size_array(1:2)))],1,delta_x,1/(2*delta_x*size_array(1)),waven_high_bound,steepness_factor_spatfilt,0,1,1,0);
else
    [nan_2d_x_bp,~,~] = bandpass_err_fcn_no_pad([nan_mask_array_2d ((-1)*ones(size_array(1:2)))],1,delta_x,1/(2*delta_x*size_array(1)),waven_high_bound,steepness_factor_spatfilt,1,1,1,0);
end
[nan_2d_spat_bp_expanded,~,~] = bandpass_err_fcn_no_pad(nan_2d_x_bp,2,delta_y,1/(2*2*delta_y*size_array(2)),waven_high_bound,steepness_factor_spatfilt,1,1,1,0);

SSH_in_range_lp = SSH_in_range_lp_expanded(:,1:1:size_array(2),:);
nan_2d_spat_bp = nan_2d_spat_bp_expanded(:,1:1:size_array(2));

clear *_expanded

SSH_in_range_lp(repmat(abs(nan_2d_spat_bp),[1 1 size_array(3)]) > 0.5) = NaN;

SSH_in_range = SSH_in_range_lp;
clear SSH_in_range_lp


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_all = sum(nan_mask_SSH.*SSH_in_range_nans_zeroed,3)./sum(nan_mask_SSH,3);



% compute EKE within specified time range

time_datenum_in_range = time_datenum(in_time_range_ind);
time_bound_datenum_in_range = time_bound(in_time_range_ind,:);
file_num_in_range_vec = file_num_vec(in_time_range_ind);
in_file_ind_in_range_vec = 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)];
% delta_t = mean(diff(time_datenum_in_range));

SSH_in_range = NaN(size_array([1 2 4]));
for file_num_ind = min(file_num_vec):max(file_num_vec)
%     curr_model_yr = time_range_start(1) + (file_num_ind - 1) - year_offset;
    curr_model_yr = 1978 + (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);
    time_range_ind_span = max(curr_in_time_range_ind) - min(curr_in_time_range_ind) + 1;
    
    lat_range_ind_span = max(in_lat_range_ind) - min(in_lat_range_ind) + 1;
    
    if min(diff(in_lon_range_ind)) < -0.5
        lon_range_ind_span_1 = max(in_lon_range_ind) - in_lon_range_ind(1) + 1;
        start_vec = [in_lon_range_ind(1) min(in_lat_range_ind) 1 min(curr_in_time_range_ind)];
        count_vec = [lon_range_ind_span_1 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]))));
        curr_SSH_in_range(abs(curr_SSH_in_range) > 1e15) = 0;
        SSH_in_range(1:1:count_vec(1),:,curr_time_subset_ind) = curr_SSH_in_range(:,:,curr_in_time_range_ind - min(curr_in_time_range_ind) + 1);
        lon_range_ind_span_2 = in_lon_range_ind(length(in_lon_range_ind)) - min(in_lon_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_2 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]))));
        curr_SSH_in_range(abs(curr_SSH_in_range) > 1e15) = 0;
        SSH_in_range(lon_range_ind_span_1 + (1:1:count_vec(1)),:,curr_time_subset_ind) = curr_SSH_in_range(:,:,curr_in_time_range_ind - min(curr_in_time_range_ind) + 1);
    else
        lon_range_ind_span = max(in_lon_range_ind) - min(in_lon_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]))));
        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


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;
steepness_factor_spatfilt = 5;
waven_high_bound = (1/(2*altimetry_grid_res))*(exp(erfinv((2^(1/2)) - 1)/steepness_factor_spatfilt));
waven_x_low_bound = 1/(2*delta_x*size(SSH_in_range,1));
waven_y_low_bound = 1/(4*delta_y*size(SSH_in_range,2));


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 = find(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)));
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_SSH_in_range_x_bp,~,~] = bandpass_err_fcn_no_pad([SSH_in_range(:,:,curr_t_min:curr_t_max) flip(SSH_in_range([(((size_array(1)/2) + 1):1:size_array(1))'; (1:1:(size_array(1)/2))'],:,curr_t_min:curr_t_max),2)],1,delta_x,waven_x_low_bound,waven_high_bound,steepness_factor_spatfilt,0,1,1,0);

    [curr_SSH_in_range_lp_expanded,~,~] = bandpass_err_fcn_no_pad(curr_SSH_in_range_x_bp,2,delta_y,waven_y_low_bound,waven_high_bound,steepness_factor_spatfilt,0,1,1,0);

    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

[nan_2d_x_bp,~,~] = bandpass_err_fcn_no_pad([nan_mask_array_2d flip(nan_mask_array_2d([(((size_array(1)/2) + 1):1:size_array(1))'; (1:1:(size_array(1)/2))'],:),2)],1,delta_x,1/(2*delta_x*size(nan_mask_array_2d,1)),waven_high_bound,steepness_factor_spatfilt,0,1,1,0);
[nan_2d_spat_bp_expanded,~,~] = bandpass_err_fcn_no_pad(nan_2d_x_bp,2,delta_y,waven_y_low_bound,1/(2*delta_y*size(nan_2d_x_bp,2)),steepness_factor_spatfilt,0,1,1,0);

SSH_in_range_lp = SSH_in_range_lp_expanded(:,1:1:size_array(2),:);
nan_2d_spat_bp = nan_2d_spat_bp_expanded(:,1:1:size_array(2));

clear *_expanded

SSH_in_range_lp(repmat(abs(nan_2d_spat_bp),[1 1 size_array(3)]) > 0.5) = NaN;

SSH_in_range = SSH_in_range_lp;
clear SSH_in_range_lp


SSH_nomean_in_range = SSH_in_range - repmat(SSH_mean_all,[1 1 size(SSH_in_range,3)]);


lon_in_range_midpts = lon_in_range(2:size_array(1),:) - (diff(lon_in_range,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) 1])]; NaN([1 size_array(2)])];
lat_in_range_midpts = lat_in_range(2:size_array(1),:) - (diff(lat_in_range,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) 1])]; NaN([1 size_array(2)])];
dx = 111100*[(abs((cosd(lat_in_range_midpts(1:(size_array(1) - 1),:)).*diff([(lon_in_range(:,2:size_array(2)) - (diff(lon_in_range,1,2)/2)) NaN([size_array(1) 1])],1,1)) + (1i*diff([(lat_in_range(:,2:size_array(2)) - (diff(lat_in_range,1,2)/2)) NaN([size_array(1) 1])],1,1)))); NaN([1 size_array(2)])];
dy = 111100*[(abs((cosd(lat_in_range_midpts(:,1:(size_array(2) - 1))).*diff([(lon_in_range(2:size_array(1),:) - (diff(lon_in_range,1,1)/2)); NaN([1 size_array(2)])],1,2)) + (1i*diff([(lat_in_range(2:size_array(1),:) - (diff(lat_in_range,1,1)/2)); NaN([1 size_array(2)])],1,2)))) NaN([size_array(1) 1])];

diff_x_SSH_nomean_in_r = diff(SSH_nomean_in_range,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([1 (size_array(2) - 1) size_array(3)])] 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:size(diff_y_SSH_nomean_in_r,1),:,:) - (diff(diff_y_SSH_nomean_in_r,1,1)/2)) NaN([(size_array(1) - 1) 1 size_array(3)])]; NaN([1 size_array(2: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(abs(lat_mask_3d) < 1e-5) = NaN;


% compute zonally-averaged time series

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;
sum_nan_mask = sum(nan_mask_EKE,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_EKE,3)]);
nan_mask_EKE = nan_mask_EKE & nan_mask_not_enough_3D;

% EKE_mean_in_range = sum(nan_mask_EKE.*EKE_in_range_nans_zeroed,3)./sum(nan_mask_EKE,3);

curr_lat_bin_bounds = lat_transect + (2*grid_res*[-1 1]);
    
in_curr_bin_ind = find((lon_in_range_midpts >= longwest_zonavg) & (lon_in_range_midpts <= longeast_zonavg) & (lat_in_range_midpts >= curr_lat_bin_bounds(1)) & (lat_in_range_midpts <= curr_lat_bin_bounds(2)));
in_region_mask = zeros(size(lon_in_range_midpts));
in_region_mask(in_curr_bin_ind) = 1;
    
EKE_in_range_zonavg_POP = squeeze(sum(sum(repmat(in_region_mask,[1 1 size(nan_mask_EKE,3)]).*nan_mask_EKE.*EKE_in_range_nans_zeroed,2),1)./(sum(sum(repmat(in_region_mask,[1 1 size_array(3)]).*nan_mask_EKE,2),1)));

save(['EKE_zonavg_',SSH_file_id,'_',tseries_plot_id,'_',num2str(time_range_start(1)),'_',num2str(time_range_end(1) - 1),'.mat'],'time_datenum_in_range','EKE_in_range_zonavg_POP')
% load(['EKE_zonavg_',SSH_file_id,'_',tseries_plot_id,'_',num2str(time_range_start(1)),'_',num2str(time_range_end(1) - 1),'.mat'],'time_datenum_in_range','EKE_in_range_zonavg_POP')


% % load altimetry EKE averages for comparison
% load(['../SSH/EKE_tmean_zonavg_SSALTO_DUACS_',region_id,'_1993_2016.mat'])
% lat_in_range_midpts_alt = lat_in_range_midpts;
% EKE_mean_in_range_zonavg_alt = EKE_mean_in_range_zonavg;


% apply temporal filters to time series

t_seasoncyc_opt = 0;    % 0 = remove seasonal cycle, 1 = retain seasonal cycle

t_bandpass_opt = 1;    % 0 = no bandpassing in time, 1 = bandpassing in time
low_freq_bound = 1/(4*mean(diff(time_datenum_in_range))*length(time_datenum_in_range));
high_freq_bound = 1/426;
half_power_adj = exp(erfinv((2^(1/2)) - 1)/5);   % adjustment factor to set bounds at half-power (rather than half-amplitude)


if abs(t_seasoncyc_opt) < 1e-5
    
    % remove seasonal cycle
    
%     nan_mask_zonavg = ones(size(EKE_in_range_zonavg_POP));
%     nan_mask_zonavg((isnan(EKE_in_range_zonavg_POP) == 1) | (abs(EKE_in_range_zonavg_POP) < 1e-5)) = 0;
%     EKE_in_range_zonavg_POP(abs(nan_mask_zonavg) < 1e-5) = 0;
    
    n_bins = round(365/delta_t);
    EKE_in_range_zonavg_POP_noseason = NaN(size(EKE_in_range_zonavg_POP));
    for bin_ind = 1:n_bins
        curr_in_bin_ind = find(abs(mod(time_datenum_in_range - time_datenum_in_range(bin_ind) + (365/2),365) - (365/2)) < (0.4*delta_t));
        good_ind_in_bin = find(isnan(EKE_in_range_zonavg_POP(curr_in_bin_ind)) == 0);
        if length(good_ind_in_bin) > 0.8*length(curr_in_bin_ind)
            EKE_in_range_zonavg_POP_noseason(curr_in_bin_ind) = EKE_in_range_zonavg_POP(curr_in_bin_ind) - (sum(EKE_in_range_zonavg_POP(curr_in_bin_ind(good_ind_in_bin)),1)/(length(good_ind_in_bin)));
        end
    end
    
    EKE_in_range_zonavg_POP_noseason = EKE_in_range_zonavg_POP_noseason + mean(EKE_in_range_zonavg_POP(isnan(EKE_in_range_zonavg_POP) == 0));
    
    EKE_in_range_zonavg_POP_filt = EKE_in_range_zonavg_POP_noseason;
    
    EKE_in_range_zonavg_POP_filt(isnan(EKE_in_range_zonavg_POP_filt) == 1) = 0;
    
    clear *_noseason
    
else
    EKE_in_range_zonavg_POP_filt = EKE_in_range_zonavg_POP;
end

if abs(t_bandpass_opt - 1) < 1e-5
    % filter time series for interannual & decadal frequencies
    [EKE_in_range_zonavg_POP_bpfilt,~,~] = bandpass_err_fcn(EKE_in_range_zonavg_POP_filt,1,delta_t,low_freq_bound/half_power_adj,high_freq_bound*half_power_adj,5,1,1,1,0);
    EKE_in_range_zonavg_POP_filt = EKE_in_range_zonavg_POP_bpfilt;
    clear *_bpfilt
end


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


% plot zonally-averaged time series of EKE

fig2 = figure(2);
% h = plot(time_datenum_in_range,(1e4)*EKE_in_range_zonavg_POP_filt);
yyaxis left
h1 = plot(time_datenum_in_range,(1e4)*EKE_in_range_zonavg_POP_filt);
set(gca,'YColor',[0 0 0.8])
set(h1,'Color',[0 0 0.8],'LineWidth',2)
set(get(gca,'ylabel'),'String','EKE (cm ^2 s^{-2})','FontSize',14)
yyaxis right
% h2 = plot(time_datenum_in_range,(1e-15)*T_flux_hpcross10_halfamp);
h2 = plot(time_datenum_in_range,(1e-15)*T_flux_mesoscale_fullbasin);
set(gca,'YColor',[0.8 0 0])
set(h2,'Color',[0.8 0 0],'LineWidth',2)
set(get(gca,'ylabel'),'String','Mesoscale heat flux (PW)','FontSize',14)
set(gca,'FontSize',14)
% set(h(1),'Color',[0 0 0.8],'LineStyle','-','LineWidth',2)
% % set(h(2),'Color',[0.8 0 0],'LineStyle','-','LineWidth',2)
% set(get(gca,'ylabel'),'String','EKE (cm ^2 s^{-2})','FontSize',14)
years_to_plot_ticks = ((ceil(time_range_start(1)/5)*5):5:(ceil(time_range_end(1)/5)*5))';
xtick_datenums_plot = 365*(years_to_plot_ticks - year_offset);
xtick_labels_plot = cell(length(years_to_plot_ticks),1);
for xtick_ind = 1:length(years_to_plot_ticks)
%     xtick_labels_plot{xtick_ind} = [num2str(years_to_plot_ticks(xtick_ind)),'-01-01'];
    xtick_labels_plot{xtick_ind} = num2str(years_to_plot_ticks(xtick_ind));
end
% n_ticks = (3*(time_range_end(1) - time_range_start(1) + 1)) + 1;
% xtick_datenums_plot = (365*time_range_start(1)) + (datenum([(1990*ones([n_ticks 1])) (1:4:(((n_ticks - 1)*4) + 1))' ones([n_ticks 1]) zeros([n_ticks 3])]) - datenum([1990 1 1 0 0 0])) - ceil((datenum([(1990*ones([n_ticks 1])) (1:4:(((n_ticks - 1)*4) + 1))' ones([n_ticks 1]) zeros([n_ticks 3])]) - datenum([1992 2 29 0 0 0]))/(4*365.25)); 
% xtick_labels_plot = cell(length(xtick_datenums_plot),1);
% for xtick_ind = 1:length(xtick_datenums_plot)
%     xtick_labels_plot{xtick_ind} = ['00',num2str(floor(xtick_datenums_plot(xtick_ind)/365)),datestr([1990 01 mod(xtick_datenums_plot(xtick_ind),365) 0 0 0],'-mm-dd')];
% end
set(gca,'FontSize',12,'xlim',[datenum_start datenum_end],'xtick',xtick_datenums_plot,'xticklabel',xtick_labels_plot,'xgrid','on')
hold on
% % line([min(lat_band_bounds) max(lat_band_bounds)],[0 0],[0 0],'Color',[0 0 0],'LineWidth',1,'LineStyle','-')
% line([datenum_start datenum_end],[0 0],[0 0],'Color',[0 0 0],'LineWidth',1,'LineStyle','-')
hold off
% title({'Zonally-averaged EKE,'; [tseries_title_id,', ',num2str(time_range_start(1)),'-',num2str(time_range_end(1) - 1),' POP']; ' '},'FontSize',8)
title({'Zonally-averaged EKE and mesoscale heat flux,'; [tseries_title_id,', ',num2str(time_range_start(1)),'-',num2str(time_range_end(1) - 1),' POP']; ' '},'FontSize',8)
% leg = legend('Zonally-averaged EKE',' ','location','NorthEast');

% print(fig2,['EKE_zonavg_POP_',tseries_plot_id,'_',num2str(time_range_start(1)),'_',num2str(time_range_end(1) - 1),'.pdf'],'-dpdf','-r300')
print(fig2,['EKE_zonavg_meso_hflux_POP_',tseries_plot_id,'_',num2str(time_range_start(1)),'_',num2str(time_range_end(1) - 1),'.pdf'],'-dpdf','-r300')
close(fig2)