function [regions,overlappedreg] = gensubregions(ax_range,overlap,nregions)
% Generate a set of subregion tiles to cover the input range of lat/lon
%
% [regions,overlappedregions] = gensubregions(ax_range,overlap,nregions)
%
% Input:
%       AX_RANGE [west east south north] of entire domain to be tiled
%       OVERLAP  in degrees lat/lon
%       NREGIONS [nlon nlat] number of subregions in each direction
%
% Output:
%       REGIONS matrix of region boundaries with each row in format 
%               [west east south north overlap] where column 5 is simply the
%               input OVERLAP value
%       OVERLAPPEDREGIONS (optional) as above but with overlap added (or
%       subtracted) to the limits
%
% John Wilkin 22 Jul 97

xl = linspace(ax_range(1),ax_range(2),nregions(1)+1);
yl = linspace(ax_range(3),ax_range(4),nregions(2)+1);
[x,y] = meshgrid(xl,yl);

regions = [];
if nargout == 2
  overlappedreg = [];
end
for i=1:nregions(2)
  for j=1:nregions(1)
    regions = [regions; ...
	[x(i,j) x(i+1,j+1) y(i,j) y(i+1,j+1) overlap]];
  end
end

if nargout == 2
  add = overlap*[-1 1 -1 1 0];
  add = add(ones([nregions 1]),:);
  overlappedreg = regions + add;
end