%** plot_slices.m
%
%   script to plot slices across the ice sheet

clear 
load gland-500-1-100.mat


%** preliminarires
%   select plotting indexes in row, column cooridnates
%       the default is to select one northern and one southern
%   select a starting figure number 
%   select time step to plot
%   

R1n=86;    % row # for slice through northern part of ice sheet
R2n=86;
C1n=20;    % column #
C2n=68;
R1s=65;      % row # for slice through southern part of ice sheet
R2s=65;
C1s=22;
C2s=64;

f=500;      % starting number for figures
n=6;        % time step


%******
% figures
%******

%* temperature slices
Tslice=T(n,:,R1n,C1n:C2n);
Tslice=double(squeeze(Tslice));

figure(f+11)
clf
cvT=[-20:2:0];
[cs, ha]=contourf(x1(C1n:C2n)*ones(1,length(zeta))/1e3, squeeze(b(n,R1n,C1n:C2n))*ones(1,length(zeta))+double(squeeze(h(n,R1n,C1n:C2n)))*(1-zeta'), Tslice', cvT);
%clabel(cs, ha, 'manual')
colorbar
ylabel('elevation (m)')
xlabel('distance (km)')
title('temperature in a slice through the ice sheet (N)')



Tslice=T(n,:,R1s,C1s:C2s);
Tslice=double(squeeze(Tslice));

figure(f+12)
clf
cvT=[-20:2:0];
[cs, ha]=contourf(x1(C1s:C2s)*ones(1,length(zeta))/1e3, squeeze(b(n,R1s,C1s:C2s))*ones(1,length(zeta))+double(squeeze(h(n,R1s,C1s:C2s)))*(1-zeta'), Tslice', cvT);
%clabel(cs, ha, 'manual')
colorbar
ylabel('elevation (m)')
xlabel('distance (km)')
title('temperature in a slice through the ice sheet (S)')




%* velocity slices
cv=[0:100:2000];
uxslice=double(squeeze(ux(n,:,R1n,C1n:C2n)));
uyslice=double(squeeze(ux(n,:,R1n,C1n:C2n)));
uzslice=double(squeeze(uz(n,:,R1n,C1n:C2n)));

figure(f+13)
clf
speed=sqrt(uxslice.^2+uyslice.^2);
[cs, ha]=contourf(x1(C1n:C2n)*ones(1,length(zeta))/1e3, squeeze(b(n,R1n,C1n:C2n))*ones(1,length(zeta))+double(squeeze(h(n,R1n,C1n:C2n)))*(1-zeta'), speed', cv);
colorbar
ylabel('elevation (m)')
xlabel('distance (km)')
title('horizontal speed in a slice through the ice sheet (N)')



uxslice=double(squeeze(ux(n,:,R1s,C1s:C2s)));
uyslice=double(squeeze(ux(n,:,R1s,C1s:C2s)));
uzslice=double(squeeze(uz(n,:,R1s,C1s:C2s)));

figure(f+14)
clf
speed=sqrt(uxslice.^2+uyslice.^2);
[cs, ha]=contourf(x1(C1s:C2s)*ones(1,length(zeta))/1e3, squeeze(b(n,R1s,C1s:C2s))*ones(1,length(zeta))+double(squeeze(h(n,R1s,C1s:C2s)))*(1-zeta'), speed', cv);
colorbar
ylabel('elevation (m)')
xlabel('distance (km)')
title('horizontal speed in a slice through the ice sheet (S)')