/*!\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->law=0;
	this->apply_dim = 1.;
	this->domaintype=-1;
	this->vx_input=NULL;
	this->vy_input=NULL;
	this->vz_input=NULL;
}
/*}}}*/
Friction::Friction(Element* element_in){/*{{{*/
	/* Determine the dimension according to the domain type automatically. 
	 * There are exceptions, e.g. HO, which needs the user to specify the dimension used in Friciton.*/
	this->element=element_in;

	/* Load necessary parameters */
	element_in->FindParam(&this->law,FrictionLawEnum);
	element_in->FindParam(&this->domaintype,DomainTypeEnum);

	/* Load VxBase and VyBase for this special case */
	switch(this->domaintype){
		case Domain2DhorizontalEnum: 
			this->apply_dim = 2.;
			this->vx_input = element_in->GetInput(VxBaseEnum);			_assert_(this->vx_input); 
			this->vy_input = element_in->GetInput(VyBaseEnum);			_assert_(this->vy_input);
			this->vz_input = NULL;
			break;
      case Domain2DverticalEnum:
			this->apply_dim = 2.;
			this->vx_input = element_in->GetInput(VxEnum);				_assert_(this->vx_input);
			this->vy_input = element_in->GetInput(VyEnum);				_assert_(this->vy_input);
			this->vz_input = NULL;
			break;
      case Domain3DEnum:           
			this->apply_dim = 3.;
			this->vx_input = element_in->GetInput(VxEnum);				_assert_(this->vx_input);
			this->vy_input = element_in->GetInput(VyEnum);				_assert_(this->vy_input);
			this->vz_input = element_in->GetInput(VzEnum);				_assert_(this->vz_input);
			break;
      default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
	}
}
/*}}}*/
Friction::Friction(Element* element_in,int dim) : Friction(element_in) {/*{{{*/
	this->apply_dim = reCast<IssmDouble>(dim);
}
/*}}}*/
Friction::Friction(Element* element_in,IssmDouble dim) : Friction(element_in) {/*{{{*/
	this->apply_dim = dim;
}
/*}}}*/
Friction::~Friction(){/*{{{*/
}
/*}}}*/


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

	switch(this->law){
		case 1:
			GetAlphaViscousComplement(palpha_complement,gauss);
			break;
		case 2:
			GetAlphaWeertmanComplement(palpha_complement, gauss);
			break;
		case 3:
			GetAlphaHydroComplement(palpha_complement,gauss);
			break;
		case 4:
			GetAlphaTempComplement(palpha_complement,gauss);
			break;
		case 11:
			GetAlphaSchoofComplement(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-1), 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-1)*/

	/*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);

	if(s==1.){
		/*This is to make AD happy and avoid 0^0*/
		alpha_complement=pow(Neff,r);
	}
	else{
		/*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::GetAlphaSchoofComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/

	/* Compute the complement of Schoof's law for inversion
	 * d alpha2                       
	 * -------- = |u_b|^(m-1) *(1+ (C/(Cmax N))^(1/m)|u_b|)^(-m-1)
	 *  dC                           
	*/
	/*diverse: */
	IssmDouble  m,Cmax;
	IssmDouble  C, coeff;
	IssmDouble  alpha_complement;

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

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

	/*Check to prevent dividing by zero if vmag==0*/
	if((vmag==0.) || (Neff == 0.)) {
		alpha_complement=0.;
	}
	else {
		alpha_complement= pow(vmag, m-1.)*pow((1 + pow(C/(Cmax*Neff),1./m)*vmag), -m-1.);
	}
	/*Assign output pointers:*/
	*palpha_complement=alpha_complement;
}/*}}}*/
void Friction::GetAlphaWeertmanComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/

	/* Compute the complement of Weertman's law for inversion
	 * alpha2 = C^2 * vel^(1/m-1)
	 * alpha_complement = vel^(1/m-1)
	*/
	/*diverse: */
	IssmDouble  m;
	IssmDouble  alpha_complement;

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

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

	/*Check to prevent dividing by zero if vmag==0*/
	if(vmag==0. && (1./m-1.)<0.) alpha_complement=0.;
	else alpha_complement= pow(vmag, 1.0/m-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;
		case 13:
			GetAlpha2Coulomb2(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);

	if(s==1.){
		/*This is to make AD happy and avoid 0^0*/
		alpha2=drag_coefficient*drag_coefficient*pow(Neff,r);
	}
	else{
		/*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 vmag = VelMag(gauss);

	bool ispwStochastic;
	IssmDouble Neff;
	element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum);
	if(ispwStochastic){
		/*Retrieve stochastic water pressure and compute ice pressure*/
		IssmDouble p_ice,p_water,Neff_limit;
		element->GetInputValue(&p_water,gauss,FrictionWaterPressureEnum);
		element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
		p_ice = IcePressure(gauss);
		Neff  = max(Neff_limit*p_ice, p_ice - p_water);
	}	
	else{
		/*Compute effective pressure directly*/
		Neff = EffectivePressure(gauss);
	}

	/*Check to prevent dividing by zero if vmag==0*/
	if(s==1.){
		/*This is to make AD happy and avoid 0^0*/
		alpha2=drag_coefficient*drag_coefficient*pow(Neff,r);
	}
	else{
		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);

	if(s==1.){
		/*This is to make AD happy and avoid 0^0*/
		alpha2=drag_coefficient*drag_coefficient*pow(Neff,r);
	}
	else{
		/*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^2 |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=C*C*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,coeff,Cmax,m,alpha2;

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

	/* scale C for a better inversion */
	C = coeff*coeff;

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

	/*Compute alpha^2*/
	if((ub<1e-10) ||(N==0.0)){
		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;
}/*}}}*/
void Friction::GetAlpha2Coulomb2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/

	/*This routine calculates the basal friction coefficient
	 *
	 *               C^2 |u_b|^(m-1) * (.5*N)
	 * alpha2= ___________________________________
	 *          (C^(2/m) |u_b| + (0.5*N)^(1/m) )^m
	 *
	 * */

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

	/*Recover parameters: */
	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= (pow(C,2)*pow(ub,m-1.)*(0.5*N)) / pow(pow(C,2./m)*ub + pow(0.5*N,1./m),m);
	}

	/*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);

	/*Compute ice pressure*/
	p_ice = IcePressure(gauss);

	/*From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing:*/
	switch(coupled_flag){
		case 0:{
			element->GetInputValue(&base, gauss,BaseEnum);
			element->GetInputValue(&sealevel, gauss,SealevelEnum);
			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
			p_water = rho_water*gravity*(sealevel-base);
			Neff = p_ice - p_water;
		}
			break;
		case 1:{
			p_water = 0.;
			Neff = p_ice - p_water;
		}
			break;
		case 2:{
			element->GetInputValue(&base, gauss,BaseEnum);
			element->GetInputValue(&sealevel, gauss,SealevelEnum);
			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
			p_water = max(0.,rho_water*gravity*(sealevel-base));
			Neff = p_ice - p_water;
		}
			break;
		case 3:{
			element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
		}
			break;
		case 4:{
			element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
		}
			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::IcePressure(Gauss* gauss){/*{{{*/
	/*Get ice pressure*/

	IssmDouble  thickness,p_ice;
	/*Recover Inputs and Parameters*/
	element->GetInputValue(&thickness, gauss,ThicknessEnum);
	IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum);
	IssmDouble gravity = element->FindParam(ConstantsGEnum);

	/*Compute*/
	p_ice = gravity*rho_ice*thickness;

	/*Return ice pressure*/
	return p_ice;

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

	int         coupled_flag;
	IssmDouble  base,sealevel,p_water;

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

	switch(coupled_flag){
		case 0:{
			element->GetInputValue(&base, gauss,BaseEnum);
			element->GetInputValue(&sealevel, gauss,SealevelEnum);
			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
			p_water = rho_water*gravity*(sealevel-base);
		}
			break;
		case 1:{
			p_water = 0.;
		}
			break;
		case 2:{
			element->GetInputValue(&base, gauss,BaseEnum);
			element->GetInputValue(&sealevel, gauss,SealevelEnum);
			IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum);
			IssmDouble gravity   = element->FindParam(ConstantsGEnum);
			p_water = max(0.,rho_water*gravity*(sealevel-base));
		}
			break;
		case 3:{
			_error_("water pressure not computed for coupling==3 in friction law");
		}
			break;
		case 4:{
			_error_("water pressure not computed for coupling==4 in friction law");
		}
			break;
		default:
			_error_("not supported");
	}

	/*Return water pressure*/
	return p_water;

}/*}}}*/
IssmDouble Friction::VelMag(Gauss* gauss){/*{{{*/
	/*Get the velocity magnitude as a function of flag */

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

	this->vx_input->GetInputValue(&vx, gauss);
	this->vy_input->GetInputValue(&vy, gauss);

	if ((this->vz_input == NULL) || (this->apply_dim<3.)) vz = 0.0;
	else this->vz_input->GetInputValue(&vz, gauss);

	if (this->apply_dim<2.) vy = 0.0;

	vmag = sqrt(vx*vx+vy*vy+vz*vz);
	return vmag;
}/*}}}*/
void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, Gauss* gauss){/*{{{*/

	this->vx_input->GetInputValue(pvx, gauss);
	/*Checks*/
	_assert_(!xIsNan<IssmDouble>(*pvx));
	_assert_(!xIsInf<IssmDouble>(*pvx));
}/*}}}*/
void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, Gauss* gauss){/*{{{*/

	this->vx_input->GetInputValue(pvx, gauss);
	this->vy_input->GetInputValue(pvy, gauss);
	/*Checks*/
	_assert_(!xIsNan<IssmDouble>(*pvx));
	_assert_(!xIsInf<IssmDouble>(*pvx));
	_assert_(!xIsNan<IssmDouble>(*pvy));
	_assert_(!xIsInf<IssmDouble>(*pvy));
}/*}}}*/
void Friction::GetBasalSlidingSpeeds(IssmDouble* pvx, IssmDouble* pvy, IssmDouble* pvz, Gauss* gauss){/*{{{*/

	this->vx_input->GetInputValue(pvx, gauss);
	this->vy_input->GetInputValue(pvy, gauss);
	this->vz_input->GetInputValue(pvz, gauss);
	/*Checks*/
	_assert_(!xIsNan<IssmDouble>(*pvx));
	_assert_(!xIsInf<IssmDouble>(*pvx));
	_assert_(!xIsNan<IssmDouble>(*pvy));
	_assert_(!xIsInf<IssmDouble>(*pvy));
	_assert_(!xIsNan<IssmDouble>(*pvz));
	_assert_(!xIsInf<IssmDouble>(*pvz));
}/*}}}*/
