function grad_mix(m)

%  This function performs the gradient Richardson Number relaxation
%  by mixing adjacent cells just enough to bring them to a new
%  Richardson Number.
%
% Updated to remove for loop and make it more efficient (FS-M, 4/6/09)
%

global dz g rg richard_grad
global nz z T s d
global u v
global kcount;   % iteration counter

rc 	= rg;

%  Compute the gradient Richardson Number, taking care to avoid dividing by
%  zero in the mixed layer.  The numerical values of the minimum allowable
%  density and velocity differences are entirely arbitrary, and should not
%  effect the calculations (except that on some occasions they evidently have!)

j1 = 1;
j2 = nz-1;

kcount = 0;

while 1
%  	for j = j1:j2
%  		if j <= 0
%  			keyboard
%  		end
%  		dd = (d(j+1)-d(j))/d(j);
%  		dv = (u(j+1)-u(j))^2+(v(j+1)-v(j))^2;
%  		if dv == 0
%  			r(j) = Inf;
%  		else
%  			r(j) = g*dz*dd/dv;
%  		end
%  	end

        dd = []; dv = [];
        dd(j1:j2) = (d(j1+1:j2+1)-d(j1:j2))./d(j1:j2);
        dv(j1:j2) = (u(j1+1:j2+1)-u(j1:j2)).^2 + (v(j1+1:j2+1)-v(j1:j2)).^2;

%        l = find(dv ~= 0);
%        r(j1:j2) = Inf;
        l = find(dv > 1e-6);
        r(j1:j2) = 9999;
        r(l) = g.*dz.*dd(l)./dv(l);
        
	%  Find the smallest value of r in profile
% 	rs = min(r);
% 	js = min(find(r==rs));
	[rs,js] = min(r);
    
    % keep the index of the smallest value of r for diagnostic
   
	%  Check to see whether the smallest r is critical or not.

	if rs > rc
        richard_grad(m)=js; % save the index of the minimum value of r for diagnostic
		break
	end

	%  Mix the cells js and js+1 that had the smallest Richardson Number
	stir(rc,rs,js);


        kcount = kcount + 1;
	%  Recompute the Richardson Number over the part of the profile that has changed

	j1 = js-2;
	if j1 < 1
		 j1 = 1;
	end
	j2 = js+2;
	if j2 > nz-1
		 j2 = nz-1;
	end
end


% =-------------------------------------------------------------------------=
