function rv = rvalue(h)
% function y = rvalue(h)
% ROMS stiffness parameter 
% Wilkin

[ m n ] = size(h);
if m < 2 | n < 2
  str = [ 'h must be a matrix but m = ' num2str(m) ' & n = ' ...
	num2str(n) ];
  disp(str)
end

% Check whether h contains any NaNs.  If it does then print out a
% message.
if any(isnan(h))
  warning('the height array contains NaNs')
end

% Check whether h changes sign.  If there are sign changes then store
% the relevant indices so that negative values of rv can later be stored
% in the proper positions.
minmax = nanrange(h);
if minmax(1)*minmax(2) <= 0
  str = [ 'Minimum h = ' num2str(minmax(1)) ' & maximum h = ' ...
	num2str(minmax(2)) ': these are different signs' ];
  warning(str)
end

% compute  diff(h)/2*mean(h) at each velocity grid point
dhdx_u = diff(h')';
dhdy_v = diff(h);
th_u = 2*av2(h')';
th_v = 2*av2(h);
r_u = abs(dhdx_u./th_u);
r_v = abs(dhdy_v./th_v);

% for each rho point, find the maximum rvalue at the 4 
% surrounding u,v points
n = size(r_u,2);
r_u = max(r_u(:,2:n),r_u(:,1:n-1));
n = size(r_v,1);
r_v = max(r_v(2:n,:),r_v(1:n-1,:));

% pad rows and columns to give a result the same size as the input h
r_u = r_u(:,[1 1:size(r_u,2) size(r_u,2)]);
r_v = r_v([1 1:size(r_v,1) size(r_v,1)],:);
rv = max(r_u,r_v);
