function nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff)
%VISCOSITY - compute the viscosity
%
%   This routine calculates the averaged viscosity over each 
%   element using Glen's flow law.
%
%   Usage:
%      nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff)


if (isempty(u) & isempty(v)),
	nu_bar=10^14*ones(nel,1);
else
	nu_bar=zeros(nel,1);

	power=(glen_coeff-1)./(2*glen_coeff);
	[ux ,uy ,vx, vy]=shear(index,alpha,beta,u,v);
	second_inv=(ux.^2 + vy.^2 + ((uy+vx).^2)/4 + ux.*vy);

	location=find(second_inv~=0);
	nu_bar(location)=B_bar(location)./(second_inv(location).^power(location));

	location=find(second_inv==0 & power~=0);
	nu_bar(location)=10^18; 	%arbitrary maximum viscosity to apply where there is no effective shear
							
	location=find(second_inv==0 & power==0);
	nu_bar(location)=B_bar(location);
end
