Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 20046)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 20047)
@@ -376,4 +376,5 @@
 
 	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
 	delete gauss;
 	return divergence;
@@ -1574,7 +1575,8 @@
 	IssmDouble  mantleconductivity,nusselt,dtbg,plumeradius,topplumedepth,bottomplumedepth,plumex,plumey;
 	IssmDouble  crustthickness,uppercrustthickness,uppercrustheat,lowercrustheat;
-	//IssmDouble* base     = xNew<IssmDouble>(numvertices);
-	//IssmDouble* bed      = xNew<IssmDouble>(numvertices);
+	IssmDouble  crustheat,plumeheat,dt,middleplumedepth,a,e,eprime,A0,lambda,Alambda,dAlambda;
+	IssmDouble  x,y,z,c;
 	IssmDouble* values   = xNew<IssmDouble>(numvertices);
+	IssmDouble *xyz_list = NULL;
 
 	parameters->FindParam(&mantleconductivity,BasalforcingsMantleconductivityEnum);
@@ -1591,13 +1593,26 @@
 	parameters->FindParam(&lowercrustheat,BasalforcingsLowercrustheatEnum);
 
-	//this->GetInputListOnVertices(base,BaseEnum);
-	//this->GetInputListOnVertices(bed,BedEnum);
+	this->GetVerticesCoordinates(&xyz_list);
+	c=plumeradius;
+	a=(bottomplumedepth-topplumedepth)/2.;
+	e=pow(a*a-c*c,1./2.)/a;
+	A0=(1-e*e)/(e*e*e)*(1/2*log((1+e)/(1-e))-e);
+	middleplumedepth=-topplumedepth-crustthickness-a;
 	for(int i=0;i<numvertices;i++){
-		values[i]=uppercrustheat*uppercrustthickness+ lowercrustheat*(crustthickness-uppercrustthickness);
+		y=xyz_list[i*3+0]-plumex;
+		z=xyz_list[i*3+1]-plumey;
+		x=(a+topplumedepth+crustthickness);
+		lambda=(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))/2;
+		dAlambda=(-8*a*pow(c,2)*x*(-2*pow(a,2)+2*pow(c,2)+sqrt(2)*sqrt((a-c)*(a+c))*sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))*(pow(a,4)*(pow(y,2)+pow(z,2))+pow(c,4)*(pow(y,2)+pow(z,2))+pow(pow(x,2)+pow(y,2)+pow(z,2),2)*(pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))+pow(c,2)*(pow(x,4)-pow(x,2)*(pow(y,2)+pow(z,2))-(pow(y,2)+pow(z,2))*(2*pow(y,2)+2*pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))+pow(a,2)*(-pow(x,4)+pow(x,2)*(pow(y,2)+pow(z,2))+(pow(y,2)+pow(z,2))*(-2*pow(c,2)+2*(pow(y,2)+pow(z,2))+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))))/(sqrt((a-c)*(a+c))*sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))*pow(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))),3.5)*pow(-(sqrt(2)*sqrt((a-c)*(a+c)))+sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))),2)*(sqrt(2)*sqrt((a-c)*(a+c))+sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))));
+		eprime=pow((a*a-plumeradius*plumeradius)/(a*a+lambda),1./2.);
+		Alambda=(1.-e*e)/(e*e*e)*(1./2.*log((1.+eprime)/(1.-eprime))-eprime);
+		dt=dtbg-(nusselt-1.)/(1.+A0*(nusselt-1.))*(Alambda*dtbg+middleplumedepth*dtbg*dAlambda);
+		plumeheat=mantleconductivity*dt;
+		crustheat=uppercrustheat*uppercrustthickness+lowercrustheat*(crustthickness-uppercrustthickness);
+		values[i]=crustheat+plumeheat;
 	}
 
 	this->AddInput(BasalforcingsGeothermalfluxEnum,values,P1Enum);
-	//xDelete<IssmDouble>(base);
-	//xDelete<IssmDouble>(bed);
+	xDelete<IssmDouble>(xyz_list);
 	xDelete<IssmDouble>(values);
 
@@ -2669,4 +2684,5 @@
 	/*Clean up and return*/
 	xDelete<IssmDouble>(maxprincipal);
+	xDelete<IssmDouble>(xyz_list);
 	delete gauss;
 }
@@ -3098,4 +3114,5 @@
 	/*Clean up and return*/
 	xDelete<IssmDouble>(viscousheating);
+	xDelete<IssmDouble>(xyz_list);
 	delete gauss;
 }