Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp	(revision 22468)
+++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp	(revision 22469)
@@ -29,4 +29,8 @@
 			MismipFloatingiceMeltingRatex(femmodel);
 			break;
+		case BasalforcingsPicoEnum:
+			if(VerboseSolution())_printf0_(" call Pico Floating melting rate module\n");
+			PicoFloatingiceMeltingRatex(femmodel);
+			break;
 		default:
 			_error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");
@@ -49,4 +53,89 @@
 		element->MismipFloatingiceMeltingRate();
 	}
+}
+/*}}}*/
 
+void PicoFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/
+
+	int i,k,p,numvertices;
+	int num_basins,maxbox,basinid,boxid;
+	IssmDouble dist_max,area,earth_grav,rhoi,rhow,rho_star,nu,latentheat;
+	IssmDouble c_p_ocean,lambda,a,b,c,alpha,Beta,gamma_T,overturning_coeff;
+	IssmDouble g1,g2,s1,p_coeff,tavg,savg,overturning_avg;
+	IssmDouble mean_toc,mean_soc,mean_overturning,pressure,T_star;
+	IssmDouble q_coeff,potential_pressure_melting_point,T_pressure_melting,overturning;
+	IssmDouble thickness,toc_farocean,soc_farocean,area_box1,Toc,Soc,basalmeltrate_shelf;
+	int* nd=NULL;
+	IssmDouble* dmax_basin=NULL;
+	IssmDouble* t_farocean=NULL;
+	IssmDouble* s_farocean=NULL;
+	IssmDouble* distances=NULL;
+	IssmDouble* boxareas=NULL;
+	Element* element=NULL;
+
+	femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
+	femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
+	dmax_basin=xNew<IssmDouble>(num_basins);
+	femmodel->DistanceToFieldValue(MaskGroundediceLevelsetEnum,0.,DistanceToGroundinglineEnum);
+	femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0.,DistanceToCalvingfrontEnum);
+	// find maximum distance to grounding line
+	dist_max=-1.;
+	for(i=0;i<num_basins;i++){
+		dmax_basin[i]=-1;
+	}
+	for(i=0;i<femmodel->elements->Size();i++){
+		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		if(!element->IsFloating()) continue;
+		numvertices = element->GetNumberOfVertices();
+		distances=xNew<IssmDouble>(numvertices);
+		element->GetInputListOnVertices(&distances[0],DistanceToGroundinglineEnum);
+		element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+		for(k=0; k<numvertices; k++){
+			if(abs(distances[k])>dist_max){dist_max=abs(distances[k]);}
+			if(abs(distances[k])>dmax_basin[basinid]){dmax_basin[basinid]=abs(distances[k]);}
+		}
+	}
+	//1 Define pico boxes
+	nd=xNew<int>(num_basins);
+	for(i=0; i<num_basins;i++){
+		nd[i] = round(sqrt(dmax_basin[i]/dist_max)*(maxbox-1)); //0-based
+	}
+	for(int i=0;i<femmodel->elements->Size();i++){
+		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		element->PicoUpdateBoxid(nd);
+	}
+
+	//2 Get area of the boxes
+	boxareas=new Vector<IssmDouble>(num_basins*maxbox);
+	for(i=0;i<num_basins*maxbox;i++){boxareas[i]=0.;}
+	for(i=0;i<femmodel->elements->Size();i++){
+		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+		if(!element->IsFloating()) continue;
+		element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+		element->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
+		boxareas[basinid*maxbox+boxid]+=element->GetArea();
+	}
+
+	//3 Compute variables in first box
+	for(int p=0;p<femmodel->elements->Size();p++){
+		element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(p));
+		element->PicoUpdateFirstBox();
+		//Need to save the box1 average of temp, sal, and overt for input into next box.
+	}
+
+
+	//4 Compute variables of all other boxes
+
+// 	for(int i=0;i<femmodel->elements->Size();i++){
+// 		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+// 		element->PicoFloatingiceMeltingRate();
+// 	}
+	/*Cleanup and return */
+	xDelete<int>(nd);
+	xDelete<IssmDouble>(dmax_basin);
+	xDelete<IssmDouble>(distances);
+	xDelete<IssmDouble>(boxareas);
+	xDelete<IssmDouble>(t_farocean);
+   xDelete<IssmDouble>(s_farocean);
+	xDelete<IssmDouble>(element);
 }/*}}}*/