Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15868)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15869)
@@ -5059,5 +5059,109 @@
 void Penta::ComputeBasalMeltrate(void){
   /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
-  // Do nothing for now.
+
+  /*Parameters*/
+  int dim;
+  parameters->FindParam(&dim, MeshDimensionEnum);
+
+  /* Intermediaries */
+  bool isenthalpy;
+  int i, analysis_type;
+  IssmDouble xyz_list[NUMVERTICES][3];
+  IssmDouble heatflux, geothermalflux_value;
+  IssmDouble vec_heatflux[dim];
+  IssmDouble normal_base[dim], d1enthalpy[dim];
+  IssmDouble kappa;
+  IssmDouble enthalpy, enthalpyup;
+  IssmDouble pressure, pressureup;
+  IssmDouble meltrate, watercolumn;
+  IssmDouble temperature, waterfraction;
+  IssmDouble latentheat;
+  IssmDouble basalfriction, alpha2;
+  IssmDouble vx,vy,vz;
+  IssmDouble dt;
+  Friction*  friction=NULL;
+  GaussPenta* gauss=NULL;
+  GaussPenta* gaussup=NULL;
+
+  /*check that dimension is 3*/
+  if(dim!=3) return;
+  /* Only compute melt rates at the base of grounded ice*/
+  if(!IsOnBed() || IsFloating()) return;
+  /*Check wether enthalpy is activated*/
+  parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
+  if(!isenthalpy) return;
+
+  /*Fetch parameters: */
+  latentheat=matpar->GetLatentHeat();
+
+  Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
+  Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
+  Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
+  Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
+  Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
+  Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
+  Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
+
+  /*Build friction element, needed later: */
+  parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+  friction=new Friction("3d",inputs,matpar,analysis_type);
+
+  /* Start  looping on the number of gaussian points: */
+  gauss=new GaussPenta(0,1,2,2);
+  gaussup=new GaussPenta(3,4,5,2);
+  for(int ig=gauss->begin();ig<gauss->end();ig++){
+    gauss->GaussPoint(ig);
+    gaussup->GaussPoint(ig);
+    
+    watercolumn_input->GetInputValue(&watercolumn, gauss);
+    enthalpy_input->GetInputValue(&enthalpy, gauss);
+    pressure_input->GetInputValue(&pressure, gauss);
+
+    /*Calculate basal meltrate*/
+    if((watercolumn>0) && (enthalpy<matpar->PureIceEnthalpy(pressure)))
+      enthalpy=matpar->PureIceEnthalpy(pressure);
+    else if(enthalpy<matpar->PureIceEnthalpy(pressure)){
+      meltrate=0.;   //TODO: set zero meltrate and watercolumn in model
+      watercolumn=0.;
+      return;
+    }
+
+    /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice) */
+    enthalpy_input->GetInputValue(&enthalpyup, gaussup);
+    pressure_input->GetInputValue(&pressureup, gaussup);    
+    if(enthalpyup>=matpar->PureIceEnthalpy(pressureup))
+      for(i=0;i<dim;i++)
+	vec_heatflux[i]=0.;
+    else{
+      enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
+      kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
+      for(i=0;i<dim;i++)
+	vec_heatflux[i]=-1.*kappa*d1enthalpy[i];
+    }
+    /*Get normal at base*/
+    normal_base[0]=0.; normal_base[1]=0.; normal_base[2]=1.; // TODO: get normal from model geometry
+    
+    heatflux=0.;
+    for(i=0;i<dim;i++)
+      heatflux+=(vec_heatflux[i])*normal_base[i];
+    geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
+    
+    matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy,pressure);
+    
+    friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
+    vx_input->GetInputValue(&vx,gauss);
+    vy_input->GetInputValue(&vy,gauss);
+    vz_input->GetInputValue(&vz,gauss);
+    basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
+    meltrate=(basalfriction-(heatflux-geothermalflux_value))/(1-waterfraction)/latentheat; 
+    
+    /*Update water column*/
+    this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+    if(reCast<bool,IssmDouble>(dt))
+      watercolumn+=dt*meltrate;
+    else
+      watercolumn=meltrate;
+    // TODO: feed meltrate & watercolumn back to model
+  }  
 }
 /*}}}*/