Index: issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp	(revision 24002)
+++ issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp	(revision 24010)
@@ -17,5 +17,7 @@
 
 	/*First fetch data: */
+	if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
 	::CreateNodes(nodes,iomodel,GLheightadvectionAnalysisEnum,P1Enum);
+	iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
 }/*}}}*/
 int  GLheightadvectionAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
@@ -56,7 +58,10 @@
 ElementMatrix* GLheightadvectionAnalysis::CreateKMatrix(Element* element){/*{{{*/
 
+	/* Spawn basal element */
+	if(!element->IsOnBase()) return NULL;
+	Element* basalelement = element->SpawnBasalElement();
+
 	/* Check if ice in element */
-	if(!element->IsIceInElement()) return NULL;
-	//if(!element->IsFloating()) return NULL;
+	if(!basalelement->IsIceInElement()) return NULL;
 
 	/*Intermediaries */
@@ -65,8 +70,11 @@
 	IssmDouble Jdet,D_scalar,onboundary;
 	IssmDouble vel,vx,vy;
-	IssmDouble* xyz_list = NULL;
+	IssmDouble* xyz_list      = NULL;
+	Input* vx_input           = NULL;
+	Input* vy_input           = NULL;
+	Input* bc_input           = NULL;
 
 	/*Get problem dimension*/
-	element->FindParam(&domaintype,DomainTypeEnum);
+	basalelement->FindParam(&domaintype,DomainTypeEnum);
 	switch(domaintype){
 		case Domain2DhorizontalEnum: dim = 2; break;
@@ -76,8 +84,8 @@
 
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
+	int numnodes = basalelement->GetNumberOfNodes();
 
 	/*Initialize Element vector and other vectors*/
-	ElementMatrix* Ke     = element->NewElementMatrix();
+	ElementMatrix* Ke     = basalelement->NewElementMatrix();
 	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 	IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
@@ -86,22 +94,29 @@
 
 	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	//Input* vx_input=element->GetInput(BaseSlopeXEnum); _assert_(vx_input);
-	//Input* vy_input=element->GetInput(BaseSlopeYEnum); _assert_(vy_input);
-	Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
-	Input* bc_input=element->GetInput(MeshVertexonboundaryEnum); _assert_(bc_input);
-
-	IssmDouble h = element->CharacteristicLength();
+	basalelement->GetVerticesCoordinates(&xyz_list);
+	switch(domaintype){
+		case Domain2DhorizontalEnum:
+			vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
+			vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
+		   bc_input=basalelement->GetInput(MeshVertexonboundaryEnum); _assert_(bc_input);
+		break;
+		case Domain3DEnum:
+			vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
+			vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
+			bc_input=basalelement->GetInput(MeshVertexonboundaryEnum); _assert_(bc_input);
+		break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+	}
+	IssmDouble h = basalelement->CharacteristicLength();
 
 	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGauss(4);
+	Gauss* gauss=basalelement->NewGauss(4);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		element->NodalFunctions(basis,gauss);
-		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-		GetBprime(Bprime,element,dim,xyz_list,gauss);
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctions(basis,gauss);
+		basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+		GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
 
 		/*Get velocity*/
@@ -168,9 +183,9 @@
 	xDelete<IssmDouble>(Bprime);
 	delete gauss;
+	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return Ke;
 
 }/*}}}*/
 ElementVector* GLheightadvectionAnalysis::CreatePVector(Element* element){/*{{{*/
-
 	return NULL;
 }/*}}}*/
@@ -203,5 +218,4 @@
 	for(int i=0;i<femmodel->elements->Size();i++){
 		Element    *element  = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-
 		int         numnodes  = element->GetNumberOfNodes();
 		IssmDouble *mask      = xNew<IssmDouble>(numnodes);
Index: issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24002)
+++ issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24010)
@@ -2157,4 +2157,61 @@
 }
 /*}}}*/
+void       Element::Ismip6FloatingiceMeltingRate(){/*{{{*/
+
+	if(!this->IsIceInElement() || !this->IsFloating()) return;
+
+	int         numvertices = this->GetNumberOfVertices();
+	int         basinid,num_basins,M,N;
+	IssmDouble  tf,gamma0,base,delta_t_basin,mean_tf_basin,absval;
+	IssmDouble  basalmeltrate[numvertices];
+	bool        islocal;
+	IssmDouble* delta_t = NULL;
+	IssmDouble* mean_tf = NULL;
+	IssmDouble* depths  = NULL;
+
+	/*Get variables*/
+	IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum);
+	IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum);
+	IssmDouble lf   = this->FindParam(MaterialsLatentheatEnum);
+	IssmDouble cp   = this->FindParam(MaterialsMixedLayerCapacityEnum);
+
+	/*Hard code sea water density to be consistent with ISMIP6 documentation*/
+	rhow = 1028.;
+
+	/* Get parameters and inputs */
+	this->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum);
+	this->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum);
+	this->parameters->FindParam(&gamma0,BasalforcingsIsmp6Gamma0Enum);
+	this->parameters->FindParam(&delta_t,&M,BasalforcingsIsmp6DeltaTEnum);    _assert_(M==num_basins);
+	this->parameters->FindParam(&islocal,BasalforcingsIsmp6IsLocalEnum);
+	if(!islocal) {
+		this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmp6AverageTfEnum); _assert_(N==num_basins);
+	}
+	Input* tf_input=this->GetInput(BasalforcingsIsmp6TfShelfEnum);            _assert_(tf_input);
+	delta_t_basin = delta_t[basinid];
+	if(!islocal) mean_tf_basin = mean_tf[basinid];
+
+	/*Compute melt rate for Local and Nonlocal parameterizations*/
+	Gauss* gauss=this->NewGauss();
+	for(int i=0;i<numvertices;i++){
+		gauss->GaussVertex(i);
+		tf_input->GetInputValue(&tf,gauss);
+		if(!islocal) {
+			absval = mean_tf_basin+delta_t_basin;
+			if (absval<0) {absval = -1.*absval;}
+			basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval;}
+		else {basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*pow(max(tf+delta_t_basin,0.),2);}
+	}
+
+	/*Return basal melt rate*/
+	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrate,P1Enum);
+
+	/*Cleanup and return*/
+	delete gauss;
+	xDelete<IssmDouble>(delta_t);
+	xDelete<IssmDouble>(mean_tf);
+	xDelete<IssmDouble>(depths);
+
+}/*}}}*/
 bool       Element::IsWaterInElement(){/*{{{*/
 	return (this->inputs->Max(MaskOceanLevelsetEnum)>0.);
@@ -2520,4 +2577,280 @@
 }
 /*}}}*/
+void       Element::PicoUpdateBoxid(int* max_boxid_basin_list){/*{{{*/
+
+	if(!this->IsIceInElement() || !this->IsFloating()) return;
+
+	int        basin_id;
+	IssmDouble dist_gl,dist_cf;
+
+	inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum);
+	IssmDouble boxid_max=reCast<IssmDouble>(max_boxid_basin_list[basin_id])+1.;
+
+	Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
+	Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum);  _assert_(dist_cf_input);
+
+	/*Get dist_gl and dist_cf at center of element*/
+	Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
+	dist_gl_input->GetInputValue(&dist_gl,gauss);
+	dist_cf_input->GetInputValue(&dist_cf,gauss);
+	delete gauss;
+
+	/*Ensure values are positive for floating ice*/
+	dist_gl = fabs(dist_gl);
+	dist_cf = fabs(dist_cf);
+
+	/*Compute relative distance to grounding line*/
+	IssmDouble rel_dist_gl=dist_gl/(dist_gl+dist_cf);
+
+	/*Assign box numbers based on rel_dist_gl*/
+	int boxid = -1;
+	for(IssmDouble i=0.;i<boxid_max;i++){
+		IssmDouble lowbound  = 1. -sqrt((boxid_max-i   )/boxid_max);
+		IssmDouble highbound = 1. -sqrt((boxid_max-i-1.)/boxid_max);
+		if(rel_dist_gl>=lowbound && rel_dist_gl<=highbound){
+			boxid=reCast<int>(i);
+			break;
+		}
+	}
+	if(boxid==-1) _error_("No boxid found for element " << this->Sid() << "!");
+
+	this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid));
+	
+}/*}}}*/
+void       Element::PicoUpdateBox(int loopboxid){/*{{{*/
+
+	if(!this->IsIceInElement() || !this->IsFloating()) return;
+
+	int boxid;
+	this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
+	if(loopboxid!=boxid) return;
+
+	int        NUMVERTICES = this->GetNumberOfVertices();
+	int        basinid, maxbox, num_basins, numnodes, M;
+	IssmDouble gamma_T, overturning_coeff, thickness;
+	IssmDouble pressure, T_star,p_coeff, q_coeff;
+	bool       isplume;
+
+	/*Get variables*/
+	IssmDouble rhoi       = this->FindParam(MaterialsRhoIceEnum);
+	IssmDouble rhow       = this->FindParam(MaterialsRhoSeawaterEnum);
+	IssmDouble earth_grav = this->FindParam(ConstantsGEnum);
+	IssmDouble rho_star   = 1033.;             // kg/m^3
+	IssmDouble nu         = rhoi/rhow;
+	IssmDouble latentheat = this->FindParam(MaterialsLatentheatEnum);
+	IssmDouble c_p_ocean  = this->FindParam(MaterialsMixedLayerCapacityEnum);
+	IssmDouble lambda     = latentheat/c_p_ocean;
+	IssmDouble a          = -0.0572;          // K/psu
+	IssmDouble b          = 0.0788 + this->FindParam(MaterialsMeltingpointEnum);  //K
+	IssmDouble c          = 7.77e-4;
+	IssmDouble alpha      = 7.5e-5;           // 1/K
+	IssmDouble Beta       = 7.7e-4;           // K
+
+	/* Get non-box-specific parameters and inputs */
+	this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum);
+	this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
+	this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);
+	this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
+	this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+	this->parameters->FindParam(&isplume, BasalforcingsPicoIsplumeEnum);
+	Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
+	_assert_(basinid<=num_basins);
+
+	IssmDouble* boxareas = NULL;
+	this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
+	_assert_(M==num_basins*maxbox);
+
+	IssmDouble area_boxi        = boxareas[basinid*maxbox+boxid];
+	IssmDouble g1               = area_boxi*gamma_T;
+
+	IssmDouble basalmeltrates_shelf[NUMVERTICES];
+	IssmDouble potential_pressure_melting_point[NUMVERTICES];
+	IssmDouble Tocs[NUMVERTICES];
+	IssmDouble Socs[NUMVERTICES];
+
+	/* First box calculations */
+	if(boxid==0){
+		/* Get box1 parameters and inputs */
+		IssmDouble time, toc_farocean, soc_farocean;
+		this->parameters->FindParam(&time,TimeEnum);
+		this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum);
+		this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);
+		IssmDouble s1 = soc_farocean/(nu*lambda);
+		IssmDouble overturnings[NUMVERTICES];
+
+		/* Start looping on the number of verticies and calculate ocean vars */
+		Gauss* gauss=this->NewGauss();
+		for(int i=0;i<NUMVERTICES;i++){
+			gauss->GaussVertex(i);
+			thickness_input->GetInputValue(&thickness,gauss);
+			pressure = (rhoi*earth_grav*1e-4)*thickness;
+			T_star   = a*soc_farocean+b-c*pressure-toc_farocean;
+			p_coeff  = g1/(overturning_coeff*rho_star*(Beta*s1-alpha));
+			q_coeff  = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha)));
+
+			/* To avoid negatives under the square root */
+			if((0.25*pow(p_coeff,2)-q_coeff)<0) q_coeff = 0.25*p_coeff*p_coeff;
+
+			Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff));
+			Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]);
+			potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
+			if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
+			overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i]));
+		}
+
+		if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
+		this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);
+		this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);
+		this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1Enum);
+
+		/*Cleanup and return*/
+		delete gauss;
+	}
+
+	/* Subsequent box calculations */
+	else {
+		/* Get subsequent box parameters and inputs */
+		IssmDouble* toc_weighted_avg         = NULL;
+		IssmDouble* soc_weighted_avg         = NULL;
+		IssmDouble* overturning_weighted_avg = NULL;
+		this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum);
+		this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum);
+		this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum);
+		IssmDouble mean_toc                  = toc_weighted_avg[basinid];
+		IssmDouble mean_soc                  = soc_weighted_avg[basinid];
+		IssmDouble mean_overturning          = overturning_weighted_avg[basinid];
+		IssmDouble g2                        = g1/(nu*lambda);
+
+		/* Start looping on the number of verticies and calculate ocean vars */
+		Gauss* gauss=this->NewGauss();
+		for(int i=0;i<NUMVERTICES;i++){
+			gauss->GaussVertex(i);
+			thickness_input->GetInputValue(&thickness,gauss);
+			pressure = (rhoi*earth_grav*1.e-4)*thickness;
+			T_star   = a*mean_soc+b-c*pressure-mean_toc;
+			Tocs[i]  = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc));
+			Socs[i]  = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda));
+			potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
+			if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
+		}
+
+		if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
+		this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);
+		this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);
+
+		/*Cleanup and return*/
+		xDelete<IssmDouble>(toc_weighted_avg);
+		xDelete<IssmDouble>(soc_weighted_avg);
+		xDelete<IssmDouble>(overturning_weighted_avg);
+		delete gauss;
+	}
+
+	/*Cleanup and return*/
+	xDelete<IssmDouble>(boxareas);
+	
+}/*}}}*/
+void       Element::PicoComputeBasalMelt(){/*{{{*/
+
+	if(!this->IsIceInElement() || !this->IsFloating()) return;
+
+	int        NUMVERTICES = this->GetNumberOfVertices();
+	IssmDouble E0, Cd, CdT, YT, lam1, lam2, lam3, M0, CdTS0, y1, y2, x0;
+	IssmDouble Tf_gl, YTS, CdTS, G1, G2, G3, g_alpha, M, l, X_hat, M_hat;
+	IssmDouble alpha, zgl, Toc, Soc, z_base, yts, slopex, slopey;
+
+	/*Get variables*/
+	E0    = 3.6e-2;        //Entrainment coefficient
+	Cd    = 2.5e-3;        //Drag coefficient
+	CdT   = 1.1e-3;        //Turbulent heat exchange coefficient
+	YT    = CdT/sqrt(Cd);  //Heat exchange coefficient
+	lam1  = -5.73e-2;      //Freezing point-salinity coefficient (degrees C)
+	lam2  = 8.32e-2;       //Freezing point offset (degrees C)
+	lam3  = 7.61e-4;       //Freezing point-depth coefficient (K m-1)
+	M0    = 10.;           //Melt-rate parameter (m yr-1 C-2)
+	CdTS0 = 6e-4;          //Heat exchange parameter
+	y1    = 0.545;         //Heat exchange parameter
+	y2    = 3.5e-5;        //Heat exchange parameter
+	x0    = 0.56;          //Dimentionless scaling factor
+
+	/*Define arrays*/
+	IssmDouble basalmeltrates_shelf[NUMVERTICES];  //Basal melt-rate
+
+	/*Polynomial coefficients*/
+	IssmDouble p[12];
+	p[0]  =  0.1371330075095435;
+	p[1]  =  5.527656234709359E1;
+	p[2]  = -8.951812433987858E2;
+	p[3]  =  8.927093637594877E3;
+	p[4]  = -5.563863123811898E4;
+	p[5]  =  2.218596970948727E5;
+	p[6]  = -5.820015295669482E5;
+	p[7]  =  1.015475347943186E6;
+	p[8]  = -1.166290429178556E6;
+	p[9]  =  8.466870335320488E5;
+	p[10] = -3.520598035764990E5;
+	p[11] =  6.387953795485420E4;
+
+	/*Get inputs*/
+	Input* zgl_input            = this->GetInput(GroundinglineHeightEnum);                     _assert_(zgl_input);
+	Input* toc_input            = this->GetInput(BasalforcingsPicoSubShelfOceanTempEnum);      _assert_(toc_input);
+	Input* soc_input            = this->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum);  _assert_(soc_input);
+	Input* base_input           = this->GetInput(BaseEnum);                                    _assert_(base_input);
+	Input* baseslopex_input     = this->GetInput(BaseSlopeXEnum);                              _assert_(baseslopex_input);
+	Input* baseslopey_input     = this->GetInput(BaseSlopeYEnum);                              _assert_(baseslopey_input);
+	this->FindParam(&yts, ConstantsYtsEnum);
+
+	/*Loop over the number of vertices in this element*/
+	Gauss* gauss=this->NewGauss();
+	for(int i=0;i<NUMVERTICES;i++){
+		gauss->GaussVertex(i);
+
+		/*Get inputs*/
+		zgl_input->GetInputValue(&zgl,gauss);
+		toc_input->GetInputValue(&Toc,gauss); //(K)
+		soc_input->GetInputValue(&Soc,gauss);
+		base_input->GetInputValue(&z_base,gauss);
+		baseslopex_input->GetInputValue(&slopex,gauss);
+		baseslopey_input->GetInputValue(&slopey,gauss);
+
+		/*Compute ice shelf base slope (radians)*/
+		alpha = atan(sqrt(slopex*slopex + slopey*slopey));
+		if(alpha>=M_PI) alpha = M_PI - 0.001;               //ensure sin(alpha) > 0 for meltrate calculations
+
+		/*Make necessary conversions*/
+		Toc = Toc-273.15;
+		if(zgl>z_base) zgl=z_base;
+
+		/*Low bound for Toc to ensure X_hat is between 0 and 1*/
+		if(Toc<lam1*Soc+lam2) Toc=lam1*Soc+lam2;
+
+		/*Compute parameters needed for melt-rate calculation*/
+		Tf_gl = lam1*Soc+lam2+lam3*zgl;                                              //Characteristic freezing point
+		YTS = YT*(y1+y2*(((Toc-Tf_gl)*E0*sin(alpha))/(lam3*(CdTS0+E0*sin(alpha))))); //Effective heat exchange coefficient
+		CdTS = sqrt(Cd)*YTS;                                                         //Heat exchange coefficient
+		G1 = sqrt(sin(alpha)/(Cd+E0*sin(alpha)));                                    //Geometric factor
+		G2 = sqrt(CdTS/(CdTS+E0*sin(alpha)));                                        //Geometric factor
+		G3 = (E0*sin(alpha))/(CdTS+E0*sin(alpha));                                   //Geometric factor
+		g_alpha = G1*G2*G3;                                                          //Melt scaling factor
+		M = M0*g_alpha*pow((Toc-Tf_gl),2);                                           //Melt-rate scale
+		l = ((Toc-Tf_gl)*(x0*CdTS+E0*sin(alpha)))/(lam3*x0*(CdTS+E0*sin(alpha)));    //Length scale
+		X_hat = (z_base-zgl)/l;                                                      //Dimentionless coordinate system
+
+		/*Compute polynomial fit*/
+		M_hat = 0.;                                                                  //Reset summation variable for each node
+		for(int ii=0;ii<12;ii++) {
+			M_hat += p[ii]*pow(X_hat,ii);                                               //Polynomial fit
+		}
+
+		/*Compute melt-rate*/
+		basalmeltrates_shelf[i] = (M*M_hat)/yts;                                     //Basal melt-rate (m/s)
+	}
+
+	/*Save computed melt-rate*/
+	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
+
+	/*Cleanup and return*/
+	delete gauss;
+		
+}/*}}}*/
 void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/
 
Index: issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24002)
+++ issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24010)
@@ -134,4 +134,5 @@
 		bool               IsIceInElement();
 		bool               IsLandInElement();
+		void               Ismip6FloatingiceMeltingRate();
 		bool               IsWaterInElement();
 		void               LinearFloatingiceMeltingRate();
@@ -145,4 +146,7 @@
 		ElementMatrix*     NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum);
 		ElementVector*     NewElementVector(int approximation_enum=NoneApproximationEnum);
+      void               PicoUpdateBoxid(int* pmax_boxid_basin); 
+		void               PicoUpdateBox(int loopboxid);
+		void               PicoComputeBasalMelt(); 
 		void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac);
 		void               PositiveDegreeDaySicopolis(bool isfirnwarming);
@@ -257,5 +261,4 @@
 		virtual bool       IsFaceOnBoundary(void)=0;
 		virtual bool       IsIcefront(void)=0;
-		virtual void       Ismip6FloatingiceMeltingRate(void)=0;
 		virtual bool       IsNodeOnShelfFromFlags(IssmDouble* flags)=0;
 		virtual bool       IsOnBase()=0;
@@ -299,7 +302,4 @@
 		virtual int        NumberofNodesPressure(void)=0;
 		virtual int        NumberofNodesVelocity(void)=0;
-		virtual void       PicoUpdateBoxid(int* pmax_boxid_basin)=0;
-		virtual void       PicoUpdateBox(int loopboxid)=0;
-		virtual void       PicoComputeBasalMelt(void)=0;
 		virtual void       PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
 		virtual int        PressureInterpolation()=0;
Index: issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24002)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24010)
@@ -95,5 +95,4 @@
 		IssmDouble     IceVolume(bool scaled);
 		IssmDouble     IceVolumeAboveFloatation(bool scaled);
-		void				Ismip6FloatingiceMeltingRate(void){_error_("not implemented yet");};
 		void           InputDepthAverageAtBase(int enum_type,int average_enum_type);
 		void	         InputExtrude(int enum_type,int start);
@@ -144,7 +143,4 @@
 		int            NumberofNodesPressure(void);
 		int            NumberofNodesVelocity(void);
-		void				PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");};
-		void				PicoUpdateBox(int loopboxid){_error_("not implemented yet");};
-		void				PicoComputeBasalMelt(void){_error_("not implemented yet");};
 		void           PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
 		int            PressureInterpolation();
Index: issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 24002)
+++ issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 24010)
@@ -88,5 +88,4 @@
 		bool        IsFaceOnBoundary(void){_error_("not implemented yet");};
 		bool		   IsIcefront(void);
-		void		   Ismip6FloatingiceMeltingRate(void){_error_("not implemented yet");};
 		bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
 		bool        IsOnBase(){_error_("not implemented yet");};
@@ -130,7 +129,4 @@
 		int         NumberofNodesPressure(void){_error_("not implemented yet");};
 		int         NumberofNodesVelocity(void){_error_("not implemented yet");};
-		void        PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");};
-		void			PicoUpdateBox(int loopboxid){_error_("not implemented yet");};
-		void			PicoComputeBasalMelt(void){_error_("not implemented yet");};
 		void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
 		int         PressureInterpolation(void){_error_("not implemented yet");};
Index: issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 24002)
+++ issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 24010)
@@ -88,5 +88,4 @@
 		IssmDouble  IceVolume(bool scaled){_error_("not implemented yet");};
 		IssmDouble  IceVolumeAboveFloatation(bool scaled){_error_("not implemented yet");};
-		void			Ismip6FloatingiceMeltingRate(void){_error_("not implemented yet");};
 		bool        IsFaceOnBoundary(void){_error_("not implemented yet");};
 		bool		   IsIcefront(void);
@@ -138,7 +137,4 @@
 		int         NumberofNodesPressure(void);
 		int         NumberofNodesVelocity(void);
-		void			PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");};	
-		void			PicoUpdateBox(int loopboxid){_error_("not implemented yet");};
-		void			PicoComputeBasalMelt(void){_error_("not implemented yet");};
 		void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
 		int         PressureInterpolation(void);
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24002)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24010)
@@ -2726,60 +2726,4 @@
 	return isicefront;
 }/*}}}*/
-void       Tria::Ismip6FloatingiceMeltingRate(void){/*{{{*/
-
-	if(!this->IsIceInElement() || !this->IsFloating()) return;
-
-	int         basinid,num_basins,M,N;
-	IssmDouble  tf,gamma0,base,delta_t_basin,mean_tf_basin,absval;
-	IssmDouble  basalmeltrate[NUMVERTICES];
-	bool        islocal;
-	IssmDouble* delta_t = NULL;
-	IssmDouble* mean_tf = NULL;
-	IssmDouble* depths  = NULL;
-
-	/*Get variables*/
-	IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum);
-	IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum);
-	IssmDouble lf   = this->FindParam(MaterialsLatentheatEnum);
-	IssmDouble cp   = this->FindParam(MaterialsMixedLayerCapacityEnum);
-
-	/*Hard code sea water density to be consistent with ISMIP6 documentation*/
-	rhow = 1028.;
-
-	/* Get parameters and inputs */
-	this->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum);
-	this->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum);
-	this->parameters->FindParam(&gamma0,BasalforcingsIsmp6Gamma0Enum);
-	this->parameters->FindParam(&delta_t,&M,BasalforcingsIsmp6DeltaTEnum);    _assert_(M==num_basins);
-	this->parameters->FindParam(&islocal,BasalforcingsIsmp6IsLocalEnum);
-	if(!islocal) {
-		this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmp6AverageTfEnum); _assert_(N==num_basins);
-	}
-	Input* tf_input=this->GetInput(BasalforcingsIsmp6TfShelfEnum);            _assert_(tf_input);
-	delta_t_basin = delta_t[basinid]; 
-	if(!islocal) mean_tf_basin = mean_tf[basinid];
-
-	/* Compute melt rate */
-	Gauss* gauss=this->NewGauss();
-	for(int i=0;i<NUMVERTICES;i++){
-		gauss->GaussVertex(i);
-		tf_input->GetInputValue(&tf,gauss);
-		if(!islocal) {
-			absval = mean_tf_basin+delta_t_basin;
-			if (absval<0) {absval = -1.*absval;}
-			basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval;}
-		else {basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*pow(max(tf+delta_t_basin,0.),2);}
-	}
-
-	/*Return basal melt rate*/
-	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrate,P1Enum);
-
-	/*Cleanup and return*/
-	delete gauss;
-	xDelete<IssmDouble>(delta_t);
-	xDelete<IssmDouble>(mean_tf);
-	xDelete<IssmDouble>(depths);
-
-}/*}}}*/
 bool       Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
 
@@ -3269,276 +3213,4 @@
 }
 /*}}}*/
-void       Tria::PicoUpdateBoxid(int* max_boxid_basin_list){/*{{{*/
-
-	if(!this->IsIceInElement() || !this->IsFloating()) return;
-
-	int        basin_id;
-	IssmDouble dist_gl,dist_cf;
-
-	inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum);
-	IssmDouble boxid_max=reCast<IssmDouble>(max_boxid_basin_list[basin_id])+1.;
-
-	Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
-	Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum);  _assert_(dist_cf_input);
-
-	/*Get dist_gl and dist_cf at center of element*/
-	Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
-	dist_gl_input->GetInputValue(&dist_gl,gauss);
-	dist_cf_input->GetInputValue(&dist_cf,gauss);
-	delete gauss;
-
-	/*Ensure values are positive for floating ice*/
-	dist_gl = fabs(dist_gl);
-	dist_cf = fabs(dist_cf);
-
-	/*Compute relative distance to grounding line*/
-	IssmDouble rel_dist_gl=dist_gl/(dist_gl+dist_cf);
-
-	/*Assign box numbers based on rel_dist_gl*/
-	int boxid = -1;
-	for(IssmDouble i=0.;i<boxid_max;i++){
-		IssmDouble lowbound  = 1. -sqrt((boxid_max-i   )/boxid_max);
-		IssmDouble highbound = 1. -sqrt((boxid_max-i-1.)/boxid_max);
-		if(rel_dist_gl>=lowbound && rel_dist_gl<=highbound){
-			boxid=reCast<int>(i);
-			break;
-		}
-	}
-	if(boxid==-1) _error_("No boxid found for element " << this->Sid() << "!");
-
-	this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid));	
-}/*}}}*/
-void       Tria::PicoUpdateBox(int loopboxid){/*{{{*/
-
-	if(!this->IsIceInElement() || !this->IsFloating()) return;
-
-	int boxid;
-	this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
-	if(loopboxid!=boxid) return;
-
-	int        basinid, maxbox, num_basins, numnodes, M;
-	IssmDouble gamma_T, overturning_coeff, thickness;
-	IssmDouble pressure, T_star,p_coeff, q_coeff;
-	bool       isplume;
-
-	/*Get variables*/
-	IssmDouble rhoi       = this->FindParam(MaterialsRhoIceEnum);
-	IssmDouble rhow       = this->FindParam(MaterialsRhoSeawaterEnum);
-	IssmDouble earth_grav = this->FindParam(ConstantsGEnum);
-	IssmDouble rho_star   = 1033.;             // kg/m^3
-	IssmDouble nu         = rhoi/rhow;
-	IssmDouble latentheat = this->FindParam(MaterialsLatentheatEnum);
-	IssmDouble c_p_ocean  = this->FindParam(MaterialsMixedLayerCapacityEnum);
-	IssmDouble lambda     = latentheat/c_p_ocean;
-	IssmDouble a          = -0.0572;          // K/psu
-	IssmDouble b          = 0.0788 + this->FindParam(MaterialsMeltingpointEnum);  //K
-	IssmDouble c          = 7.77e-4;
-	IssmDouble alpha      = 7.5e-5;           // 1/K
-	IssmDouble Beta       = 7.7e-4;           // K
-
-	/* Get non-box-specific parameters and inputs */
-	this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum);
-	this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
-	this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);
-	this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
-	this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
-	this->parameters->FindParam(&isplume, BasalforcingsPicoIsplumeEnum);
-	Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
-	_assert_(basinid<=num_basins);
-
-   IssmDouble* boxareas = NULL;
-	this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
-	_assert_(M==num_basins*maxbox);
-
-	IssmDouble area_boxi        = boxareas[basinid*maxbox+boxid];
-	IssmDouble g1               = area_boxi*gamma_T;
-
-	IssmDouble basalmeltrates_shelf[NUMVERTICES];
-	IssmDouble potential_pressure_melting_point[NUMVERTICES];
-	IssmDouble Tocs[NUMVERTICES];
-	IssmDouble Socs[NUMVERTICES];
-
-	/* First box calculations */
-	if(boxid==0){
-		/* Get box1 parameters and inputs */
-		IssmDouble time, toc_farocean, soc_farocean;
-		this->parameters->FindParam(&time,TimeEnum);
-		this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum);
-		this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);
-		IssmDouble s1 = soc_farocean/(nu*lambda);
-		IssmDouble overturnings[NUMVERTICES];
-
-		/* Start looping on the number of verticies and calculate ocean vars */
-		Gauss* gauss=this->NewGauss();
-		for(int i=0;i<NUMVERTICES;i++){
-			gauss->GaussVertex(i);
-			thickness_input->GetInputValue(&thickness,gauss);
-			pressure = (rhoi*earth_grav*1e-4)*thickness;
-			T_star   = a*soc_farocean+b-c*pressure-toc_farocean;
-			p_coeff  = g1/(overturning_coeff*rho_star*(Beta*s1-alpha));
-			q_coeff  = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha)));
-
-			/* To avoid negatives under the square root */
-			if((0.25*pow(p_coeff,2)-q_coeff)<0) q_coeff = 0.25*p_coeff*p_coeff;
-
-			Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff));
-			Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]);
-			potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
-			if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
-			overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i]));
-		}
-
-		if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
-		this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);
-		this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);
-		this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1Enum);
-
-		/*Cleanup and return*/
-		delete gauss;
-	}
-
-	/* Subsequent box calculations */
-	else {
-		/* Get subsequent box parameters and inputs */
-		IssmDouble* toc_weighted_avg         = NULL;
-		IssmDouble* soc_weighted_avg         = NULL;
-		IssmDouble* overturning_weighted_avg = NULL;
-		this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum);
-		this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum);
-		this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum);
-		IssmDouble mean_toc                  = toc_weighted_avg[basinid];
-		IssmDouble mean_soc                  = soc_weighted_avg[basinid];
-		IssmDouble mean_overturning          = overturning_weighted_avg[basinid];
-		IssmDouble g2                        = g1/(nu*lambda);
-
-		/* Start looping on the number of verticies and calculate ocean vars */
-		Gauss* gauss=this->NewGauss();
-		for(int i=0;i<NUMVERTICES;i++){
-			gauss->GaussVertex(i);
-			thickness_input->GetInputValue(&thickness,gauss);
-			pressure = (rhoi*earth_grav*1.e-4)*thickness;
-			T_star   = a*mean_soc+b-c*pressure-mean_toc;
-			Tocs[i]  = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc));
-			Socs[i]  = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda));
-			potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
-			if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
-		}
-
-		if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
-		this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);
-		this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);
-
-		/*Cleanup and return*/
-		xDelete<IssmDouble>(toc_weighted_avg);
-		xDelete<IssmDouble>(soc_weighted_avg);
-		xDelete<IssmDouble>(overturning_weighted_avg);
-		delete gauss;
-	}
-
-	/*Cleanup and return*/
-	xDelete<IssmDouble>(boxareas);
-}
-/*}}}*/
-void       Tria::PicoComputeBasalMelt(void){/*{{{*/
-
-	if(!this->IsIceInElement() || !this->IsFloating()) return;
-
-   IssmDouble E0, Cd, CdT, YT, lam1, lam2, lam3, M0, CdTS0, y1, y2, x0;
-	IssmDouble Tf_gl, YTS, CdTS, G1, G2, G3, g_alpha, M, l, X_hat, M_hat;
-	IssmDouble alpha, zgl, Toc, Soc, z_base, yts, slopex, slopey;
-
-	/*Get variables*/
-	E0    = 3.6e-2;        //Entrainment coefficient
-	Cd    = 2.5e-3;        //Drag coefficient
-	CdT   = 1.1e-3;        //Turbulent heat exchange coefficient
-	YT    = CdT/sqrt(Cd);  //Heat exchange coefficient
-	lam1  = -5.73e-2;      //Freezing point-salinity coefficient (degrees C)
-	lam2  = 8.32e-2;       //Freezing point offset (degrees C)
-	lam3  = 7.61e-4;       //Freezing point-depth coefficient (K m-1)
-	M0    = 10.;           //Melt-rate parameter (m yr-1 C-2)
-	CdTS0 = 6e-4;          //Heat exchange parameter
-	y1    = 0.545;         //Heat exchange parameter
-	y2    = 3.5e-5;        //Heat exchange parameter
-	x0    = 0.56;          //Dimentionless scaling factor
-
-	/*Define arrays*/
-	IssmDouble basalmeltrates_shelf[NUMVERTICES];  //Basal melt-rate
-
-	/*Polynomial coefficients*/
-	IssmDouble p[12];
-	p[0]  =  0.1371330075095435;
-	p[1]  =  5.527656234709359E1;
-	p[2]  = -8.951812433987858E2;
-	p[3]  =  8.927093637594877E3;
-	p[4]  = -5.563863123811898E4;
-	p[5]  =  2.218596970948727E5;
-	p[6]  = -5.820015295669482E5;
-	p[7]  =  1.015475347943186E6;
-	p[8]  = -1.166290429178556E6;
-	p[9]  =  8.466870335320488E5;
-	p[10] = -3.520598035764990E5;
-	p[11] =  6.387953795485420E4;
-
-	/*Get inputs*/
-   Input* zgl_input            = this->GetInput(GroundinglineHeightEnum);                     _assert_(zgl_input);
-	Input* toc_input            = this->GetInput(BasalforcingsPicoSubShelfOceanTempEnum);      _assert_(toc_input);
-	Input* soc_input            = this->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum);  _assert_(soc_input);
-	Input* base_input           = this->GetInput(BaseEnum);                                    _assert_(base_input);
-	Input* baseslopex_input     = this->GetInput(BaseSlopeXEnum);                              _assert_(baseslopex_input);
-	Input* baseslopey_input     = this->GetInput(BaseSlopeYEnum);                              _assert_(baseslopey_input);
-	this->FindParam(&yts, ConstantsYtsEnum);
-
-	/*Loop over the number of vertices in this element*/
-	Gauss* gauss=this->NewGauss();
-	for(int i=0;i<NUMVERTICES;i++){
-		gauss->GaussVertex(i);
-
-		/*Get inputs*/
-      zgl_input->GetInputValue(&zgl,gauss); 
-		toc_input->GetInputValue(&Toc,gauss); //(K)
-		soc_input->GetInputValue(&Soc,gauss);
-		base_input->GetInputValue(&z_base,gauss);
-		baseslopex_input->GetInputValue(&slopex,gauss);
-		baseslopey_input->GetInputValue(&slopey,gauss);
-
-		/*Compute ice shelf base slope (radians)*/
-	   alpha = atan(sqrt(slopex*slopex + slopey*slopey));
-		if(alpha>=M_PI) alpha = M_PI - 0.001;               //ensure sin(alpha) > 0 for meltrate calculations
-	
-		/*Make necessary conversions*/
-		Toc = Toc-273.15;
-	   if(zgl>z_base) zgl=z_base;
-
-		/*Low bound for Toc to ensure X_hat is between 0 and 1*/
-		if(Toc<lam1*Soc+lam2) Toc=lam1*Soc+lam2;
-
-		/*Compute parameters needed for melt-rate calculation*/
-		Tf_gl = lam1*Soc+lam2+lam3*zgl;                                              //Characteristic freezing point
-		YTS = YT*(y1+y2*(((Toc-Tf_gl)*E0*sin(alpha))/(lam3*(CdTS0+E0*sin(alpha))))); //Effective heat exchange coefficient
-		CdTS = sqrt(Cd)*YTS;                                                         //Heat exchange coefficient
-		G1 = sqrt(sin(alpha)/(Cd+E0*sin(alpha)));                                    //Geometric factor
-		G2 = sqrt(CdTS/(CdTS+E0*sin(alpha)));                                        //Geometric factor
-		G3 = (E0*sin(alpha))/(CdTS+E0*sin(alpha));                                   //Geometric factor
-      g_alpha = G1*G2*G3;                                                          //Melt scaling factor
-		M = M0*g_alpha*pow((Toc-Tf_gl),2);                                           //Melt-rate scale
-		l = ((Toc-Tf_gl)*(x0*CdTS+E0*sin(alpha)))/(lam3*x0*(CdTS+E0*sin(alpha)));    //Length scale
-		X_hat = (z_base-zgl)/l;                                                      //Dimentionless coordinate system
-
-		/*Compute polynomial fit*/
-		M_hat = 0.;                                                                  //Reset summation variable for each node
-		for(int ii=0;ii<12;ii++) {
-		 M_hat += p[ii]*pow(X_hat,ii);                                               //Polynomial fit
-	   }
-
-		/*Compute melt-rate*/
-		basalmeltrates_shelf[i] = (M*M_hat)/yts;                                     //Basal melt-rate (m/s)
-	}
-	  
-	/*Save computed melt-rate*/
-   this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
-
-   /*Cleanup and return*/
-	delete gauss;
-}/*}}}*/
 void       Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
 
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24002)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24010)
@@ -98,5 +98,4 @@
 		IssmDouble  IcefrontMassFluxLevelset(bool scaled);
 		IssmDouble  GroundinglineMassFlux(bool scaled);
-		void        Ismip6FloatingiceMeltingRate(void);
 		void        InputDepthAverageAtBase(int enum_type,int average_enum_type);
 		void        InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
@@ -116,7 +115,4 @@
 		int         NumberofNodesPressure(void);
 		int         NumberofNodesVelocity(void);
-		void        PicoUpdateBoxid(int* pmax_boxid_basin);
-		void        PicoUpdateBox(int loopboxid);
-		void        PicoComputeBasalMelt(void);
 		void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
 		int         PressureInterpolation();
Index: issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp	(revision 24002)
+++ issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp	(revision 24010)
@@ -12,20 +12,20 @@
 	bool isplume;
 
-   /*First, reset all melt to 0 */
-   for(int i=0;i<femmodel->elements->Size();i++){
-	      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-	      int numvertices = element->GetNumberOfVertices();
-	      IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
-	      element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
-	      xDelete<IssmDouble>(values);
-	   }
+	/*First, reset all melt to 0 */
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		int numvertices = element->GetNumberOfVertices();
+		IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
+		element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
+		xDelete<IssmDouble>(values);
+	}
 
 	/*PICO melt rate parameterization (Reese et al., 2018)*/
-   femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
-   UpdateBoxIdsPico(femmodel);
-   ComputeBoxAreasPico(femmodel);
-   for(int i=0;i<maxbox;i++){
-	      UpdateBoxPico(femmodel,i);
-	      ComputeAverageOceanvarsPico(femmodel,i);
+	femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
+	UpdateBoxIdsPico(femmodel);
+	ComputeBoxAreasPico(femmodel);
+	for(int i=0;i<maxbox;i++){
+		UpdateBoxPico(femmodel,i);
+		ComputeAverageOceanvarsPico(femmodel,i);
 	}
 
@@ -103,5 +103,5 @@
 void ComputeBoxAreasPico(FemModel* femmodel){/*{{{*/
 
-	int num_basins,maxbox,basinid,boxid;
+	int num_basins,maxbox,basinid,boxid,domaintype;
 	IssmDouble dist_max,area;
 
@@ -113,8 +113,12 @@
 	for(int i=0;i<femmodel->elements->Size();i++){
 		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		if(!element->IsIceInElement() || !element->IsFloating()) continue;
-		element->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
-		element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
-		boxareas[basinid*maxbox+boxid]+=element->GetHorizontalSurfaceArea();
+		if(!element->IsOnBase()) continue;
+		Element* basalelement = element->SpawnBasalElement();
+		if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
+		basalelement->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
+		basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+		boxareas[basinid*maxbox+boxid]+=basalelement->GetHorizontalSurfaceArea();
+		basalelement->FindParam(&domaintype,DomainTypeEnum);
+		if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	}
 
@@ -140,5 +144,5 @@
 void ComputeAverageOceanvarsPico(FemModel* femmodel, int boxid){/*{{{*/
 
-	int num_basins, basinid, maxbox, M;
+	int num_basins, basinid, maxbox, M, domaintype;
 	IssmDouble area, toc, soc, overturning;
 	IssmDouble* boxareas=NULL;
@@ -157,20 +161,24 @@
 	for(int i=0;i<femmodel->elements->Size();i++){
 		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		if(!element->IsIceInElement() || !element->IsFloating()) continue;
+		if(!element->IsOnBase()) continue;
+		Element* basalelement = element->SpawnBasalElement();
+		if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
 		int el_boxid;
-		element->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
+		basalelement->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
 		if(el_boxid!=boxid) continue;
 
-		Input* tocs_input=element->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(tocs_input); 
-		Input* socs_input=element->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(socs_input);
-
-		element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
-		Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
+		Input* tocs_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(tocs_input); 
+		Input* socs_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(socs_input);
+
+		basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+		Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
 		tocs_input->GetInputValue(&toc,gauss);
 		socs_input->GetInputValue(&soc,gauss);
 		delete gauss;
-		area=element->GetHorizontalSurfaceArea();
+		area=basalelement->GetHorizontalSurfaceArea();
 		toc_weighted_avg[basinid]+=toc*area;
 		soc_weighted_avg[basinid]+=soc*area;
+		basalelement->FindParam(&domaintype,DomainTypeEnum);
+		if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	}
 
@@ -194,17 +202,21 @@
 		for(int i=0;i<femmodel->elements->Size();i++){
 			Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-			if(!element->IsIceInElement() || !element->IsFloating()) continue;
+			if(!element->IsOnBase()) continue;
+			Element* basalelement = element->SpawnBasalElement();
+			if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
 			int el_boxid;
-			element->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
+			basalelement->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
 			if(el_boxid!=boxid) continue;
 
-	     	Input* overturnings_input=element->GetInput(BasalforcingsPicoSubShelfOceanOverturningEnum); _assert_(overturnings_input);
-
-			element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
-			Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
+	     	Input* overturnings_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanOverturningEnum); _assert_(overturnings_input);
+
+			basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+			Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
 			overturnings_input->GetInputValue(&overturning,gauss);
 			delete gauss;
-			area=element->GetHorizontalSurfaceArea();
+			area=basalelement->GetHorizontalSurfaceArea();
 			overturning_weighted_avg[basinid]+=overturning*area;
+			basalelement->FindParam(&domaintype,DomainTypeEnum);
+			if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 		}
 
Index: issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp
===================================================================
--- issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp	(revision 24002)
+++ issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp	(revision 24010)
@@ -73,5 +73,5 @@
 void FloatingiceMeltingRateIsmip6x(FemModel* femmodel){/*{{{*/
 
-	int         num_basins, basinid,num_depths;
+	int         num_basins, basinid, num_depths, domaintype;
 	IssmDouble  area, tf, base, time;
 	bool        islocal;
@@ -154,13 +154,19 @@
 		for(int i=0;i<femmodel->elements->Size();i++){
 			Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-			if(!element->IsIceInElement() || !element->IsFloating()) continue;
-			Input* tf_input=element->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input);
-			element->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum);
-			Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
+			if(!element->IsOnBase()) continue;
+			/*Spawn basal element if on base to compute element area*/
+			Element* basalelement = element->SpawnBasalElement();
+			if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
+			Input* tf_input=basalelement->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input);
+			basalelement->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum);
+			Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
 			tf_input->GetInputValue(&tf,gauss);
 			delete gauss;
-			area=element->GetHorizontalSurfaceArea();
+			area=basalelement->GetHorizontalSurfaceArea();
 			tf_weighted_avg[basinid]+=tf*area;
 			areas_summed[basinid]   +=area;	
+			/*Delete spawned element if we are in 3D*/
+			basalelement->FindParam(&domaintype,DomainTypeEnum);
+			if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 		}
 
@@ -174,4 +180,10 @@
 	}
 
+   /*Compute meltrates*/
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element->Ismip6FloatingiceMeltingRate();
+	}
+
 	/*Cleanup and return */
 	xDelete<IssmDouble>(tf_weighted_avg);
@@ -180,10 +192,4 @@
 	xDelete<IssmDouble>(areas_summed_cpu);
 	xDelete<IssmDouble>(tf_depths);
-
-   /*Compute meltrates*/
-	for(int i=0;i<femmodel->elements->Size();i++){
-		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-		element->Ismip6FloatingiceMeltingRate();
-	}
 }
 /*}}}*/
Index: issm/trunk-jpl/src/m/classes/basalforcingsismip6.m
===================================================================
--- issm/trunk-jpl/src/m/classes/basalforcingsismip6.m	(revision 24002)
+++ issm/trunk-jpl/src/m/classes/basalforcingsismip6.m	(revision 24010)
@@ -19,6 +19,7 @@
 		function self = extrude(self,md) % {{{
 			self.basin_id=project3d(md,'vector',self.basin_id,'type','element','layer',1);
-			self.tf=project3d(md,'vector',self.tf,'type','element','layer',1);
-			self.delta_t=project3d(md,'vector',self.delta_t,'type','element','layer',1);
+			%self.tf=project3d(md,'vector',self.tf,'type','element','layer',1);
+			%self.delta_t=project3d(md,'vector',self.delta_t,'type','element','layer',1);
+			self.tf=project3d(md,'vector',self.tf,'type','node');
 			self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','element','layer',1); %bedrock only gets geothermal flux
 			self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1);
Index: issm/trunk-jpl/src/m/classes/basalforcingspico.m
===================================================================
--- issm/trunk-jpl/src/m/classes/basalforcingspico.m	(revision 24002)
+++ issm/trunk-jpl/src/m/classes/basalforcingspico.m	(revision 24010)
@@ -13,5 +13,5 @@
 		farocean_temperature      = NaN;
 		farocean_salinity         = NaN;
-		isplume                   = NaN;
+		isplume                   = 0;
 		geothermalflux            = NaN;
 		groundedice_melting_rate  = NaN;
@@ -68,5 +68,5 @@
 				md = checkfield(md,'fieldname','basalforcings.overturning_coeff','numel',1,'NaN',1,'Inf',1,'>',0);
 				md = checkfield(md,'fieldname','basalforcings.gamma_T','numel',1,'NaN',1,'Inf',1,'>',0);
-				md = checkfield(md,'fieldname','basalforcings.farocean_temperature','NaN',1,'Inf',1,'>',0,'size',[md.basalforcings.num_basins+1 NaN]);
+				md = checkfield(md,'fieldname','basalforcings.farocean_temperature','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins+1 NaN]);
 				md = checkfield(md,'fieldname','basalforcings.farocean_salinity','NaN',1,'Inf',1,'>',0,'size',[md.basalforcings.num_basins+1 NaN]);
 				md = checkfield(md,'fieldname','basalforcings.isplume','values',[0 1]);
