close all
clear all

t_start   = tic;
file_path = pwd;                  % Change path if data not in local director
file_dat  = fullfile(file_path, 'CSIRO-320km_RCP26_R01_2009-2100.mat');
fprintf(1,'Loading data file "%s"\n', file_dat);
load(file_dat)

x_grid = map_western_edge-0.5*dx: dx : map_eastern_edge+0.5*dx;
y_grid = map_northern_edge+0.5*dy:-dy: map_southern_edge-0.5*dy;

figure(1)
imagesc(x_grid, y_grid, B), axis equal	
colorbar
xlabel('Easting distance (m)')
ylabel('Northing distance (m)')
title('Bed surface topography (m)')

% The year 2060 corresponds to the array element year(51). 
% To extract the surface topography at 2060 (year(51))

S_yr_map = squeeze(S(51,:,:));

figure(2)
imagesc(x_grid, y_grid, S_yr_map), axis equal
colorbar
xlabel('Easting distance (m)')
ylabel('Northing distance (m)')
title('Surface topography (m)')

H = S_yr_map-B;

figure(3)
imagesc(x_grid, y_grid, H), axis equal
colorbar
xlabel('Easting distance (m)')
ylabel('Northing distance (m)')
title('Ice thickness (m)')
	
I = zeros(size(H));
I(H>0) = 1;

figure(4)
imagesc(x_grid, y_grid, I), axis equal
colorbar
xlabel('Easting distance (m)')
ylabel('Northing distance (m)')
title('Ice mask')

figure(5)
imagesc(x_grid, y_grid, R_mask), axis equal
colorbar
xlabel('Easting distance (m)')
ylabel('Northing distance (m)')
title('Region mask')
  		
% Apply region mask to ice thickness data and plot masked result

H_masked = zeros(size(H));
H_masked(R_mask==1) = H(R_mask==1);

figure(6)
imagesc(x_grid, y_grid, H_masked), axis equal
colorbar
xlabel('Easting distance (m)')
ylabel('Northing distance (m)')
title('Region-masked ice thickness (m)')        

fprintf(1,'\nALL DONE. Elapsed time %.2f sec\n', toc(t_start));

