% oceanrise.m
%
% ess345 problem 3: calculate water volume in oceans
% and amount of sea level rise if Antarctic and Greenland ice melts
% this is my second attempt (first script is lost, buh)
% CGB, 25 Feb 2014
 clear all; close all; format compact;
 
% 1. load elevation data (1x1 degree Matlab internal file)
%    show simple world map (no specific projection), add coastline
 load topo.mat
 lat=[-89.5:1:89.5]; lon=[0.5:1:359.5];
 figure; pcolor(topo); axis image; shading flat; colorbar;
 hold on;
 %[Lon,Lat]=meshgrid(lon,lat);
 contour(topo,[0 0],'k')
 
 
% 2. determine gridpoints with negative elevations (ie, oceans)
 [iy,ix]=find(topo>0);
 %plot(ix,iy,'m.');
 
% 3. create new elevation array, with continents at 0 elevation
 ocean=topo;
 for ii=1:size(ix,1)
    ocean(iy(ii),ix(ii))=0; 
 end
 figure; pcolor(ocean); shading flat; axis image; colorbar;
 title('ocean depth')
 
% 4. calculate volume 
%    a) surface area of 1x1 degree grid
 R = 6371;           % radius of Earth in km
 dlon = 2*pi*R/360;  % distance along 1 degree of meridian
 for indx=1:size(lon,2)
    for indy=1:size(lat,2)
       dlat=lat(indy);
       pix(indy,indx) = dlon * 2*pi*R*cosd(lat(indy))/360;
    end
 end
%    b) volume at each grid point is area times depth
 vol = -.001*ocean .* pix;
%    c) total volume
 VOL = sum(sum(vol))
 
 
% 5. volume of meltwater
%    a) Antarctica above sea level (blow 60 deg S)
 ant=topo(1:30,:);
 [ixA,iyA]=find(ant < 0);
 for ii=1:size(ixA,1)
    ant(ixA(ii),iyA(ii))=0; 
 end
 %figure; pcolor(ant); shading flat; axis image; colorbar;
 volA = ant/1000 .* pix(1:30,:);
 VolA = sum(sum(volA)) * 0.9   % meltwater has about 90% of original ice volume
%    b) Greenland above sea level 
 figure(1)
 %[xx,yy]=ginput   % selected boundary of Greenland by hand
 xx =[ 297.3443 287.1840 303.6341 321.0519 352.9845];
 yy =[ 177.1051 166.4609 157.2682 139.3666 175.6536];
 xx=[xx,xx(1)]; yy=[yy,yy(1)];   % to close polygon
 plot(xx,yy,'r')
 [X,Y]=meshgrid([1:360],[1:180]);
 in = inpolygon(X,Y,xx,yy);
 plot(X(in),Y(in),'m.')
 eleG = in.*topo;
 [ixG,iyG]=find(eleG < 0);
 for ii=1:size(ixG,1)
    eleG(ixG(ii),iyG(ii))=0; 
 end
 volG = eleG/1000 .* pix;
 VolG = sum(sum(volG)) * 0.9
%  c) total volume
 meltwater = VolG + VolA
 
% 6. now need to determine rise of ocean
%    a) check by how much ocean volume would increase
 increase = 10:10:160;  % 10 meter steps
 for ind=1:size(increase,2)
     ocean=topo-increase(ind);
     [iy,ix]=find(topo>increase(ind));
     for ii=1:size(ix,1)
         ocean(iy(ii),ix(ii))=0; 
     end
     %figure; pcolor(ocean); shading flat; axis image; colorbar
%    b) volume at each grid point is area times depth
     vol = -.001*ocean .* pix;
%    c) total volume
     VOLm(ind) = sum(sum(vol)) - VOL;
 end
%    d) show result; estimated rise is 80 m 
 figure; plot(increase,VOLm,'ko-'); hold on; grid on;
 plot([0 max(increase)],[meltwater meltwater],'r')
 xlabel('sea level rise [m]')
 ylabel('volume of meltwater added to oceans [km^3]','interpreter','none')
 title('sea-level rise (red: expected volume)')
%    e) show now coastline on a world map (equal-area projection) 
 figure; pcolor(lon,lat,topo); shading flat; axis image; 
 colormap('gray'), colorbar; hold on
 contourf(lon,lat,topo,[ 0 0],'b')
 contour(lon,lat,topo,[75 75],'r')
 title('current (blue) and projected (red) coastline')