Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 22469)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 22470)
@@ -811,7 +811,7 @@
 
 		/*Update signed distance*/
 		dmin = sqrt(dmin);
-		if(dmin>10000) dmin=10000;
+		// if(dmin>10000) dmin=10000;
 		if(ls[j]>0){
 			ls[j] = dmin;
 		}
@@ -2736,6 +2736,99 @@
 	return TriaRef::NumberofNodes(this->VelocityInterpolation());
 }
 /*}}}*/
+void       Tria::PicoUpdateBoxid(int* pmax_boxid_basin){/*{{{*/
+
+	if(!this->IsIceInElement() || !this->IsFloating()) return;
+
+	//Load inputs and params
+	int i, boxid;
+	int basin_id, boxid_max;
+	IssmDouble dist_gl;
+	IssmDouble dist_cf;
+	IssmDouble rel_dist_gl;
+
+	inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum);
+	boxid_max=pmax_boxid_basin[basin_id]; //0-based
+
+	Input* dist_gl_input=NULL;
+	Input* dist_cf_input=NULL;
+	dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
+	dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum); _assert_(dist_cf_input);
+
+	Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
+	dist_gl_input->GetInputValue(&dist_gl,gauss);
+	dist_cf_input->GetInputValue(&dist_cf,gauss);
+
+	rel_dist_gl=dist_gl/(dist_gl+dist_cf);
+
+	boxid=-1;
+	for(i=0;i<boxid_max;i++){
+		if(rel_dist_gl>=(1-sqrt((boxid_max-i)/boxid_max)) && rel_dist_gl<=(1-sqrt((boxid_max-i-1)/boxid_max))){
+			boxid=i; break;
+		}
+	}
+	if(boxid==-1){_error_("no boxid found for element" << this->Sid());}
+
+	this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid));	
+
+	/*Cleanup & return: */
+	delete gauss;
+}
+/*}}}*/
+void       Tria::PicoUpdateFirstBox(){/*{{{*/
+
+	IssmDouble rhoi, rhow, earth_grav, rho_star, nu, latentheat, c_p_ocean, lambda, a, b, c;
+	IssmDouble alpha, Beta, gamma_T, overturning_coeff, t_farocean, s_farocean;
+	IssmDouble basinid, boxid, thickness, toc_farocean, soc_farocean, area_box1;
+	IssmDouble pressure, T_star, g1, s1, p_coeff, q_coeff, Toc, Soc, potential_pressure_melting_point;
+	IssmDouble basalmeltrate_shelf, overturning, maxbox;
+
+	//Get variables
+	rhoi			= this->GetMaterialParameter(MaterialsRhoIceEnum);
+	rhow			= this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
+	earth_grav  = this->GetMaterialParameter(ConstantsGEnum);
+	rho_star		= 1033;             //kg/m^3
+	nu				= rhoi/rhow;
+	latentheat	= this->GetMaterialParameter(MaterialsLatentheatEnum);
+	c_p_ocean	= this->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
+	lambda		= latentheat/c_p_ocean;
+	a				= -0.0572;          //K/psu
+	b				= 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum);  //K
+	c				= this->GetMaterialParameter(MaterialsBetaEnum);
+	alpha			= 0.000075;         //1/K
+	Beta			= 0.00077;          //K
+
+	this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
+	this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);
+	this->parameters->FindParam(&t_farocean,BasalforcingsPicoFarOceanTemperatureEnum);
+	this->parameters->FindParam(&s_farocean,BasalforcingsPicoFarOceanSalinityEnum);
+
+	this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+	this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
+	this->inputs->GetInputValue(&maxbox,BasalforcingsPicoMaxboxcountEnum);
+	this->inputs->GetInputValue(&thickness,ThicknessEnum);
+
+	toc_farocean = t_farocean[basinid];
+	soc_farocean = s_farocean[basinid];
+	area_box1 = areas[basinid*maxbox+boxid];
+	pressure = (rhoi*earth_grav)*thickness;
+	T_star = a*soc_farocean+b-c*pressure-toc_farocean;
+	g1 = area_box1*gamma_T;
+	s1 = soc_farocean/(nu*lambda);
+	p_coeff = g1/(overturning_coeff*rho_star*(Beta*s1-alpha));
+	q_coeff = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha)));
+	if ((0.25*(p_coeff)^(2) - q_coeff) < 0)
+	  {q_coeff = (0.25*(p_coeff)^(2))}
+	Toc = toc_farocean-(-0.5*p_coeff+sqrt(0.25*p_coeff^(2)-q_coeff));
+	Soc = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Toc);
+	potential_pressure_melting_point = a*Soc+b-c*pressure;
+	basalmeltrate_shelf = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point-Toc);
+	overturning = overturning_coeff*rho_star*(Beta*(soc_farocean-Soc)-alpha*(toc_farocean-Toc));
+
+	this->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,this->GetElementType());
+
+}
+/*}}}*/
 void       Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
 
 	IssmDouble  h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
@@ -2967,7 +3060,7 @@
 	/*Get out if this is not an element input*/
 	if(!IsInput(control_enum)) return;
 
-	/*Prepare index list*/
+	/*hrepare index list*/
 	GradientIndexing(&vertexpidlist[0],control_index);
 
 	/*Get values on vertices*/