Index: /issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp	(revision 26074)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp	(revision 26075)
@@ -38,6 +38,6 @@
 
 	/*Initialize sea level cumulated sea level loads :*/
-	InputUpdateFromConstantx(inputs,elements,0.,AccumulatedDeltaTwsEnum);
-	InputUpdateFromConstantx(inputs,elements,0.,OldAccumulatedDeltaTwsEnum);
+	iomodel->ConstantToInput(inputs,elements,0.,AccumulatedDeltaTwsEnum,P1Enum);
+	iomodel->ConstantToInput(inputs,elements,0.,OldAccumulatedDeltaTwsEnum,P1Enum);
 	iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",WatercolumnEnum);
 
@@ -79,12 +79,4 @@
 	if(!element->IsOnBase()) return;
 
-	/*deal with logic of accumulating thickness if we are coupled to a 
-	 * sea level core:*/
-	int frequency,count,isgrd;
-	element->FindParam(&isgrd,SolidearthSettingsGRDEnum);
-	if(isgrd){
-		element->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
-		element->FindParam(&count,SealevelchangeRunCountEnum);
-	}	
 	/*Fetch dof list and allocate solution vector*/
 	int *doflist = NULL;
@@ -106,50 +98,4 @@
 	xDelete<IssmDouble>(watercolumn);
 
-	/*Get basal element*/
-	int domaintype; element->FindParam(&domaintype,DomainTypeEnum);
-	Element* basalelement=element;
-	if(domaintype!=Domain2DhorizontalEnum) basalelement = element->SpawnBasalElement();
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numvertices = basalelement->GetNumberOfVertices();
-
-	/*Now, we need to do some "processing"*/
-	watercolumn  = xNew<IssmDouble>(numvertices);
-	IssmDouble* oldwatercolumn      = xNew<IssmDouble>(numvertices);
-	IssmDouble* cumdeltawatercolumn = xNew<IssmDouble>(numvertices);
-	IssmDouble* oldcumdeltawatercolumn = xNew<IssmDouble>(numvertices);
-	IssmDouble* deltawatercolumn = xNew<IssmDouble>(numvertices);
-
-	/*Get previous base, tws, surfac and current sealevel and bed:*/
-	basalelement->GetInputListOnVertices(&watercolumn[0],WatercolumnEnum);
-	basalelement->GetInputListOnVertices(&oldwatercolumn[0],WaterColumnOldEnum);
-	basalelement->GetInputListOnVertices(&cumdeltawatercolumn[0],AccumulatedDeltaTwsEnum);
-
-	/*What is the delta tws forcing the sea-level change core: cumulated over time, hence the +=:*/
-	if(isgrd){
-		for(int i=0;i<numvertices;i++){
-			cumdeltawatercolumn[i] += watercolumn[i]-oldwatercolumn[i];
-		}
-	}
-
-	/*Add input to the element: */
-	if(isgrd){
-		element->AddBasalInput(AccumulatedDeltaTwsEnum,cumdeltawatercolumn,P1Enum);
-		if(count==frequency){
-			basalelement->GetInputListOnVertices(&oldcumdeltawatercolumn[0],OldAccumulatedDeltaTwsEnum);
-			element->AddBasalInput(OldAccumulatedDeltaTwsEnum,cumdeltawatercolumn,P1Enum);
-			for(int i=0;i<numvertices;i++)deltawatercolumn[i]=cumdeltawatercolumn[i]-oldcumdeltawatercolumn[i];
-			element->AddBasalInput(DeltaTwsEnum,deltawatercolumn,P1Enum);
-		}
-	}
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(watercolumn);
-	xDelete<IssmDouble>(deltawatercolumn);
-	xDelete<IssmDouble>(cumdeltawatercolumn);
-	xDelete<IssmDouble>(oldcumdeltawatercolumn);
-	xDelete<IssmDouble>(oldwatercolumn);
-	xDelete<int>(doflist);
-	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 }/*}}}*/
 void           HydrologyTwsAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 26074)
+++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp	(revision 26075)
@@ -9,6 +9,7 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
-
-void hydrology_core(FemModel* femmodel){
+void SolidEarthWaterUpdates(FemModel* femmodel);
+
+void hydrology_core(FemModel* femmodel){ /*{{{*/
 
 	/*Start profiler*/
@@ -43,4 +44,7 @@
 		/*transfer water column thickness to old water column thickness: */
 		InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
+		/*solid earth considerations:*/
+		SolidEarthWaterUpdates(femmodel);
+
 
 	}
@@ -59,4 +63,8 @@
 		GetVectorFromInputsx(&ug,femmodel,HydrologyTwsSpcEnum,VertexPIdEnum);
 		InputUpdateFromSolutionx(femmodel,ug); 
+
+		/*solid earth considerations:*/
+		SolidEarthWaterUpdates(femmodel);
+
 		delete ug;
 	
@@ -219,3 +227,32 @@
 	/*End profiler*/
 	femmodel->profiler->Stop(HYDROLOGYCORE);
-}
+} /*}}}*/
+void SolidEarthWaterUpdates(FemModel* femmodel){ /*{{{*/
+
+	int isgrd;
+	int grdmodel;
+	IssmDouble time;
+	int frequency,count;
+
+	/*retrieve parameters:*/
+	femmodel->parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
+	femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
+	femmodel->parameters->FindParam(&time,TimeEnum);
+	femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
+	femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
+	
+	/*early return?:*/
+	if(!isgrd)return;
+
+	/*From old and new thickness, create delta thickness  and accumulate:*/
+	femmodel->inputs->AXPY(-1, WaterColumnOldEnum,WatercolumnEnum,DeltaTwsEnum);
+	femmodel->inputs->AXPY(+1, DeltaTwsEnum,AccumulatedDeltaTwsEnum,DummyEnum);
+	femmodel->inputs->DuplicateInput(DummyEnum,AccumulatedDeltaTwsEnum);
+
+	/*compute total water column change between two sea-level solver time steps, ie. every frequency*dt:*/
+	if(count==frequency){
+		femmodel->inputs->AXPY(-1, OldAccumulatedDeltaTwsEnum,AccumulatedDeltaTwsEnum,DeltaTwsEnum);
+		femmodel->inputs->DuplicateInput(AccumulatedDeltaTwsEnum,OldAccumulatedDeltaTwsEnum);
+	}
+	return;
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/cores/masstransport_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/masstransport_core.cpp	(revision 26074)
+++ /issm/trunk-jpl/src/c/cores/masstransport_core.cpp	(revision 26075)
@@ -10,6 +10,5 @@
 #include "../solutionsequences/solutionsequences.h"
 #include "../classes/Inputs/TransientInput.h"
-void SolidEarthUpdates(FemModel* femmodel);
-
+void SolidEarthIceUpdates(FemModel* femmodel);
 void masstransport_core(FemModel* femmodel){ /*{{{*/
 
@@ -77,5 +76,5 @@
 			//}
 		}
-		SolidEarthUpdates(femmodel);
+		SolidEarthIceUpdates(femmodel);
 		
 		femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
@@ -100,5 +99,5 @@
 	femmodel->profiler->Stop(MASSTRANSPORTCORE);
 } /*}}}*/
-void SolidEarthUpdates(FemModel* femmodel){ /*{{{*/
+void SolidEarthIceUpdates(FemModel* femmodel){ /*{{{*/
 
 	int isgrd;
Index: /issm/trunk-jpl/test/NightlyRun/test2001.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2001.m	(revision 26074)
+++ /issm/trunk-jpl/test/NightlyRun/test2001.m	(revision 26075)
@@ -50,8 +50,5 @@
 
 %Fields and tolerances to track changes
-field_names     ={'UGia','UGiaRate'};
-field_tolerances={1e-13,1e-13};
-field_values={...
-	(md.results.TransientSolution(2).SealevelUGrd),
-	(md.results.TransientSolution(2).SealevelUGrd)
-	};
+field_names     ={'UGrd'};
+field_tolerances={1e-13};
+field_values={md.results.TransientSolution(2).SealevelUGrd};