/*!\file Friction.c
 * \brief: implementation of the Friction object
 */

/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "../classes.h"
#include "shared/shared.h"
/*}}}*/

/*Constructors/destructors*/
Friction::Friction(){/*{{{*/
	this->element=NULL;
	this->dim=0;
	this->law=0;

}
/*}}}*/
Friction::Friction(Element* element_in,int dim_in){/*{{{*/

	this->element=element_in;
	this->dim=dim_in;
	element_in->FindParam(&this->law,FrictionLawEnum);
}
/*}}}*/
Friction::~Friction(){/*{{{*/
}
/*}}}*/

/*methods: */
void Friction::Echo(void){/*{{{*/
	_printf_("Friction:\n");
	_printf_("   dim: " << this->dim<< "\n");
}
/*}}}*/
void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/

	switch(this->law){
		case 1:
			GetAlphaViscousComplement(palpha_complement,gauss);
			break;
		case 3:
			GetAlphaHydroComplement(palpha_complement,gauss);
			break;
		case 4:
			GetAlphaTempComplement(palpha_complement,gauss);
			break;
	  default:
			_error_("not supported");
	}


	/*Checks*/
	_assert_(!xIsNan<IssmDouble>(*palpha_complement));
	_assert_(!xIsInf<IssmDouble>(*palpha_complement));

}/*}}}*/
void Friction::GetAlphaHydroComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/

	/*diverse: */
	IssmDouble  q_exp;
	IssmDouble  C_param;
	IssmDouble  As;
	IssmDouble  n;
	IssmDouble  alpha;
	IssmDouble  Chi,Gamma;
	IssmDouble  alpha_complement;

	/*Recover parameters: */
	element->GetInputValue(&q_exp,gauss,FrictionQEnum);
	element->GetInputValue(&C_param,gauss,FrictionCEnum);
	element->GetInputValue(&As,gauss,FrictionAsEnum);
	element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);

	/*Get effective pressure and velocity magnitude*/
	IssmDouble Neff = EffectivePressure(gauss);
	IssmDouble vmag = VelMag(gauss);

	if (q_exp==1){
		alpha=1;
	}
	else{
		alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
	}
	Chi   = vmag/(pow(C_param,n)*pow(Neff,n)*As);
	Gamma = (Chi/(1.+alpha*pow(Chi,q_exp)));
	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0.) alpha_complement=0.;
	else	if(Neff==0.) alpha_complement=0.;
	else	alpha_complement=-(C_param*Neff/(n*vmag)) *
					pow(Gamma,((1.-n)/n)) *
					(Gamma/As - (alpha*q_exp*pow(Chi,q_exp-1.)* Gamma * Gamma/As));

	/*Assign output pointers:*/
	*palpha_complement=alpha_complement;
}
/*}}}*/
void Friction::GetAlphaTempComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
	/*Here, we want to parameterize the friction as a function of temperature
	 *
	 * alpha2 = alpha2_viscous * 1/f(T)
	 *
	 * where f(T) = exp((T-Tpmp)/gamma)
	 */

	/*Intermediaries: */
	IssmDouble  f,T,pressure,Tpmp,gamma;
	IssmDouble  alpha_complement;

	/*Get viscous part*/
	this->GetAlphaViscousComplement(&alpha_complement,gauss);

	/*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
	element->GetInputValue(&T,gauss,TemperatureEnum);
	element->GetInputValue(&pressure,gauss,PressureEnum);
	Tpmp = element->TMeltingPoint(pressure);

	/*Compute scaling parameter*/
	element->parameters->FindParam(&gamma,FrictionGammaEnum);
	alpha_complement = alpha_complement/ (exp((T-Tpmp)/gamma)+1e-3);

	/*Assign output pointers:*/
	*palpha_complement=alpha_complement;
}/*}}}*/
void Friction::GetAlphaViscousComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/

	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p.
	 * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
	 * alpha_complement= Neff ^r * vel ^s*/

	/*diverse: */
	IssmDouble  r,s;
	IssmDouble  drag_p,drag_q;
	IssmDouble  drag_coefficient;
	IssmDouble  alpha_complement;

	/*Recover parameters: */
	element->GetInputValue(&drag_p,gauss,FrictionPEnum);
	element->GetInputValue(&drag_q,gauss,FrictionQEnum);
	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);

	/*compute r and q coefficients: */
	r=drag_q/drag_p;
	s=1./drag_p;

	/*Get effective pressure*/
	IssmDouble Neff = EffectivePressure(gauss);
	IssmDouble vmag = VelMag(gauss);

	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
	else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));

	/*Assign output pointers:*/
	*palpha_complement=alpha_complement;
}
/*}}}*/
void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	switch(this->law){
		case 1:
			GetAlpha2Viscous(palpha2,gauss);
			break;
		case 2:
			GetAlpha2Weertman(palpha2,gauss);
			break;
		case 3:
			GetAlpha2Hydro(palpha2,gauss);
			break;
		case 4:
			GetAlpha2Temp(palpha2,gauss);
			break;
		case 5:
			GetAlpha2WaterLayer(palpha2,gauss);
			break;
		case 6:
			GetAlpha2WeertmanTemp(palpha2,gauss);
			break;
		case 7:
			GetAlpha2Coulomb(palpha2,gauss);
			break;
		case 8:
			GetAlpha2Shakti(palpha2,gauss);
			break;
		case 9:
			GetAlpha2Josh(palpha2,gauss);
			break;
		case 10:
			GetAlpha2PISM(palpha2,gauss);
			break;
		case 11:
			GetAlpha2Schoof(palpha2,gauss);
			break;
		case 12:
			GetAlpha2Tsai(palpha2,gauss);
			break;
	  default:
			_error_("Friction law "<< this->law <<" not supported");
	}

	/*Checks*/
	_assert_(!xIsNan<IssmDouble>(*palpha2));
	_assert_(!xIsInf<IssmDouble>(*palpha2));
	_assert_(*palpha2>=0);

}/*}}}*/
void Friction::GetAlpha2Coulomb(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/*This routine calculates the basal friction coefficient
	  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/

	/*diverse: */
	IssmDouble  r,s;
	IssmDouble  drag_p, drag_q;
	IssmDouble  thickness,base,floatation_thickness,sealevel;
	IssmDouble  drag_coefficient,drag_coefficient_coulomb;
	IssmDouble  alpha2,alpha2_coulomb;

	/*Recover parameters: */
	element->GetInputValue(&drag_p,gauss,FrictionPEnum);
	element->GetInputValue(&drag_q,gauss,FrictionQEnum);
	element->GetInputValue(&thickness, gauss,ThicknessEnum);
	element->GetInputValue(&base, gauss,BaseEnum);
	element->GetInputValue(&sealevel, gauss,SealevelEnum);
	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
	element->GetInputValue(&drag_coefficient_coulomb, gauss,FrictionCoefficientcoulombEnum);
	IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
	IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
	IssmDouble gravity   = element->FindParam(ConstantsGEnum);

	//compute r and q coefficients: */
	r=drag_q/drag_p;
	s=1./drag_p;

	/*Get effective pressure*/
	IssmDouble Neff = EffectivePressure(gauss);
	IssmDouble vmag = VelMag(gauss);

	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0. && (s-1.)<0.){
		alpha2=0.;
	}
	else{
		alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
	}

	floatation_thickness=0;
	if(base<0) floatation_thickness=-(rho_water/rho_ice)*base;
	if(vmag==0.){
		alpha2_coulomb=0.;
	}
	else{
		alpha2_coulomb=drag_coefficient_coulomb*drag_coefficient_coulomb*rho_ice*gravity*max(0.,thickness-floatation_thickness)/vmag;
	}

	if(alpha2_coulomb<alpha2) alpha2=alpha2_coulomb;

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlpha2Hydro(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/*This routine calculates the basal friction coefficient
		Based on Gagliardini 2007, needs a good effective pressure computation
		Not tested so far so use at your own risks
	  alpha2= NeffC[Chi/(1+alpha*Chi^q)]^(1/n)*|vel|^(-1)  with
		 Chi=|vel|/(C^n*Neff^n*As)
		 alpha=(q-1)^(q-1)/q^q */

	/*diverse: */
	IssmDouble  q_exp;
	IssmDouble  C_param;
	IssmDouble  As;
	IssmDouble  n;
	IssmDouble  alpha;
	IssmDouble  Chi,Gamma;
	IssmDouble  alpha2;

	/*Recover parameters: */
	element->GetInputValue(&q_exp,gauss,FrictionQEnum);
	element->GetInputValue(&C_param,gauss,FrictionCEnum);
	element->GetInputValue(&As,gauss,FrictionAsEnum);
	element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);

	/*Get effective pressure*/
	IssmDouble Neff = EffectivePressure(gauss);
	IssmDouble vmag = VelMag(gauss);

	//compute alpha and Chi coefficients: */
	if (q_exp==1){
		alpha=1;
	}
	else{
		alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
	}
	Chi=vmag/(pow(C_param,n)*pow(Neff,n)*As);
	Gamma=(Chi/(1. + alpha * pow(Chi,q_exp)));
	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0.) alpha2=0.;
	else	if (Neff==0) alpha2=0.0;
	else	alpha2=Neff * C_param * pow(Gamma,1./n) * pow(vmag,-1);

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlpha2Shakti(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/* FrictionGetAlpha2 computes alpha2= drag^2 * Neff, with Neff=rho_ice*g*thickness+rho_ice*g*(head-base)*/

	/*diverse: */
	IssmDouble  pressure_ice,pressure_water;
	IssmDouble  Neff;
	IssmDouble  drag_coefficient;
	IssmDouble  base,thickness,head,sealevel;
	IssmDouble  alpha2;

	/*Recover parameters: */
	element->GetInputValue(&thickness, gauss,ThicknessEnum);
	element->GetInputValue(&base, gauss,BaseEnum);
	element->GetInputValue(&head, gauss,HydrologyHeadEnum);
	element->GetInputValue(&sealevel, gauss,SealevelEnum);
	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
	IssmDouble rho_water   = element->FindParam(MaterialsRhoFreshwaterEnum);
	IssmDouble rho_ice     = element->FindParam(MaterialsRhoIceEnum);
	IssmDouble gravity     = element->FindParam(ConstantsGEnum);

	//From base and thickness, compute effective pressure when drag is viscous:
	pressure_ice   = rho_ice*gravity*thickness;
	pressure_water = rho_water*gravity*(head-base+sealevel);
	Neff=pressure_ice-pressure_water;
	if(Neff<0.) Neff=0.;

	alpha2=drag_coefficient*drag_coefficient*Neff;

	/*Assign output pointers:*/
	*palpha2=alpha2;
}
/*}}}*/
void Friction::GetAlpha2Temp(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
	/*Here, we want to parameterize the friction as a function of temperature
	 *
	 * alpha2 = alpha2_viscous * 1/f(T)
	 *
	 * where f(T) = exp((T-Tpmp)/gamma)
	 */

	/*Intermediaries: */
	IssmDouble  f,T,pressure,Tpmp,gamma;
	IssmDouble  alpha2;

	/*Get viscous part*/
	this->GetAlpha2Viscous(&alpha2,gauss);

	/*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
	element->GetInputValue(&T,gauss,TemperatureEnum);
	element->GetInputValue(&pressure,gauss,PressureEnum);
	Tpmp = element->TMeltingPoint(pressure);

	/*Compute scaling parameter*/
	element->parameters->FindParam(&gamma,FrictionGammaEnum);
	alpha2 = alpha2 / (exp((T-Tpmp)/gamma)+1e-3);

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlpha2Josh(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
	/*Here, we want to parameterize the friction as a function of temperature
	 *
	 * alpha2 = alpha2_viscous * 1/f(T)
	 *
	 * where f(T) = exp((T-Tpmp)/gamma)
	 */

	/*Intermediaries: */
	IssmDouble  T,Tpmp,deltaT,deltaTref,pressure,diff,drag_coefficient;
	IssmDouble  alpha2,time,gamma,ref,alp_new,alphascaled;
	const IssmDouble yts = 365*24*3600.;

	/*Get viscous part*/
	this->GetAlpha2Viscous(&alpha2,gauss);

	/*Get delta Refs*/
	element->GetInputValue(&deltaTref,gauss,FrictionPressureAdjustedTemperatureEnum);
	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
	/*New*/
	/*element->GetInputValue(&deltaTrefsfc,gauss,FrictionSurfaceTemperatureEnum);
	 *    element->GetInputValue(&Tpdd,gauss,TemperaturePDDEnum);
	 *       */

	/*Compute delta T*/
	element->GetInputValue(&T,gauss,TemperatureEnum);
	element->GetInputValue(&pressure,gauss,PressureEnum);
	Tpmp = element->TMeltingPoint(pressure);
	deltaT = T-Tpmp;


	/*Compute gamma*/
	element->parameters->FindParam(&time,TimeEnum);
	element->parameters->FindParam(&gamma,FrictionGammaEnum);

	ref = exp(deltaTref/gamma);
	alp_new = ref/exp(deltaT/gamma);

	alphascaled = sqrt(alp_new)*drag_coefficient;
	if (alphascaled > 300) alp_new = (300/drag_coefficient)*(300/drag_coefficient);

	alp_new=alp_new*alpha2;

	/*Assign output pointers:*/
	*palpha2=alp_new;
}/*}}}*/
void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/*This routine calculates the basal friction coefficient
	  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/

	/*diverse: */
	IssmDouble  r,s;
	IssmDouble  drag_p, drag_q;
	IssmDouble  drag_coefficient;
	IssmDouble  alpha2;

	/*Recover parameters: */
	element->GetInputValue(&drag_p,gauss,FrictionPEnum);
	element->GetInputValue(&drag_q,gauss,FrictionQEnum);
	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);

	/*compute r and q coefficients: */
	r=drag_q/drag_p;
	s=1./drag_p;

	/*Get effective pressure and basal velocity*/
	IssmDouble Neff = EffectivePressure(gauss);
	IssmDouble vmag = VelMag(gauss);

	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0. && (s-1.)<0.) alpha2=0.;
	else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlpha2WaterLayer(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/*This routine calculates the basal friction coefficient
	  alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/

	/*diverse: */
	IssmDouble  r,s;
	IssmDouble  drag_p, drag_q;
	IssmDouble  Neff,F;
	IssmDouble  thickness,base,sealevel;
	IssmDouble  drag_coefficient,water_layer;
	IssmDouble  alpha2;

	/*Recover parameters: */
	element->parameters->FindParam(&F,FrictionFEnum);
	element->GetInputValue(&drag_p,gauss,FrictionPEnum);
	element->GetInputValue(&drag_q,gauss,FrictionQEnum);
	element->GetInputValue(&thickness, gauss,ThicknessEnum);
	element->GetInputValue(&base, gauss,BaseEnum);
	element->GetInputValue(&sealevel, gauss,SealevelEnum);
	element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
	element->GetInputValue(&water_layer, gauss,FrictionWaterLayerEnum);
	IssmDouble rho_water   = element->FindParam(MaterialsRhoSeawaterEnum);
	IssmDouble rho_ice     = element->FindParam(MaterialsRhoIceEnum);
	IssmDouble gravity     = element->FindParam(ConstantsGEnum);

	//compute r and q coefficients: */
	r=drag_q/drag_p;
	s=1./drag_p;

	//From base and thickness, compute effective pressure when drag is viscous:
	if(base>0) base=0;
	if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*(base-sealevel);
	else if(water_layer>0) Neff=gravity*rho_ice*thickness*F;
	else _error_("negative water layer thickness");
	if(Neff<0) Neff=0;

	IssmDouble vmag = VelMag(gauss);

	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0. && (s-1.)<0.) alpha2=0.;
	else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */

	/*diverse: */
	IssmDouble  C,m;
	IssmDouble  alpha2;

	/*Recover parameters: */
	element->GetInputValue(&C,gauss,FrictionCEnum);
	element->GetInputValue(&m,gauss,FrictionMEnum);

	/*Get velocity magnitude*/
	IssmDouble vmag = VelMag(gauss);

	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0. && (1./m-1.)<0.) alpha2=0.;
	else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.));

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlpha2WeertmanTemp(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
	/*Here, we want to parameterize the friction as a function of temperature
	 *
	 * alpha2 = alpha2_weertman * 1/f(T)
	 *
	 * where f(T) = exp((T-Tpmp)/gamma)
	 */

	/*Intermediaries: */
	IssmDouble  f,T,pressure,Tpmp,gamma;
	IssmDouble  alpha2;

	/*Get viscous part*/
	this->GetAlpha2Weertman(&alpha2,gauss);

	/*Get pressure melting point (Tpmp) for local pressure and get current temperature*/
	element->GetInputValue(&T,gauss,TemperatureEnum);
	element->GetInputValue(&pressure,gauss,PressureEnum);
	Tpmp = element->TMeltingPoint(pressure);

	/*Compute scaling parameter*/
	element->parameters->FindParam(&gamma,FrictionGammaEnum);
	alpha2 = alpha2 / exp((T-Tpmp)/gamma);

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlpha2PISM(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
	/*Here, we want to parameterize the friction using a pseudoplastic friction law,
	 * computing the basal shear stress as
	 *
	 * alpha2 = tau_c (u_b/(abs(u_b)^(1-q)*u_0^q))
	 *
	 * The yield stress tau_c is a function of the effective pressure N
	 * using a Mohr-Coloumb criterion, so that
	 * tau_c = tan(phi)*N,
	 * where phi is the till friction angle, representing sediment strength
	 *
	 * The effective pressure is given by Eq. (5) in Aschwanden et al. 2016:
	 *
	 * N = delta * P0 * 10^((e_0/Cc)(1-(W/Wmax)))
	 *
	 * W is calculated by a non-conserving hydrology model in HydrologyPismAnalysis.cpp
	 *
	 * see Aschwanden et al. 2016 and Bueler and Brown, 2009 for more details
	 */

	/*compute ice overburden pressure P0*/
	IssmDouble thickness,base,P0;
	element->GetInputValue(&thickness, gauss,ThicknessEnum);
	element->GetInputValue(&base, gauss,BaseEnum);
	//element->GetInputValue(&sealevel, gauss,SealevelEnum);
	IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
	IssmDouble gravity   = element->FindParam(ConstantsGEnum);
	P0 = gravity*rho_ice*thickness;

	/*Compute effective pressure*/
	IssmDouble  N,delta,W,Wmax,e0,Cc;
	element->parameters->FindParam(&delta,FrictionDeltaEnum);
	element->parameters->FindParam(&e0,FrictionVoidRatioEnum);
	element->GetInputValue(&Cc,gauss,FrictionSedimentCompressibilityCoefficientEnum);
	element->GetInputValue(&W,gauss,WatercolumnEnum);
	element->GetInputValue(&Wmax,gauss,HydrologyWatercolumnMaxEnum);

 	/*Check that water column height is within 0 and upper bound, correct if needed*/
 	if(W>Wmax) W=Wmax;
 	if(W<0)    W=0.;

	N = delta*P0*pow(10.,(e0/Cc)*(1.-W/Wmax));

	/*Get till friction angles, defined by user [deg]*/
	IssmDouble phi;
	element->GetInputValue(&phi,gauss,FrictionTillFrictionAngleEnum);

	/*Convert till friction angle from user-defined deg to rad, which Matlab uses*/
	phi = phi*PI/180.;

	/*Compute yield stress following a Mohr-Colomb criterion*/
	IssmDouble tau_c = N*tan(phi);

	/*Compute basal speed*/
	IssmDouble ub;
	element->GetInputValue(&ub,gauss,VelEnum);

	/*now compute alpha^2*/
	IssmDouble u0,q;
	element->parameters->FindParam(&u0,FrictionThresholdSpeedEnum);
	element->parameters->FindParam(&q,FrictionPseudoplasticityExponentEnum);
	IssmDouble alpha2 = tau_c/(pow(ub+1.e-10,1.-q)*pow(u0,q));

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlpha2Schoof(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/*This routine calculates the basal friction coefficient
	 *
	 *               C |u_b|^(m-1)
	 * alpha2= __________________________
	 *          (1+(C/(Cmax N))^1/m |u_b| )^m
	 *
	 * */

	/*diverse: */
	IssmDouble  C,Cmax,m,alpha2;

	/*Recover parameters: */
	element->GetInputValue(&Cmax,gauss,FrictionCmaxEnum);
	element->GetInputValue(&C,gauss,FrictionCEnum);
	element->GetInputValue(&m,gauss,FrictionMEnum);

	/*Get effective pressure and velocity magnitude*/
	IssmDouble N  = EffectivePressure(gauss);
	IssmDouble ub = VelMag(gauss);

	/*Compute alpha^2*/
	if(ub<1e-10){
		alpha2 = 0.;
	}
	else{
		alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*N),1./m)*ub,m);
	}

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/
void Friction::GetAlpha2Tsai(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/*This routine calculates the basal friction coefficient
	 *
	 * alpha2= min(C |ub|^m , f N ) / |ub|
	 *
	 * */

	/*diverse: */
	IssmDouble  C,f,m,alpha2;

	/*Recover parameters: */
	element->GetInputValue(&f,gauss,FrictionfEnum);
	element->GetInputValue(&C,gauss,FrictionCEnum);
	element->GetInputValue(&m,gauss,FrictionMEnum);

	/*Get effective pressure and velocity magnitude*/
	IssmDouble N  = EffectivePressure(gauss);
	IssmDouble ub = VelMag(gauss);

	/*Compute alpha^2*/
	if(ub<1e-10){
		alpha2 = 0.;
	}
	else{
		alpha2= C*pow(ub,m);

		if(alpha2>f*N) alpha2 = f*N;

		alpha2 = alpha2/ub;
	}

	/*Assign output pointers:*/
	*palpha2=alpha2;
}/*}}}*/

IssmDouble Friction::EffectivePressure(Gauss* gauss){/*{{{*/
	/*Get effective pressure as a function of  flag */

	/*diverse: */
	int         coupled_flag;
	IssmDouble  thickness,base,sealevel;
	IssmDouble  p_ice,p_water;
	IssmDouble  Neff,Neff_limit;

	/*Recover parameters: */
	element->parameters->FindParam(&coupled_flag,FrictionCouplingEnum);
	element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);

	/*From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing:*/
	switch(coupled_flag){
		case 0:{
			element->GetInputValue(&thickness, gauss,ThicknessEnum);
			element->GetInputValue(&base, gauss,BaseEnum);
			element->GetInputValue(&sealevel, gauss,SealevelEnum);
			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
			p_ice   = gravity*rho_ice*thickness;
			p_water = rho_water*gravity*(sealevel-base);
			Neff = p_ice - p_water;
		}
			break;
		case 1:{
			element->GetInputValue(&thickness, gauss,ThicknessEnum);
			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
			p_ice   = gravity*rho_ice*thickness;
			p_water = 0.;
			Neff = p_ice - p_water;
		}
			break;
		case 2:{
			element->GetInputValue(&thickness, gauss,ThicknessEnum);
			element->GetInputValue(&base, gauss,BaseEnum);
			element->GetInputValue(&sealevel, gauss,SealevelEnum);
			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
			p_ice   = gravity*rho_ice*thickness;
			p_water = max(0.,rho_water*gravity*(sealevel-base));
			Neff = p_ice - p_water;
		}
			break;
		case 3:{
			element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
			element->GetInputValue(&thickness, gauss,ThicknessEnum);
			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
			p_ice   = gravity*rho_ice*thickness;
		}
			break;
		case 4:{
			element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
			element->GetInputValue(&thickness, gauss,ThicknessEnum);
			IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
			p_ice   = gravity*rho_ice*thickness;
		}
			break;
		default:
			_error_("not supported");
	}

	/*Make sure Neff is positive*/
	if(Neff<Neff_limit*p_ice) Neff=Neff_limit*p_ice;

	/*Return effective pressure*/
	return Neff;

}/*}}}*/
IssmDouble Friction::VelMag(Gauss* gauss){/*{{{*/
	/*Get effective pressure as a function of  flag */


	/*diverse*/
	IssmDouble vx,vy,vz,vmag;

	switch(dim){
		case 1:
			element->GetInputValue(&vx,gauss,VxEnum);
			vmag=sqrt(vx*vx);
			break;
		case 2:
			element->GetInputValue(&vx,gauss,VxEnum);
			element->GetInputValue(&vy,gauss,VyEnum);
			vmag=sqrt(vx*vx+vy*vy);
			break;
		case 3:
			element->GetInputValue(&vx,gauss,VxEnum);
			element->GetInputValue(&vy,gauss,VyEnum);
			element->GetInputValue(&vz,gauss,VzEnum);
			vmag=sqrt(vx*vx+vy*vy+vz*vz);
			break;
		default:
			_error_("not supported");
	}

	return vmag;

}/*}}}*/
