function [h] = paintmask(twelfths,height)
% [h] = paintmask(twelfths,height)
% paints a land mask on the existing axes
%
% Optional input TWELFTHS is the resolution of the mask in 12ths of a
% degree.  The default is TWELFTHS=4 => 4/12 = 1/3 degree lat/lon.
% 
% Optional input HEIGHT determines the land/sea discriminating altitude, 
% in metres above sea level, and of course defaults to 0. 
% eg height = -200  masks above the 200m isobath.
%
% The mask is created by reading the etopo5 bathymetry file using 
% function LANDMASK. 
%
% If the axis limits and TWELFTHS don't change, then PAINTMASK will 
% use the mask created on a previous call (stored as global variables).
% In this case it's quite fast.  It is reading and creating the mask that 
% is slow.  
% 
% Optional output H is the handle to the mask object generated by PATCH.
%
% This function should work happily in conjunction with a pcolor plot.
%
% Note that ETOPO5, the source of this data, is flaky; for example, most of
% Timor and all of Sumba are recorded as being deep ocean!
%
% John Wilkin 30 Oct 96

if nargin < 2
  height = 0;
end
if nargin == 0
  % default resolution is 3/12 = 1/4 degree
  twelfths = 3;
end

global paintmaskXX paintmaskYY paintmaskXlim paintmaskYlim paintmasktwelfths
global paintmaskheight

xlim = get(gca,'xlim');
ylim = get(gca,'ylim');

if isempty(paintmaskXX) ...
    | max([abs(xlim-paintmaskXlim) abs(ylim-paintmaskYlim)]) ~= 0 ...
    | paintmasktwelfths ~= twelfths ...
    | paintmaskheight ~= height
 
  [mask,lon,lat] = landmask(xlim,ylim,twelfths,1,NaN,height);
  [lon,lat] = meshgrid(lon,lat);

  water = find(isnan(mask));
  lon(water) = [];
  lat(water) = [];
  
  dx = twelfths/24;
  dy = twelfths/24;
  lon = lon(:)';
  lat = lat(:)';
  ll = length(lon);
  
  paintmaskXX = [-dx*ones(1,ll); dx*ones(2,ll); -dx*ones(1,ll)] + ...
      [lon;lon;lon;lon];
  paintmaskYY = [-dy*ones(2,ll); dy*ones(2,ll)] + [lat;lat;lat;lat];

  paintmaskXlim = xlim;
  paintmaskYlim = ylim;

  paintmasktwelfths = twelfths;
  paintmaskheight = height;
end

han = patch(paintmaskXX,paintmaskYY,'w');
set(han,'edgecolor','none','facecolor',0.8*[1 1 1])

if nargout > 0
  h = han;
end
t