function [r]=rfactor(h,rmask);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2001 Rutgers University.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
% function [r]=rfactor(h,rmask)                                             %
%                                                                           %
% This function computes the bathymetry slope from a SCRUM NetCDF file.     %
%                                                                           %
% On Input:                                                                 %
%                                                                           %
%    h           bathymetry at RHO-points.                                  %
%    rmask       Land/Sea masking at RHO-points.                            %
%                                                                           %
% On Output:                                                                %
%                                                                           %
%    r           R-factor.                                                  %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Lp,Mp]=size(h);
L=Lp-1;
M=Mp-1;

if (nargin < 2),
  rmask=ones(size(h));
end,

%  Land/Sea mask on U-points.

for j=1:Mp,
  for i=2:Lp,
    umask(i-1,j)=rmask(i,j)*rmask(i-1,j);
  end,
end,

%  Land/Sea mask on V-points.

for j=2:Mp,
  for i=1:Lp,
    vmask(i,j-1)=rmask(i,j)*rmask(i,j-1);
  end,
end,

%----------------------------------------------------------------------------
%  Compute R-factor.
%----------------------------------------------------------------------------

hx(1:L,1:Mp)=abs(h(2:Lp,1:Mp)-h(1:L,1:Mp))./(h(2:Lp,1:Mp)+h(1:L,1:Mp));
hy(1:Lp,1:M)=abs(h(1:Lp,2:Mp)-h(1:Lp,1:M))./(h(1:Lp,2:Mp)+h(1:Lp,1:M));

hx=hx.*umask;
hy=hy.*vmask;

r(1:L,1:M)=max(max(hx(1:L,1:M),hx(1:L,2:Mp)), ...
               max(hy(1:L,1:M),hy(2:Lp,1:M)));

rmin=min(min(r));
rmax=max(max(r));
ravg=mean(mean(r));
rmed=median(median(r));

disp('  ');
disp(['Minimum r-value = ', num2str(rmin)]);
disp(['Maximum r-value = ', num2str(rmax)]);
disp(['Mean    r-value = ', num2str(ravg)]);
disp(['Median  r-value = ', num2str(rmed)]);

return
