Index: issm/trunk/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- issm/trunk/src/c/analyses/EnthalpyAnalysis.cpp	(revision 21729)
+++ issm/trunk/src/c/analyses/EnthalpyAnalysis.cpp	(revision 22758)
@@ -20,6 +20,8 @@
 	/*Output*/
 	IssmDouble *spcvector  = NULL;
+	IssmDouble *spcvectorstatic  = NULL;
 	IssmDouble* times=NULL;
 	IssmDouble* values=NULL;
+	IssmDouble* issurface = NULL;
 
 	/*Fetch parameters: */
@@ -32,13 +34,7 @@
 
 	/*Fetch data: */
+	iomodel->FetchData(&issurface,&M,&N,"md.mesh.vertexonsurface"); _assert_(N>0); _assert_(M==iomodel->numberofvertices);
 	iomodel->FetchData(&spcvector,&M,&N,"md.thermal.spctemperature");
-
-	/*Convert spcs from temperatures to enthalpy*/
-	_assert_(N>0); _assert_(M>=iomodel->numberofvertices);
-	for(int i=0;i<iomodel->numberofvertices;i++){
-		for(int j=0;j<N;j++){
-			spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
-		}
-	}
+	iomodel->FetchData(&spcvectorstatic,&M,&N,"md.thermal.spctemperature");
 
 	/*Specific case for PDD, we want the constaints to be updated by the PDD scheme itself*/
@@ -51,6 +47,27 @@
 	}
 
+	/*Convert spcs from temperatures to enthalpy*/
+	_assert_(N>0); _assert_(M>=iomodel->numberofvertices);
+	for(int i=0;i<iomodel->numberofvertices;i++){
+		for(int j=0;j<N;j++){
+			if (isdynamic){
+				if (issurface[i]==1){
+					spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
+					spcvectorstatic[i*N+j] = NAN;
+				}
+				else{
+					spcvector[i*N+j] = NAN;
+					spcvectorstatic[i*N+j] = heatcapacity*(spcvectorstatic[i*N+j]-referencetemperature);
+				}
+			}
+			else{
+				spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
+			}
+		}
+	}
+
 	if(isdynamic){
-		IoModelToDynamicConstraintsx(constraints,iomodel,spcvector,M,N,EnthalpyAnalysisEnum,finiteelement);
+		IoModelToDynamicConstraintsx(constraints,iomodel,spcvector,iomodel->numberofvertices,1,EnthalpyAnalysisEnum,finiteelement);
+		IoModelToConstraintsx(constraints,iomodel,spcvectorstatic,M,N,EnthalpyAnalysisEnum,finiteelement);
 	}
 	else{
@@ -60,4 +77,6 @@
 	/*Free ressources:*/
 	iomodel->DeleteData(spcvector,"md.thermal.spctemperature");
+	iomodel->DeleteData(spcvectorstatic,"md.thermal.spctemperature");
+	iomodel->DeleteData(issurface,"md.mesh.vertexonsurface");
 	xDelete<IssmDouble>(times);
 	xDelete<IssmDouble>(values);
@@ -174,7 +193,11 @@
 	switch(frictionlaw){
 		case 1:
+			iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
 			iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
 			iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
 			iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
+			if (FrictionCoupling==1){
+				iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
+			}
 			break;
 		case 2:
@@ -187,5 +210,5 @@
 			iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum);
 			iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
-			if (FrictionCoupling==0){
+			if (FrictionCoupling==1){
 				iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
 			}
@@ -210,4 +233,14 @@
 			iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum);
 			break;
+		case 7:
+			iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
+			iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
+			iomodel->FetchDataToInput(elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum);
+			iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
+			iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
+			if (FrictionCoupling==1){
+				iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
+			}
+			break;
 		case 9:
 			iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
@@ -217,5 +250,5 @@
 			break;
 		default:
-			_error_("not supported");
+			_error_("friction law not supported");
 	}
 
@@ -244,7 +277,14 @@
 	int frictionlaw;
 	iomodel->FindConstant(&frictionlaw,"md.friction.law");
-	if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
-	if(frictionlaw==3) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
-	if(frictionlaw==9) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+	if(frictionlaw==4 || frictionlaw==6){
+		parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+	}
+	if(frictionlaw==3 || frictionlaw==1 || frictionlaw==7){
+		parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
+	}
+	if(frictionlaw==9){
+		parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+		parameters->AddObject(new IntParam(FrictionCouplingEnum,0));
+	}
 }/*}}}*/
 
@@ -302,4 +342,8 @@
 		ComputeBasalMeltingrate(element);
 	}
+
+	/*extrude inputs*/
+	femmodel->parameters->SetParam(BasalforcingsGroundediceMeltingRateEnum,InputToExtrudeEnum);
+	extrudefrombase_core(femmodel);
 }/*}}}*/
 void           EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/
@@ -525,7 +569,5 @@
 	IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
 	IssmDouble*    dbasis   = xNew<IssmDouble>(3*numnodes);
-	IssmDouble*    B        = xNew<IssmDouble>(3*numnodes);
 	IssmDouble*    Bprime   = xNew<IssmDouble>(3*numnodes);
-	IssmDouble     D[3][3]  = {0.};
 	IssmDouble     K[3][3];
 
@@ -556,39 +598,39 @@
 
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
+		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
 		D_scalar=gauss->weight*Jdet;
 		if(dt!=0.) D_scalar=D_scalar*dt;
 
 		/*Conduction: */
-		GetBConduct(B,element,xyz_list,gauss); 
-		D[0][0]=D_scalar*kappa/rho_ice;
-		D[1][1]=D_scalar*kappa/rho_ice;
-		D[2][2]=D_scalar*kappa/rho_ice;
-		TripleMultiply(B,3,numnodes,1,
-					&D[0][0],3,3,0,
-					B,3,numnodes,0,
-					&Ke->values[0],1);
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += D_scalar*kappa/rho_ice*(
+							dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
+							);
+			}
+		}
 
 		/*Advection: */
-		GetBAdvec(B,element,xyz_list,gauss); 
-		GetBAdvecprime(Bprime,element,xyz_list,gauss); 
 		vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
 		vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
 		vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
-		D[0][0]=D_scalar*vx;
-		D[1][1]=D_scalar*vy;
-		D[2][2]=D_scalar*vz;
-		TripleMultiply(B,3,numnodes,1,
-					&D[0][0],3,3,0,
-					Bprime,3,numnodes,0,
-					&Ke->values[0],1);
+		for(int i=0;i<numnodes;i++){
+			for(int j=0;j<numnodes;j++){
+				Ke->values[i*numnodes+j] += D_scalar*(
+							vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i] +vz*dbasis[2*numnodes+j]*basis[i]
+							);
+			}
+		}
 
 		/*Transient: */
 		if(dt!=0.){
 			D_scalar=gauss->weight*Jdet;
-			element->NodalFunctions(basis,gauss);
-			TripleMultiply(basis,numnodes,1,0,
-						&D_scalar,1,1,0,
-						basis,1,numnodes,0,
-						&Ke->values[0],1);
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
+				}
+			}
 			D_scalar=D_scalar*dt;
 		}
@@ -634,5 +676,4 @@
 	xDelete<IssmDouble>(basis);
 	xDelete<IssmDouble>(dbasis);
-	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(Bprime);
 	delete gauss;
@@ -1063,6 +1104,6 @@
 		xDelete<IssmDouble>(watercolumn);
 		xDelete<IssmDouble>(drainage_int);
-	}
-	xDelete<int>(basalnodeindices);
+		xDelete<int>(basalnodeindices);
+	}
 }/*}}}*/
 void				EnthalpyAnalysis::DrainageUpdateEnthalpy(FemModel* femmodel){/*{{{*/
@@ -1510,5 +1551,5 @@
 		/*Update Rheology only if converged (we must make sure that the temperature is below melting point
 		 * otherwise the rheology could be negative*/
-		element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+		rheology_law=element->GetIntegerMaterialParameter(MaterialsRheologyLawEnum);
 		element->GetInputListOnNodes(&surface[0],SurfaceEnum);
 		switch(rheology_law){
