Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 23651)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 23652)
@@ -227,4 +227,5 @@
 					./modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp\
 					./modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp\
+					./modules/FrontalForcingsx/FrontalForcingsx.cpp\
 					./modules/ConfigureObjectsx/ConfigureObjectsx.cpp\
 					./modules/SpcNodesx/SpcNodesx.cpp\
Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 23651)
+++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp	(revision 23652)
@@ -58,24 +58,18 @@
 		case DefaultCalvingEnum:
 			iomodel->FetchDataToInput(elements,"md.calving.calvingrate",CalvingCalvingrateEnum);
-			iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum);
 			break;
 		case CalvingLevermannEnum:
 			iomodel->FetchDataToInput(elements,"md.calving.coeff",CalvinglevermannCoeffEnum);
-			iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvinglevermannMeltingrateEnum);
 			break;
 		case CalvingVonmisesEnum:
 			iomodel->FetchDataToInput(elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum);
 			iomodel->FetchDataToInput(elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum);
-			iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum);
 			break;
 		case CalvingMinthicknessEnum:
-			iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum);
 			break;
 		case CalvingHabEnum:
-			iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum);
 			iomodel->FetchDataToInput(elements,"md.calving.flotation_fraction",CalvingHabFractionEnum);
 			break;
 		case CalvingCrevasseDepthEnum:
-			iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum);
 			iomodel->FetchDataToInput(elements,"md.calving.water_height",WaterheightEnum);
 			break;
@@ -83,8 +77,23 @@
 			iomodel->FetchDataToInput(elements,"md.calving.stress_threshold_groundedice",CalvingStressThresholdGroundediceEnum);
 			iomodel->FetchDataToInput(elements,"md.calving.stress_threshold_floatingice",CalvingStressThresholdFloatingiceEnum);
-			iomodel->FetchDataToInput(elements,"md.calving.meltingrate",CalvingMeltingrateEnum);
 			break;
 		default:
 			_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
+	}
+
+	/*Get frontal melt parameters*/
+	int melt_parameterization;
+	iomodel->FindConstant(&melt_parameterization,"md.frontalforcings.melt_parameterization");
+	switch(melt_parameterization){
+		case 0:
+			iomodel->FetchDataToInput(elements,"md.frontalforcings.meltingrate",CalvingMeltingrateEnum);
+			break;
+		case 1:
+			iomodel->FetchDataToInput(elements,"md.frontalforcings.basin",FrontalForcingsBasinIdEnum);
+			iomodel->FetchDataToInput(elements,"md.frontalforcings.subglacial_discharge",FrontalForcingsSubglacialDischargeEnum);
+			iomodel->FetchDataToInput(elements,"md.frontalforcings.thermalforcing",FrontalForcingsThermalForcingEnum);
+		break;
+		default:
+			_error_("Frontal forcing model not supported yet");
 	}
 }
@@ -94,4 +103,6 @@
 	parameters->AddObject(iomodel->CopyConstantObject("md.levelset.reinit_frequency",LevelsetReinitFrequencyEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.levelset.calving_max",CalvingMaxEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.melt_parameterization",FrontalForcingsParamEnum));
+	
 	int  calvinglaw;
 	iomodel->FindConstant(&calvinglaw,"md.calving.law");
@@ -114,4 +125,11 @@
 		default:
 			_error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
+	}
+
+	/*Get frontal melt parameters*/
+	int melt_parameterization;
+	iomodel->FindConstant(&melt_parameterization,"md.frontalforcings.melt_parameterization");
+	if(melt_parameterization==1){
+		parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.numberofbasins",FrontalForcingsNumberofBasinsEnum));
 	}
 	return;
@@ -250,5 +268,5 @@
 				default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
 			}
-			meltingrate_input = basalelement->GetInput(CalvinglevermannMeltingrateEnum);     _assert_(meltingrate_input);
+			meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
 			break;
 		case CalvingMinthicknessEnum:
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 23651)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 23652)
@@ -223,4 +223,5 @@
 		virtual void       GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0;
 		virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0;
+		virtual IssmDouble GetIcefrontArea(){_error_("not implemented");};
 		virtual void       GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0;
 		virtual void       GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0;
@@ -300,4 +301,5 @@
 		virtual void       ResetFSBasalBoundaryCondition()=0;
 		virtual void       ResetHooks()=0;
+		virtual void       RignotMeltParameterization(void){_error_("not implemented yet");};
 		virtual void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M)=0;
 		virtual void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 23651)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 23652)
@@ -1397,4 +1397,103 @@
 
 	return phi;
+}
+/*}}}*/
+IssmDouble Tria::GetIcefrontArea(){/*{{{*/
+	
+	IssmDouble  bed[NUMVERTICES]; //basinId[NUMVERTICES];
+	IssmDouble	Haverage,frontarea;
+	IssmDouble  x1,y1,x2,y2,distance;
+	IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
+	int* indices=NULL;
+	IssmDouble* H=NULL;;
+	int nrfrontbed,numiceverts;
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&bed[0],BedEnum);
+	GetInputListOnVertices(&surfaces[0],SurfaceEnum);
+	GetInputListOnVertices(&bases[0],BaseEnum);
+	GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+
+	if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
+
+	nrfrontbed=0;
+	for(int i=0;i<NUMVERTICES;i++){
+		/*Find if bed<0*/
+		if(bed[i]<0.) nrfrontbed++;
+	}
+
+	if(nrfrontbed==3){
+		/*2. Find coordinates of where levelset crosses 0*/
+		int         numiceverts;
+		IssmDouble  s[2],x[2],y[2];
+		int        *indices = NULL;
+		this->GetLevelsetIntersection(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
+		_assert_(numiceverts); 
+
+		/*3 Write coordinates*/
+		IssmDouble  xyz_list[NUMVERTICES][3];
+		::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+		int counter = 0;
+		if((numiceverts>0) && (numiceverts<NUMVERTICES)){
+			for(int i=0;i<numiceverts;i++){
+				for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
+					x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
+					y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
+					counter++;
+				}
+			}
+		}
+		else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
+
+			for(int i=0;i<NUMVERTICES;i++){
+				if(lsf[indices[i]]==0.){
+					x[counter]=xyz_list[indices[i]][0];
+					y[counter]=xyz_list[indices[i]][1];
+					counter++;
+				}
+				if(counter==2) break;
+			}
+			if(counter==1){
+				/*We actually have only 1 vertex on levelset, write a single point as a segment*/
+				x[counter]=x[0];
+				y[counter]=y[0];
+				counter++;
+			}
+		}
+		else{
+			_error_("not sure what's going on here...");
+		}
+		x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
+		distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
+	}
+	else return 0;	
+	
+	IssmDouble s[2]; // s:fraction of intersected triangle edges, that lies inside ice
+	this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);	
+	int numthk=numiceverts+2;
+	H=xNew<IssmDouble>(numthk);
+	for(int iv=0;iv<NUMVERTICES;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
+
+	switch(numiceverts){
+		case 1: // average over triangle
+			H[0]=Haux[0];
+			H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
+			H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
+			Haverage=(H[1]+H[2])/2;
+			break;
+		case 2: // average over quadrangle
+			H[0]=Haux[0];
+			H[1]=Haux[1];
+			H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
+			H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
+			Haverage=(H[2]+H[3])/2;
+			break;
+		default:
+			_error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
+			break;
+	}
+	frontarea=distance*Haverage;
+	_assert_(frontarea>0);
+	return frontarea;
 }
 /*}}}*/
@@ -3262,4 +3361,60 @@
 }
 /*}}}*/
+void       Tria::RignotMeltParameterization(){/*{{{*/
+
+   IssmDouble A, B, alpha, beta;
+	IssmDouble bed,qsg,qsg_basin,TF,yts;
+	int numbasins;
+	IssmDouble basinid[NUMVERTICES];
+	IssmDouble* basin_icefront_area=NULL;
+
+	/* Coefficients */
+	A    = 3e-4;        //
+	B    = 0.15;        //
+	alpha = 0.39;
+	beta = 1.18;
+	
+	/*Get inputs*/
+	Input* bed_input = this->GetInput(BedEnum);                     _assert_(bed_input);
+	Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum);		 _assert_(qsg_input);
+	Input* TF_input  = this->GetInput(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
+	GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);	
+	
+	this->FindParam(&yts, ConstantsYtsEnum);
+	this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
+	this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum);
+
+	IssmDouble meltrates[NUMVERTICES];  //frontal melt-rate
+	
+	/* Start looping on the number of vertices: */
+	GaussTria* gauss=new GaussTria();
+	for(int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get variables */
+		bed_input->GetInputValue(&bed,gauss);
+		qsg_input->GetInputValue(&qsg,gauss);
+		TF_input->GetInputValue(&TF,gauss);
+
+		if(basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
+		else{
+			/* change the unit of qsg (m^3/d -> m/d) with ice front area */
+			qsg_basin=qsg/basin_icefront_area[reCast<int>(basinid[iv])-1];
+
+			/* calculate melt rates */
+			meltrates[iv]=(A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta)/86400; //[m/s]
+		}	
+
+		if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
+		if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new TriaInput(CalvingMeltingrateEnum,&meltrates[0],P1Enum));
+   
+	/*Cleanup and return*/
+	delete gauss;
+}
+/*}}}*/
 void       Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N, int M){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 23651)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 23652)
@@ -79,4 +79,5 @@
 		void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
 		IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
+		IssmDouble  GetIcefrontArea();
 		void	      GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
 		void	      GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
@@ -120,4 +121,5 @@
 		void        ResetHooks();
 		void        ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);
+		void        RignotMeltParameterization();
 		void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N,int M);
 		void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23651)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 23652)
@@ -1496,4 +1496,36 @@
 
 }/*}}}*/
+void FemModel::IcefrontAreax(){/*{{{*/
+
+	int numvertices      = this->GetElementsWidth();
+	int numbasins;
+	IssmDouble* BasinId   = xNew<IssmDouble>(numvertices);
+	this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
+	IssmDouble* basin_icefront_area           = xNewZeroInit<IssmDouble>(numbasins);
+
+	for(int basin=1;basin<numbasins+1;basin++){
+		IssmDouble local_icefront_area = 0;
+		IssmDouble total_icefront_area;
+
+		for(int i=0;i<this->elements->Size();i++){
+			Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
+			element->GetInputListOnVertices(BasinId,FrontalForcingsBasinIdEnum);
+			for(int j=0;j<numvertices;j++){
+				if(BasinId[j]==basin){
+					local_icefront_area+=element->GetIcefrontArea();
+					break;
+				}
+			}
+		}
+		ISSM_MPI_Reduce(&local_icefront_area,&total_icefront_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
+		ISSM_MPI_Bcast(&total_icefront_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	
+		basin_icefront_area[basin-1]=total_icefront_area;
+	}
+	
+	this->parameters->AddObject(new DoubleVecParam(FrontalForcingsBasinIcefrontAreaEnum,basin_icefront_area,numbasins));
+	
+	xDelete<IssmDouble>(basin_icefront_area);
+}/*}}}*/
 void FemModel::IceMassx(IssmDouble* pM, bool scaled){/*{{{*/
 
@@ -2390,4 +2422,12 @@
 }
 /*}}}*/
+void FemModel::RignotMeltParameterizationx(){/*{{{*/
+
+	for(int i=0;i<elements->Size();i++){
+		Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->RignotMeltParameterization();
+	}
+}
+/*}}}*/
 void FemModel::StrainRateparallelx(){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 23651)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 23652)
@@ -100,4 +100,5 @@
 		void GetLocalVectorWithClonesNodes(IssmDouble** plocal_vector,Vector<IssmDouble> *vector);
 		void GroundedAreax(IssmDouble* pV, bool scaled);
+		void IcefrontAreax();
 		void IceMassx(IssmDouble* pV, bool scaled);
 		void IceVolumex(IssmDouble* pV, bool scaled);
@@ -122,4 +123,5 @@
 		void StrainRateeffectivex();
 		void StressIntensityFactorx();
+		void RignotMeltParameterizationx();
 		void TotalFloatingBmbx(IssmDouble* pFbmb, bool scaled);
 		void TotalGroundedBmbx(IssmDouble* pGbmb, bool scaled);
Index: /issm/trunk-jpl/src/c/cores/movingfront_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/movingfront_core.cpp	(revision 23651)
+++ /issm/trunk-jpl/src/c/cores/movingfront_core.cpp	(revision 23652)
@@ -36,4 +36,5 @@
 	/* start the work from here */
 	Calvingx(femmodel);
+	FrontalForcingsx(femmodel);
 	if(VerboseSolution()) _printf0_("   computing level set transport\n");
 
Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp	(revision 23652)
+++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp	(revision 23652)
@@ -0,0 +1,26 @@
+/*!\file FrontalForcingsx
+ * \brief: compute ice frontal melting rate
+ */
+
+#include "./FrontalForcingsx.h"
+#include "../../shared/shared.h"
+#include "../../toolkits/toolkits.h"
+
+void FrontalForcingsx(FemModel* femmodel){
+
+	/*Recover melt_parameterization*/
+	int melt_parameterization;
+	femmodel->parameters->FindParam(&melt_parameterization,FrontalForcingsParamEnum);
+	
+	/*Calculate melting rate*/
+	switch(melt_parameterization){
+		case 0:
+			break;
+		case 1:
+			femmodel->IcefrontAreax();
+			femmodel->RignotMeltParameterizationx();
+			break;
+		default:
+			_error_("Frontal forcings parameterization not supported yet");
+	}
+}
Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h	(revision 23652)
+++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.h	(revision 23652)
@@ -0,0 +1,10 @@
+#ifndef _FRONTALFORCINGSX_H
+#define _FRONTALFORCINGSX_H
+
+#include "../../classes/classes.h"
+#include "../../analyses/analyses.h"
+
+/* local prototypes: */
+void FrontalForcingsx(FemModel* femmodel);
+
+#endif
Index: /issm/trunk-jpl/src/c/modules/modules.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/modules.h	(revision 23651)
+++ /issm/trunk-jpl/src/c/modules/modules.h	(revision 23652)
@@ -38,4 +38,5 @@
 #include "./Gradjx/Gradjx.h"
 #include "./GroundinglineMigrationx/GroundinglineMigrationx.h"
+#include "./FrontalForcingsx/FrontalForcingsx.h"
 #include "./InputDepthAverageAtBasex/InputDepthAverageAtBasex.h"
 #include "./InputDuplicatex/InputDuplicatex.h"
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 23651)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 23652)
@@ -134,4 +134,8 @@
 	FrictionDeltaEnum,
 	FrictionVoidRatioEnum,
+	FrontalForcingsBasinIcefrontAreaEnum,
+	FrontalForcingsBasinIdEnum,
+	FrontalForcingsNumberofBasinsEnum,
+	FrontalForcingsParamEnum,
 	GiaCrossSectionShapeEnum,
 	GroundinglineMigrationEnum,
@@ -444,5 +448,4 @@
 	CalvingHabFractionEnum,
 	CalvinglevermannCoeffEnum,
-	CalvinglevermannMeltingrateEnum,
 	CalvingMeltingrateEnum,
 	CalvingratexAverageEnum,
@@ -495,4 +498,6 @@
 	FrictionQEnum,
 	FrictionWaterLayerEnum,
+	FrontalForcingsSubglacialDischargeEnum,
+	FrontalForcingsThermalForcingEnum,
 	HydrologyWatercolumnMaxEnum,
 	FrictionTillFrictionAngleEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 23651)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 23652)
@@ -142,4 +142,8 @@
 		case FrictionDeltaEnum : return "FrictionDelta";
 		case FrictionVoidRatioEnum : return "FrictionVoidRatio";
+		case FrontalForcingsBasinIcefrontAreaEnum : return "FrontalForcingsBasinIcefrontArea";
+		case FrontalForcingsBasinIdEnum : return "FrontalForcingsBasinId";
+		case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins";
+		case FrontalForcingsParamEnum : return "FrontalForcingsParam";
 		case GiaCrossSectionShapeEnum : return "GiaCrossSectionShape";
 		case GroundinglineMigrationEnum : return "GroundinglineMigration";
@@ -450,5 +454,4 @@
 		case CalvingHabFractionEnum : return "CalvingHabFraction";
 		case CalvinglevermannCoeffEnum : return "CalvinglevermannCoeff";
-		case CalvinglevermannMeltingrateEnum : return "CalvinglevermannMeltingrate";
 		case CalvingMeltingrateEnum : return "CalvingMeltingrate";
 		case CalvingratexAverageEnum : return "CalvingratexAverage";
@@ -501,4 +504,6 @@
 		case FrictionQEnum : return "FrictionQ";
 		case FrictionWaterLayerEnum : return "FrictionWaterLayer";
+		case FrontalForcingsSubglacialDischargeEnum : return "FrontalForcingsSubglacialDischarge";
+		case FrontalForcingsThermalForcingEnum : return "FrontalForcingsThermalForcing";
 		case HydrologyWatercolumnMaxEnum : return "HydrologyWatercolumnMax";
 		case FrictionTillFrictionAngleEnum : return "FrictionTillFrictionAngle";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 23651)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 23652)
@@ -145,4 +145,8 @@
 	      else if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum;
 	      else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum;
+	      else if (strcmp(name,"FrontalForcingsBasinIcefrontArea")==0) return FrontalForcingsBasinIcefrontAreaEnum;
+	      else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum;
+	      else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum;
+	      else if (strcmp(name,"FrontalForcingsParam")==0) return FrontalForcingsParamEnum;
 	      else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum;
 	      else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
@@ -256,12 +260,12 @@
 	      else if (strcmp(name,"ModelId")==0) return ModelIdEnum;
 	      else if (strcmp(name,"Nodes")==0) return NodesEnum;
-	      else if (strcmp(name,"NumModels")==0) return NumModelsEnum;
+         else stage=3;
+   }
+   if(stage==3){
+	      if (strcmp(name,"NumModels")==0) return NumModelsEnum;
 	      else if (strcmp(name,"OceanGridNx")==0) return OceanGridNxEnum;
 	      else if (strcmp(name,"OceanGridNy")==0) return OceanGridNyEnum;
 	      else if (strcmp(name,"OceanGridX")==0) return OceanGridXEnum;
-         else stage=3;
-   }
-   if(stage==3){
-	      if (strcmp(name,"OceanGridY")==0) return OceanGridYEnum;
+	      else if (strcmp(name,"OceanGridY")==0) return OceanGridYEnum;
 	      else if (strcmp(name,"OutputBufferPointer")==0) return OutputBufferPointerEnum;
 	      else if (strcmp(name,"OutputBufferSizePointer")==0) return OutputBufferSizePointerEnum;
@@ -379,12 +383,12 @@
 	      else if (strcmp(name,"StressbalanceNumRequestedOutputs")==0) return StressbalanceNumRequestedOutputsEnum;
 	      else if (strcmp(name,"StressbalancePenaltyFactor")==0) return StressbalancePenaltyFactorEnum;
-	      else if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum;
+         else stage=4;
+   }
+   if(stage==4){
+	      if (strcmp(name,"StressbalanceReltol")==0) return StressbalanceReltolEnum;
 	      else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum;
 	      else if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum;
 	      else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
-         else stage=4;
-   }
-   if(stage==4){
-	      if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
+	      else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
 	      else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
 	      else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
@@ -459,5 +463,4 @@
 	      else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
 	      else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
-	      else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
 	      else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
 	      else if (strcmp(name,"CalvingratexAverage")==0) return CalvingratexAverageEnum;
@@ -503,14 +506,16 @@
 	      else if (strcmp(name,"FrictionC")==0) return FrictionCEnum;
 	      else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum;
-	      else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
-	      else if (strcmp(name,"FrictionEffectivePressure")==0) return FrictionEffectivePressureEnum;
-	      else if (strcmp(name,"FrictionM")==0) return FrictionMEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
+	      if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
+	      else if (strcmp(name,"FrictionEffectivePressure")==0) return FrictionEffectivePressureEnum;
+	      else if (strcmp(name,"FrictionM")==0) return FrictionMEnum;
+	      else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
 	      else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
 	      else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
 	      else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
+	      else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum;
+	      else if (strcmp(name,"FrontalForcingsThermalForcing")==0) return FrontalForcingsThermalForcingEnum;
 	      else if (strcmp(name,"HydrologyWatercolumnMax")==0) return HydrologyWatercolumnMaxEnum;
 	      else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
@@ -624,13 +629,13 @@
 	      else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum;
 	      else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum;
-	      else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
 	      else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
 	      else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
 	      else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
 	      else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
-         else stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
+	      else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
 	      else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
 	      else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
@@ -747,13 +752,13 @@
 	      else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
 	      else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
-	      else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
+         else stage=7;
+   }
+   if(stage==7){
+	      if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
 	      else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum;
 	      else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum;
 	      else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
 	      else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
-         else stage=7;
-   }
-   if(stage==7){
-	      if (strcmp(name,"AutodiffJacobian")==0) return AutodiffJacobianEnum;
+	      else if (strcmp(name,"AutodiffJacobian")==0) return AutodiffJacobianEnum;
 	      else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
 	      else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
@@ -870,13 +875,13 @@
 	      else if (strcmp(name,"Gset")==0) return GsetEnum;
 	      else if (strcmp(name,"Gsl")==0) return GslEnum;
-	      else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
+         else stage=8;
+   }
+   if(stage==8){
+	      if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
 	      else if (strcmp(name,"HOFSApproximation")==0) return HOFSApproximationEnum;
 	      else if (strcmp(name,"Hook")==0) return HookEnum;
 	      else if (strcmp(name,"HydrologyBasalFlux")==0) return HydrologyBasalFluxEnum;
 	      else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
-         else stage=8;
-   }
-   if(stage==8){
-	      if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
+	      else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
 	      else if (strcmp(name,"HydrologydcEplThicknessStacked")==0) return HydrologydcEplThicknessStackedEnum;
 	      else if (strcmp(name,"HydrologydcEplThickness")==0) return HydrologydcEplThicknessEnum;
@@ -993,13 +998,13 @@
 	      else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
 	      else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
-	      else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
+         else stage=9;
+   }
+   if(stage==9){
+	      if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
 	      else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
 	      else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
 	      else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
 	      else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
-         else stage=9;
-   }
-   if(stage==9){
-	      if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
+	      else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
 	      else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
 	      else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
@@ -1116,13 +1121,13 @@
 	      else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
 	      else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
-	      else if (strcmp(name,"SealevelInertiaTensorXZ")==0) return SealevelInertiaTensorXZEnum;
+         else stage=10;
+   }
+   if(stage==10){
+	      if (strcmp(name,"SealevelInertiaTensorXZ")==0) return SealevelInertiaTensorXZEnum;
 	      else if (strcmp(name,"SealevelInertiaTensorYZ")==0) return SealevelInertiaTensorYZEnum;
 	      else if (strcmp(name,"SealevelInertiaTensorZZ")==0) return SealevelInertiaTensorZZEnum;
 	      else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum;
 	      else if (strcmp(name,"SealevelriseAnalysis")==0) return SealevelriseAnalysisEnum;
-         else stage=10;
-   }
-   if(stage==10){
-	      if (strcmp(name,"SealevelriseSolution")==0) return SealevelriseSolutionEnum;
+	      else if (strcmp(name,"SealevelriseSolution")==0) return SealevelriseSolutionEnum;
 	      else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum;
 	      else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum;
Index: /issm/trunk-jpl/src/m/classes/calving.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/calving.m	(revision 23651)
+++ /issm/trunk-jpl/src/m/classes/calving.m	(revision 23652)
@@ -30,5 +30,4 @@
 		function self = extrude(self,md) % {{{
 			self.calvingrate=project3d(md,'vector',self.calvingrate,'type','node');
-			self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
 		end % }}}
 		function self = setdefaultparameters(self) % {{{
@@ -40,10 +39,8 @@
 
 			md = checkfield(md,'fieldname','calving.calvingrate','>=',0,'timeseries',1,'NaN',1,'Inf',1);
-			md = checkfield(md,'fieldname','calving.meltingrate','>=',0,'timeseries',1,'NaN',1,'Inf',1);
 		end % }}}
 		function disp(self) % {{{
 			disp(sprintf('   Calving parameters:'));
 			fielddisplay(self,'calvingrate','calving rate at given location [m/a]');
-			fielddisplay(self,'meltingrate','melting rate at given location [m/a]');
 		end % }}}
 		function marshall(self,prefix,md,fid) % {{{
@@ -51,5 +48,4 @@
 			WriteData(fid,prefix,'name','md.calving.law','data',1,'format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','calvingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts);
-			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts);
 		end % }}}
 		function savemodeljs(self,fid,modelname) % {{{
Index: /issm/trunk-jpl/src/m/classes/calvingcrevassedepth.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/calvingcrevassedepth.m	(revision 23651)
+++ /issm/trunk-jpl/src/m/classes/calvingcrevassedepth.m	(revision 23652)
@@ -30,5 +30,4 @@
 		end % }}}
 		function self = extrude(self,md) % {{{
-			self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
 		end % }}}
 		function self = setdefaultparameters(self) % {{{
@@ -43,5 +42,4 @@
 			md = checkfield(md,'fieldname','calving.crevasse_opening_stress','numel',[1],'values',[0,1]);
 			md = checkfield(md,'fieldname','calving.water_height','NaN',1,'Inf',1,'timeseries',1,'>=',0);
-			md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0);
 		end % }}}
 		function disp(self) % {{{
@@ -49,5 +47,4 @@
 			fielddisplay(self,'crevasse_opening_stress','0: stress only in the ice-flow direction, 1: max principal');
 			fielddisplay(self,'water_height','water height in the crevasse [m]');
-			fielddisplay(self,'meltingrate','melting rate at given location [m/a]');
 
 		end % }}}
@@ -57,5 +54,4 @@
 			WriteData(fid,prefix,'object',self,'fieldname','crevasse_opening_stress','format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','water_height','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
-			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts);
 		end % }}}
 	end
Index: /issm/trunk-jpl/src/m/classes/calvinglevermann.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/calvinglevermann.m	(revision 23651)
+++ /issm/trunk-jpl/src/m/classes/calvinglevermann.m	(revision 23652)
@@ -30,5 +30,4 @@
 		function self = extrude(self,md) % {{{
 			self.coeff=project3d(md,'vector',self.coeff,'type','node');
-			self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
 		end % }}}
 		function self = setdefaultparameters(self) % {{{
@@ -42,10 +41,8 @@
 
 			md = checkfield(md,'fieldname','calving.coeff','>',0,'size',[md.mesh.numberofvertices 1]);
-			md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0);
 		end % }}}
 		function disp(self) % {{{
 			disp(sprintf('   Calving Levermann parameters:'));
 			fielddisplay(self,'coeff','proportionality coefficient in Levermann model');
-			fielddisplay(self,'meltingrate','melting rate at given location [m/a]');
 
 		end % }}}
@@ -54,5 +51,4 @@
 			WriteData(fid,prefix,'name','md.calving.law','data',3,'format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','coeff','format','DoubleMat','mattype',1);
-			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts);
 		end % }}}
 	end
Index: /issm/trunk-jpl/src/m/classes/calvingminthickness.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/calvingminthickness.m	(revision 23651)
+++ /issm/trunk-jpl/src/m/classes/calvingminthickness.m	(revision 23652)
@@ -29,5 +29,4 @@
 		end % }}}
 		function self = extrude(self,md) % {{{
-			self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
 		end % }}}
 		function self = setdefaultparameters(self) % {{{
@@ -41,10 +40,8 @@
 
 			md = checkfield(md,'fieldname','calving.min_thickness','>',0,'NaN',1,'Inf',1);
-			md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>=',0);
 		end % }}}
 		function disp(self) % {{{
 			disp(sprintf('   Calving Minimum thickness:'));
 			fielddisplay(self,'min_thickness','minimum thickness below which no ice is allowed');
-			fielddisplay(self,'meltingrate','melting rate at given location [m/a]');
 
 		end % }}}
@@ -53,5 +50,4 @@
 			WriteData(fid,prefix,'name','md.calving.law','data',4,'format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','min_thickness','format','Double');
-			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'scale',1./yts);
 		end % }}}
 	end
Index: /issm/trunk-jpl/src/m/classes/calvingvonmises.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/calvingvonmises.m	(revision 23651)
+++ /issm/trunk-jpl/src/m/classes/calvingvonmises.m	(revision 23652)
@@ -8,5 +8,5 @@
 		stress_threshold_groundedice = 0.;
 		stress_threshold_floatingice = 0.;
-		meltingrate   = NaN;
+		meltingrate=NaN;
 	end
 	methods
@@ -30,5 +30,4 @@
 		end % }}}
 		function self = extrude(self,md) % {{{
-			self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
 		end % }}}
 		function self = setdefaultparameters(self) % {{{
@@ -44,5 +43,4 @@
 			md = checkfield(md,'fieldname','calving.stress_threshold_groundedice','>',0,'nan',1,'Inf',1);
 			md = checkfield(md,'fieldname','calving.stress_threshold_floatingice','>',0,'nan',1,'Inf',1);
-			md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0);
 		end % }}}
 		function disp(self) % {{{
@@ -50,5 +48,4 @@
 			fielddisplay(self,'stress_threshold_groundedice','sigma_max applied to grounded ice only [Pa]');
 			fielddisplay(self,'stress_threshold_floatingice','sigma_max applied to floating ice only [Pa]');
-			fielddisplay(self,'meltingrate','melting rate at given location [m/a]');
 
 		end % }}}
@@ -58,5 +55,4 @@
 			WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_groundedice','format','DoubleMat','mattype',1);
 			WriteData(fid,prefix,'object',self,'fieldname','stress_threshold_floatingice','format','DoubleMat','mattype',1);
-			WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts);
 		end % }}}
 	end
Index: /issm/trunk-jpl/src/m/classes/frontalforcings.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/frontalforcings.m	(revision 23652)
+++ /issm/trunk-jpl/src/m/classes/frontalforcings.m	(revision 23652)
@@ -0,0 +1,87 @@
+%FRONTAL FORCINGS class definition
+%
+%   Usage:
+%      frontalforcings=frontalforcings();
+
+classdef frontalforcings
+	properties (SetAccess=public) 
+		meltingrate   = NaN;
+		melt_parameterization =0;
+		basin= NaN;
+		numberofbasins=0;
+		subglacial_discharge= NaN;
+		thermalforcing=NaN;
+	end
+	methods
+		function self = frontalforcings(varargin) % {{{
+			switch nargin
+				case 0
+					self=setdefaultparameters(self);
+				case 1
+					inputstruct=varargin{1};
+					list1 = properties('frontalforcings');
+					list2 = fieldnames(inputstruct);
+					for i=1:length(list1)
+						fieldname = list1{i};
+						if ismember(fieldname,list2),
+							self.(fieldname) = inputstruct.(fieldname);
+						end
+					end
+				otherwise
+					error('constructor not supported');
+			end
+		end % }}}
+		function self = extrude(self,md) % {{{
+			self.meltingrate=project3d(md,'vector',self.meltingrate,'type','node');
+		end % }}}
+		function self = setdefaultparameters(self) % {{{
+
+			meltingrate   = NaN;
+			melt_parameterization =0;
+			basin=NaN;
+			numberofbasins=0;
+			subglacial_discharge=NaN;
+			thermalforcing=NaN;
+		end % }}}
+		function md = checkconsistency(self,md,solution,analyses) % {{{
+			%Early return
+			if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end
+			
+			md = checkfield(md,'fieldname','frontalforcings.melt_parameterization','numel',[1],'values',[0,1]);
+			if self.melt_parameterization==0,
+				md = checkfield(md,'fieldname','frontalforcings.meltingrate','NaN',1,'Inf',1,'timeseries',1,'>=',0);
+			elseif self.melt_parameterization==1,
+				md = checkfield(md,'fieldname','frontalforcings.basin','>',0,'nan',1,'Inf',1);
+				md = checkfield(md,'fieldname','frontalforcings.numberofbasins','numel',[1]);
+				md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'nan',1,'Inf',1,'timeseries',1);
+				md = checkfield(md,'fieldname','frontalforcings.thermalforcing','nan',1,'Inf',1,'timeseries',1);
+			end
+
+		end % }}}
+		function disp(self) % {{{
+			disp(sprintf('   Frontalforcings parameters:'));
+			fielddisplay(self,'melt_parameterization','0: no, 1: Rignot melt parameterization (Rignot et al.,2016)');
+			if self.melt_parameterization==0,
+				fielddisplay(self,'meltingrate','melting rate at given location [m/a]');
+			elseif self.melt_parameterization==1,
+				fielddisplay(self,'basin','basin ID for vertices');
+				fielddisplay(self,'numberofbasins', 'number of basins');
+				fielddisplay(self,'subglacial_discharge','sum of subglacial discharge for each basin [m/d]');
+				fielddisplay(self,'thermalforcing','thermal forcing [∘C]');
+			end
+		end % }}}
+		function marshall(self,prefix,md,fid) % {{{
+			yts=md.constants.yts;
+			WriteData(fid,prefix,'object',self,'fieldname','melt_parameterization','format','Integer');
+			switch self.melt_parameterization
+				case 0
+					WriteData(fid,prefix,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1./yts);
+				case 1
+					WriteData(fid,prefix,'object',self,'fieldname','basin','format','DoubleMat','mattype',1);
+					WriteData(fid,prefix,'object',self,'fieldname','numberofbasins','format','Integer');
+					WriteData(fid,prefix,'object',self,'fieldname','subglacial_discharge','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
+					WriteData(fid,prefix,'object',self,'fieldname','thermalforcing','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
+			end
+		end % }}}
+	end
+end
Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 23651)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 23652)
@@ -40,4 +40,5 @@
 		levelset			  = 0;
 		calving          = 0;
+		frontalforcings  = 0;
 		love			     = 0;
 		gia				  = 0;
@@ -153,5 +154,11 @@
 				md.settings.sb_coupling_frequency=1;
 			end
-
+			%2019 Jan..
+			if isa(md.frontalforcings,'double');
+				if(~isnan(md.calving.meltingrate))
+					disp('Warning: md.calving.meltingrate is now in md.frontalforcings');
+				end
+				md.frontalforcings=frontalforcings(md.calving); 
+			end
 		end% }}}
 	end
@@ -307,6 +314,6 @@
 				md.calving.coeff=project2d(md,md.calving.coeff,1); 
 			end
-			if isprop(md.calving,'meltingrate') & ~isnan(md.calving.meltingrate),
-				md.calving.meltingrate=project2d(md,md.calving.meltingrate,1); 
+			if isprop(md.frontalforcings,'meltingrate') & ~isnan(md.frontalforcings.meltingrate),
+				md.frontalforcings.meltingrate=project2d(md,md.frontalforcings.meltingrate,1); 
 			end
 
@@ -840,4 +847,5 @@
 			md.levelset=extrude(md.levelset,md);
 			md.calving=extrude(md.calving,md);
+			md.frontalforcings=extrude(md.frontalforcings,md);
 			md.hydrology = extrude(md.hydrology,md);
 			md.slr = extrude(md.slr,md);
@@ -1148,5 +1156,4 @@
 				md.smb=SMBhenning(structmd.surfaceforcings);
 			end
-
 		end% }}}
 		function md = setdefaultparameters(md) % {{{
@@ -1182,4 +1189,5 @@
 			md.levelset         = levelset();
 			md.calving          = calving();
+			md.frontalforcings  = frontalforcings();
 			md.gia				  = giaivins();
 			md.esa              = esa();
@@ -1357,4 +1365,5 @@
 			disp(sprintf('%19s: %-22s -- %s','levelset'        ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)'));
 			disp(sprintf('%19s: %-22s -- %s','calving'         ,['[1x1 ' class(self.calving) ']'],'parameters for calving'));
+			disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings'));
 			disp(sprintf('%19s: %-22s -- %s','gia'             ,['[1x1 ' class(self.gia) ']'],'parameters for gia solution'));
 			disp(sprintf('%19s: %-22s -- %s','esa'             ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
Index: /issm/trunk-jpl/test/NightlyRun/test540.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test540.m	(revision 23651)
+++ /issm/trunk-jpl/test/NightlyRun/test540.m	(revision 23652)
@@ -10,5 +10,5 @@
 md.mask.ice_levelset = 1e4*(md.mask.ice_levelset + 0.5);
 md.calving=calvingvonmises();
-md.calving.meltingrate = zeros(md.mesh.numberofvertices,1);
+md.frontalforcings.meltingrate = zeros(md.mesh.numberofvertices,1);
 md.transient.ismovingfront = 1;
 md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1);
Index: /issm/trunk-jpl/test/NightlyRun/test804.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test804.m	(revision 23651)
+++ /issm/trunk-jpl/test/NightlyRun/test804.m	(revision 23652)
@@ -16,5 +16,5 @@
 
 md.calving.calvingrate=1000.*ones(md.mesh.numberofvertices,1);
-md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
+md.frontalforcings.meltingrate=zeros(md.mesh.numberofvertices,1);
 
 md=solve(md,'Transient');
Index: /issm/trunk-jpl/test/NightlyRun/test805.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test805.m	(revision 23651)
+++ /issm/trunk-jpl/test/NightlyRun/test805.m	(revision 23652)
@@ -23,5 +23,5 @@
 
 md.calving.calvingrate=1000.*ones(md.mesh.numberofvertices,1);
-md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
+md.frontalforcings.meltingrate=zeros(md.mesh.numberofvertices,1);
 md.groundingline.melt_interpolation='SubelementMelt1';
 md.levelset.stabilization=2;
Index: /issm/trunk-jpl/test/NightlyRun/test806.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test806.m	(revision 23651)
+++ /issm/trunk-jpl/test/NightlyRun/test806.m	(revision 23652)
@@ -27,5 +27,5 @@
 md.calving=calvinglevermann();
 md.calving.coeff=4.89e13*ones(md.mesh.numberofvertices,1);
-md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
+md.frontalforcings.meltingrate=zeros(md.mesh.numberofvertices,1);
 md.levelset.spclevelset=NaN(md.mesh.numberofvertices,1);
 
Index: /issm/trunk-jpl/test/NightlyRun/test807.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test807.m	(revision 23651)
+++ /issm/trunk-jpl/test/NightlyRun/test807.m	(revision 23652)
@@ -26,5 +26,5 @@
 
 md.calving.calvingrate=zeros(md.mesh.numberofvertices,1);
-md.calving.meltingrate=10000*ones(md.mesh.numberofvertices,1);
+md.frontalforcings.meltingrate=10000*ones(md.mesh.numberofvertices,1);
 md.levelset.spclevelset=NaN(md.mesh.numberofvertices,1);
 
Index: /issm/trunk-jpl/test/NightlyRun/test808.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test808.m	(revision 23651)
+++ /issm/trunk-jpl/test/NightlyRun/test808.m	(revision 23652)
@@ -27,5 +27,5 @@
 md.calving=calvingminthickness();
 md.calving.min_thickness=400;
-md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
+md.frontalforcings.meltingrate=zeros(md.mesh.numberofvertices,1);
 md.levelset.spclevelset=NaN(md.mesh.numberofvertices,1);
 md.levelset.reinit_frequency=1;
Index: /issm/trunk-jpl/test/NightlyRun/test809.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test809.m	(revision 23651)
+++ /issm/trunk-jpl/test/NightlyRun/test809.m	(revision 23652)
@@ -20,5 +20,5 @@
 md.calving.crevasse_opening_stress=1;
 md.calving.water_height=50*ones(md.mesh.numberofvertices,1);
-md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
+md.frontalforcings.meltingrate=zeros(md.mesh.numberofvertices,1);
 md.levelset.spclevelset=NaN(md.mesh.numberofvertices,1);
 md.levelset.reinit_frequency=1;
Index: /issm/trunk-jpl/test/Par/ValleyGlacierShelf.par
===================================================================
--- /issm/trunk-jpl/test/Par/ValleyGlacierShelf.par	(revision 23651)
+++ /issm/trunk-jpl/test/Par/ValleyGlacierShelf.par	(revision 23652)
@@ -75,5 +75,5 @@
 %Masstransport;
 md.calving.calvingrate = 0.*ones(md.mesh.numberofvertices,1);
-md.calving.meltingrate = 0.*ones(md.mesh.numberofvertices,1);
+md.frontalforcings.meltingrate = 0.*ones(md.mesh.numberofvertices,1);
 md.levelset.spclevelset=NaN(md.mesh.numberofvertices,1);
 md.masstransport.stabilization = 1.;