Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 16887)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 16888)
@@ -185,5 +185,196 @@
 /*Finite Element Analysis*/
 ElementMatrix* EnthalpyAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixVolume(element);
+	ElementMatrix* Ke2=CreateKMatrixShelf(element);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+
+	/*clean-up and return*/
+	delete Ke1;
+	delete Ke2;
+	return Ke;
+}/*}}}*/
+ElementMatrix* EnthalpyAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	int         stabilization;
+	IssmDouble  Jdet,dt,u,v,w,um,vm,wm,vel;
+	IssmDouble  h,hx,hy,hz,vx,vy,vz;
+	IssmDouble  tau_parameter,diameter;
+	IssmDouble  D_scalar;
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes    = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and other vectors*/
+	ElementMatrix* Ke       = element->NewElementMatrix();
+	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];
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	element->FindParam(&stabilization,ThermalStabilizationEnum);
+	IssmDouble  rho_water           = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+	IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble  gravity             = element->GetMaterialParameter(ConstantsGEnum);
+	IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
+	IssmDouble  thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
+	Input* vx_input  = element->GetInput(VxEnum);     _assert_(vx_input);
+	Input* vy_input  = element->GetInput(VyEnum);     _assert_(vy_input);
+	Input* vz_input  = element->GetInput(VzEnum);     _assert_(vz_input);
+	Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input);
+	Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
+	Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input);
+	if(stabilization==2) diameter=element->MinEdgeLength(xyz_list);
+
+	/*Enthalpy diffusion parameter*/
+	IssmDouble kappa=this->EnthalpyDiffusionParameterVolume(element); _assert_(kappa>0.);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,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);
+
+		/*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);
+
+		/*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);
+			D_scalar=D_scalar*dt;
+		}
+
+		/*Artifficial diffusivity*/
+		if(stabilization==1){
+			element->ElementSizes(&hx,&hy,&hz);
+			vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
+			h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
+			K[0][0]=h/(2.*vel)*fabs(vx*vx);  K[0][1]=h/(2.*vel)*fabs(vx*vy); K[0][2]=h/(2.*vel)*fabs(vx*vz);
+			K[1][0]=h/(2.*vel)*fabs(vy*vx);  K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz);
+			K[2][0]=h/(2.*vel)*fabs(vz*vx);  K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
+			for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
+
+			GetBAdvecprime(Bprime,element,xyz_list,gauss); 
+
+			TripleMultiply(Bprime,3,numnodes,1,
+						&K[0][0],3,3,0,
+						Bprime,3,numnodes,0,
+						&Ke->values[0],1);
+		}
+		else if(stabilization==2){
+			element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa/rho_ice);
+			for(int i=0;i<numnodes;i++){
+				for(int j=0;j<numnodes;j++){
+					Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
+					  ((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i])*((u-um)*dbasis[0*3+j]+(v-vm)*dbasis[1*3+j]+(w-wm)*dbasis[2*3+j]);
+				}
+			}
+			if(dt!=0.){
+				for(int i=0;i<numnodes;i++){
+					for(int j=0;j<numnodes;j++){
+						Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i]);
+					}
+				}
+			}
+		}
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(dbasis);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
+	delete gauss;
+	return Ke;
+}/*}}}*/
+ElementMatrix* EnthalpyAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/
+
+	/*Initialize Element matrix and return if necessary*/
+	if(!element->IsOnBed() || !element->IsFloating()) return NULL;
+
+	IssmDouble  dt,Jdet,D;
+	IssmDouble *xyz_list_base = NULL;
+
+	/*Get basal element*/
+	if(!element->IsOnBed() || !element->IsFloating()) return NULL;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize vectors*/
+	ElementMatrix* Ke    = element->NewElementMatrix();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	IssmDouble  gravity             = element->GetMaterialParameter(ConstantsGEnum);
+	IssmDouble  rho_water           = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+	IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
+	IssmDouble  mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
+	IssmDouble  thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussBase(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
+		if(reCast<bool,IssmDouble>(dt)) D=dt*D;
+		TripleMultiply(basis,numnodes,1,0,
+					&D,1,1,0,
+					basis,1,numnodes,0,
+					&Ke->values[0],1);
+
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list_base);
+	return Ke;
 }/*}}}*/
 ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/
@@ -289,13 +480,97 @@
 }/*}}}*/
 ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/
-	return NULL;
+
+	/* Geothermal flux on ice sheet base and basal friction */
+	if(!element->IsOnBed() || element->IsFloating()) return NULL;
+
+	IssmDouble  dt,Jdet,enthalpy,pressure,watercolumn,geothermalflux,vx,vy,vz;
+	IssmDouble  enthalpyup,pressureup,alpha2,scalar,basalfriction,heatflux;
+	IssmDouble *xyz_list_base = NULL;
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize vectors*/
+	ElementVector* pe    = element->NewElementVector();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	Input* vx_input             = element->GetInput(VxEnum);                          _assert_(vx_input);
+	Input* vy_input             = element->GetInput(VyEnum);                          _assert_(vy_input);
+	Input* vz_input             = element->GetInput(VzEnum);                          _assert_(vz_input);
+	Input* enthalpy_input       = element->GetInput(EnthalpyPicardEnum);              _assert_(enthalpy_input);
+	Input* pressure_input       = element->GetInput(PressureEnum);                    _assert_(pressure_input);
+	Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
+	Input* watercolumn_input    = element->GetInput(WatercolumnEnum);                 _assert_(watercolumn_input);
+	IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
+
+	/*Build friction element, needed later: */
+	Friction* friction=new Friction(element,3);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss   = element->NewGaussBase(2);
+	Gauss* gaussup = element->NewGaussTop(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		element->NodalFunctions(basis,gauss);
+
+		enthalpy_input->GetInputValue(&enthalpy,gauss);
+		pressure_input->GetInputValue(&pressure,gauss);
+		watercolumn_input->GetInputValue(&watercolumn,gauss);
+
+		if((watercolumn<=0.) && (enthalpy < PureIceEnthalpy(element,pressure))){
+			/* the above check is equivalent to 
+			 NOT ((watercolumn>0.) AND (enthalpy<PIE)) AND (enthalpy<PIE)*/
+			geothermalflux_input->GetInputValue(&geothermalflux,gauss);
+
+			friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
+			vx_input->GetInputValue(&vx,gauss);
+			vy_input->GetInputValue(&vy,gauss);
+			vz_input->GetInputValue(&vz,gauss);
+			basalfriction = alpha2*(vx*vx + vy*vy + vz*vz);
+			heatflux      = (basalfriction+geothermalflux)/(rho_ice);
+
+			scalar = gauss->weight*Jdet*heatflux;
+			if(dt!=0.) scalar=dt*scalar;
+
+			for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
+		}
+		else if(enthalpy >= PureIceEnthalpy(element,pressure)){
+			/* check positive thickness of temperate basal ice layer */
+			enthalpy_input->GetInputValue(&enthalpyup,gaussup);
+			pressure_input->GetInputValue(&pressureup,gaussup);
+			if(enthalpyup >= PureIceEnthalpy(element,pressureup)){
+				// TODO: temperate ice has positive thickness: grad enthalpy*n=0.
+			}
+			else{
+				// only base temperate, set Dirichlet BCs in Penta::UpdateBasalConstraintsEnthalpy()
+			}
+		}
+		else{
+			// base cold, but watercolumn positive. Set base to h_pmp.
+		}
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	delete gaussup;
+	delete friction;
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list_base);
+	return pe;
+
 }/*}}}*/
 ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/
+
+	/*Get basal element*/
+	if(!element->IsOnBed() || !element->IsFloating()) return NULL;
 
 	IssmDouble  h_pmp,dt,Jdet,scalar_ocean,pressure;
 	IssmDouble *xyz_list_base = NULL;
-
-	/*Get basal element*/
-	if(!element->IsOnBed() || !element->IsFloating()) return NULL;
 
 	/*Fetch number of nodes for this finite element*/
@@ -339,4 +614,91 @@
 	xDelete<IssmDouble>(xyz_list_base);
 	return pe;
+}/*}}}*/
+void EnthalpyAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 
+	 * For node i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_conduct=[ dh/dx ]
+	 *                  [ dh/dy ]
+	 *                  [ dh/dz ]
+	 * where h is the interpolation function for node i.
+	 *
+	 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		B[numnodes*0+i] = dbasis[0*numnodes+i];
+		B[numnodes*1+i] = dbasis[1*numnodes+i];
+		B[numnodes*2+i] = dbasis[2*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+}/*}}}*/
+void EnthalpyAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 
+	 * For node i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_advec =[ h ]
+	 *                 [ h ]
+	 *                 [ h ]
+	 * where h is the interpolation function for node i.
+	 *
+	 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	element->NodalFunctions(basis,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		B[numnodes*0+i] = basis[i];
+		B[numnodes*1+i] = basis[i];
+		B[numnodes*2+i] = basis[i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
+}/*}}}*/
+void EnthalpyAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 
+	 * For node i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Biprime_advec=[ dh/dx ]
+	 *                     [ dh/dy ]
+	 *                     [ dh/dz ]
+	 * where h is the interpolation function for node i.
+	 *
+	 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		B[numnodes*0+i] = dbasis[0*numnodes+i];
+		B[numnodes*1+i] = dbasis[1*numnodes+i];
+		B[numnodes*2+i] = dbasis[2*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 void EnthalpyAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
@@ -422,2 +784,88 @@
 	xDelete<int>(doflist);
 }/*}}}*/
+
+/*Intermediaries*/
+IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
+
+	IssmDouble heatcapacity             = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
+	IssmDouble temperateiceconductivity = element->GetMaterialParameter(MaterialsTemperateiceconductivityEnum);
+	IssmDouble thermalconductivity      = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
+
+	if(enthalpy < PureIceEnthalpy(element,pressure)){
+		return thermalconductivity/heatcapacity;
+	}
+	else{
+		return temperateiceconductivity/heatcapacity;
+	}
+}/*}}}*/
+IssmDouble EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element){/*{{{*/
+
+	int         iv;
+	IssmDouble  lambda;                   /* fraction of cold ice    */
+	IssmDouble  kappa   ,kappa_c,kappa_t; /* enthalpy conductivities */
+	IssmDouble  Hc,Ht;
+
+
+	/*Get pressures and enthalpies on vertices*/
+	int         numvertices = element->GetNumberOfVertices();
+	IssmDouble* pressures   = xNew<IssmDouble>(numvertices);
+	IssmDouble* enthalpies  = xNew<IssmDouble>(numvertices);
+	IssmDouble* PIE         = xNew<IssmDouble>(numvertices);
+	IssmDouble* dHpmp       = xNew<IssmDouble>(numvertices);
+	element->GetInputListOnVertices(pressures,PressureEnum);
+	element->GetInputListOnVertices(enthalpies,EnthalpyEnum);
+	for(iv=0;iv<numvertices;iv++){
+		PIE[iv]   = PureIceEnthalpy(element,pressures[iv]);
+		dHpmp[iv] = enthalpies[iv]-PIE[iv];
+	}
+
+	bool allequalsign = true;
+	if(dHpmp[0]<0.){
+		for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0.));
+	}
+	else{
+		for(iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0.));
+	}
+
+	if(allequalsign){
+		kappa = EnthalpyDiffusionParameter(element,enthalpies[0],pressures[0]);
+	}
+	else{
+		/* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
+			cf Patankar 1980, pp44 */
+		kappa_c = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)-1.,0.);
+		kappa_t = EnthalpyDiffusionParameter(element,PureIceEnthalpy(element,0.)+1.,0.);
+		Hc=0.; Ht=0.;
+		for(iv=0; iv<numvertices;iv++){
+			if(enthalpies[iv]<PIE[iv])
+			 Hc+=(PIE[iv]-enthalpies[iv]);
+			else
+			 Ht+=(enthalpies[iv]-PIE[iv]);
+		}
+		_assert_((Hc+Ht)>0.);
+		lambda = Hc/(Hc+Ht);
+		kappa  = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(PIE);
+	xDelete<IssmDouble>(dHpmp);
+	xDelete<IssmDouble>(pressures);
+	xDelete<IssmDouble>(enthalpies);
+	return kappa;
+}
+/*}}}*/
+IssmDouble EnthalpyAnalysis::PureIceEnthalpy(Element* element,IssmDouble pressure){/*{{{*/
+
+	IssmDouble heatcapacity         = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
+	IssmDouble referencetemperature = element->GetMaterialParameter(ConstantsReferencetemperatureEnum);
+
+	return heatcapacity*(TMeltingPoint(element,pressure)-referencetemperature);
+}/*}}}*/
+IssmDouble EnthalpyAnalysis::TMeltingPoint(Element* element,IssmDouble pressure){/*{{{*/
+
+	IssmDouble meltingpoint = element->GetMaterialParameter(MaterialsMeltingpointEnum);
+	IssmDouble beta         = element->GetMaterialParameter(MaterialsBetaEnum);
+
+	return meltingpoint-beta*pressure;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 16887)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h	(revision 16888)
@@ -22,10 +22,21 @@
 		/*Finite element Analysis*/
 		ElementMatrix* CreateKMatrix(Element* element);
+		ElementMatrix* CreateKMatrixVolume(Element* element);
+		ElementMatrix* CreateKMatrixShelf(Element* element);
 		ElementVector* CreatePVector(Element* element);
 		ElementVector* CreatePVectorVolume(Element* element);
 		ElementVector* CreatePVectorSheet(Element* element);
 		ElementVector* CreatePVectorShelf(Element* element);
+		void GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
 		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
+
+		/*Intermediaries*/
+		IssmDouble EnthalpyDiffusionParameter(Element* element,IssmDouble enthalpy,IssmDouble pressure);
+		IssmDouble EnthalpyDiffusionParameterVolume(Element* element);
+		IssmDouble PureIceEnthalpy(Element* element,IssmDouble pressure);
+		IssmDouble TMeltingPoint(Element* element,IssmDouble pressure);
 };
 #endif