%** this script produces maps of earth surface elevation for the present
%   and as reconstructed for the Last Glacial Maximum 
%
% globaltopo.mat    Earth surface elevation, meters
%   topoE       present-day 180 x 360 matrix, 1 degree grid
%   geoidE      present-day geoid height  180 x 360 matrix, 1 degree grid
%   IceAge      last glacial maximum as computed by Peltier (ICEG4)
%               180 x 360 matrix, 1 degree grid
%   maplat maplon:    latitude and longitude, degrees
%
% reference
% Peltier, W. R., Ice Age Paleotopography, 1994, Science, v. 265, 
%   pp. 195-201.


load globaltopo
load coast                          % Matlab native file

wind=[min(maplon) max(maplon) min(maplat) max(maplat)];     % plot window 
Cmint=-9000;                        % minimum surface elevation for color
Cmaxt=7000;                         % maximum surface elevation


%******
%  Earth surace elevation
%******

figure(1)
clf
pcolor(maplon, maplat, topoE)
shading flat
axis equal 
axis(wind)
hold on
caxis([Cmint Cmaxt])
plot(long, lat, 'k-')
colorbar('horiz')
title('Earth surface elevation (m) on a 1x1 grid')
xlabel('longitude')
ylabel('latitude')


figure(2)
clf
pcolor(maplon, maplat, IceAge)
shading flat
axis equal
axis(wind)
hold on
caxis([Cmint Cmaxt])
plot(long, lat, 'k-')
colorbar('horiz')
title('Peltier Earth surface elevation at LGM (m) on a 1x1 grid')
xlabel('longitude')
ylabel('latitude')


figure(3)
clf
pcolor(maplon, maplat, IceAge-topoE)
shading flat
axis equal
axis(wind)
hold on
plot(long, lat, 'k-')
colorbar('horiz')
title('change in surface elevation at LGM to present-day (m)')
xlabel('longitude')
ylabel('latitude')


figure(4)
clf
pcolor(maplon, maplat, geoidE)
shading flat
axis equal 
axis(wind)
hold on
title('geoid height (m)')
colorbar('horiz')
hold on
plot(long, lat, 'k-')