function alpha2=Getalpha2(parameters)
%GETALPHA2 - compute the friction parameter alpha2
%
%   This routine calculates the basal friction coefficient
%   o alpha2= drag^2 * Neff ^r * vel ^s
%   with 
%   o Neff=rho_ice*g*thickness+rho_ice*g*bed, 
%   o r=q/p and s=1/p
%
%   Usage:
%      alpha2=Getalpha2(parameters)

%recover parameters 
gravity=parameters.g;
rho_ice=parameters.rho_ice;
rho_water=parameters.rho_water;
K=parameters.k;
bed=parameters.b;
thickness=parameters.h;
velocity_x=parameters.velocities(:,1);
velocity_y=parameters.velocities(:,2);
if size(parameters.velocities,2)==3, %we have to consider vz if Stokes model
	velocity_z=parameters.velocities(:,3);
end

%From bed and thickness, compute effective pressure when drag is viscous:
Neff=gravity*(rho_ice*thickness+rho_water*bed);

%If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because 
%the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour 
%for this should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should 
%replace it by Neff=0 (ie, equival it to an ice shelf)*/
pos=find(Neff<0);
Neff(pos)=0;

%recover p and q in the drag law u = k * sigma ^ p * Neff ^ q :
pcoeff=parameters.p;
qcoeff=parameters.q;
rcoeff=qcoeff./pcoeff;
scoeff=1./pcoeff;

%We need the velocity magnitude to evaluate the basal stress:
if size(parameters.velocities,2)==2,
	velocity_mag=sqrt(velocity_x.^2+velocity_y.^2);
elseif size(parameters.velocities,2)==3,
	velocity_mag=sqrt(velocity_x.^2+velocity_y.^2+velocity_z.^2);
end

alpha2=(K.^2).*(Neff.^rcoeff).*velocity_mag.^(scoeff-1);
