function y = rv(h, resc_r, h_change)
% function y = rv(h, resc_r, h_change)

% Given a matrix of bathymetric heights return the corresponding matrix
% of stability parameters rv as used in numerical models - notably
% spem.  For depths less than h_change rv is as discussed in the spem
% manual and is approximately half of the slope divided by the mean
% height.  This slope factor is discussed by Mellor et al and rv = 0.2
% corresponds to their slope factor = 0.4.  For the deeper oceans rv is
% taken to be proportional to the cube of the difference between h at
% adjacent grid points divided by the mean depth.  The constant of
% proportionality is chosen so that rv is continuous if there is
% marginal stability (i.e., rv = r_stable) at h = h_change.  The
% variable resc_r = 1./(2*r_stable*h_change).^2 is used to calculate
% this proportionality constant and so is only relevant if h_change is
% less than the greatest ocean depth.

% h must be a matrix and the returned matrix, y, is of the same size.
% To generate y at any point find the difference between h at that point
% and h at the point with the same column index but row index one less.
% Then divide by the sum of h at these two points and take the absolute
% value.  Now do the same again but interchange the roles of the rows and
% columns.  The returned value at the point is the maximum of of these
% two calculated values.  If h is the bottom height in a numerical model
% then the value of y at a given point is a stability parameter; it is
% approximately half of the slope divided by the mean height.  Note that
% this parameter is meaningless if h changes sign anywhere.  Accordingly,
% I check for that and in the relevant positions a very large positive
% value is returned.
% I have changed rv so that NaNs in the height field do not contribute
% towards the stability.  This allows me to mark islands with NaNs and
% not have them affect the stability calculations.

% Copyright J. V. Mansbridge, CSIRO, august 17 1992
%$Id: rv.m,v 1.4 1995/05/15 09:12:32 mansbrid Exp $

% First check that h is a matrix.

[ 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))
  disp('rv: the height array contains NaNs')
%  return
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 = [ 'rv: Minimum h = ' num2str(minmax(1)) ' & maximum h = ' ...
      num2str(minmax(2)) ': these are different signs' ];
  disp(str)
%  return
  tmp = ones(m, n);
  tmp(2:m, :) = h(2:m, :).*h(1:(m-1), :);
  index_x = find ( tmp <= 0 );
  tmp(:, 2:n) = h(:, 2:n).*h(:, 1:(n-1));
  index_y = find ( tmp <= 0 );
  set_neg_1 = 1;
else
  set_neg_1 = 0;
end

% Find the normalised slopes in each direction.

%dhdxx = abs((h(2:m, :) - h(1:(m-1), :))./(h(2:m, :) + h(1:(m-1), :)));
%dhdyy = abs((h(:, 2:n) - h(:, 1:(n-1)))./(h(:, 2:n) + h(:, 1:(n-1))));
%dhdxx = abs( diff(h)./(h(2:m, :) + h(1:(m-1), :)) );
%dhdyy = abs( diff(h')'./(h(:, 2:n) + h(:, 1:(n-1))) );

diff_x = diff(h);
h_sum_x = h(2:m, :) + h(1:(m-1), :);
diff_y = diff(h')';
h_sum_y = h(:, 2:n) + h(:, 1:(n-1));

bad_x = find(h_sum_x == 0);
if length(bad_x) > 0
  disp(['There are ' num2str(length(bad_x)) ' x shoreline points'])
  h_sum_x(bad_x) = -0.1*ones(size(bad_x));
  set_neg_x = 1;
else
  set_neg_x = 0;
end

bad_y = find(h_sum_y == 0);
if length(bad_y) > 0
  disp(['There are ' num2str(length(bad_y)) ' y shoreline points'])
  h_sum_y(bad_y) = -0.1*ones(size(bad_y));
  set_neg_y = 1;
else
  set_neg_y = 0;
end

dhdxx = abs(diff_x./h_sum_x);
dhdxx_deep = resc_r*abs((diff_x.^3)./h_sum_x);
ff = find(h_sum_x > 2*h_change);
if length(ff) > 0
  dhdxx(ff) = dhdxx_deep(ff);
end

dhdyy = abs(diff_y./h_sum_y);
dhdyy_deep = resc_r*abs((diff_y.^3)./h_sum_y);
ff = find(h_sum_y > 2*h_change);
if length(ff) > 0
  dhdyy(ff) = dhdyy_deep(ff);
end

% Find the matrix in which each value is the maximum of the 2
% corresponding values in dhdxx and dhdyy.

%y = nanmax(dhdxx, dhdyy);
y1 = nanmax( dhdxx(:, 2:n), dhdxx(:, 1:(n-1)) );
y2 = nanmax( dhdyy(2:m, :), dhdyy(1:(m-1), :) );
y = nanmax(y1, y2);
y = [zeros(m-1,1) y];
y = [zeros(1, n); y];

%Set the unacceptable stability values of points that are above the
%water or so close to the shoreline that they have a mean depth of zero.

if set_neg_1
  if length(index_x) > 0
    y(index_x) = 10000*ones(size(index_x));
  end
  if length(index_y) > 0
    y(index_y) = 10000*ones(size(index_y));
  end
end

if set_neg_x
  y(bad_x) = 10000*ones(size(bad_x));
end

if set_neg_y
  y(bad_y) = 10000*ones(size(bad_y));
end
