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



function bulk_mix(ml_index,m)

global g rb bulk
global nz z d
global u v

rvc = rb;
for j = ml_index+1:nz
	h 	= z(j);
	dd 	= (d(j)-d(1))/d(1);
	dv 	= (u(j)-u(1))^2+(v(j)-v(1))^2;
	if dv == 0
		rv = +Inf;
	else
		rv = g*h*dd/dv;
	end
	if rv > rvc
        bulk(m)=j; % save the index of the bulk mixing for diagnostic 
		break
    else
		mix5(j);
	end
end


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