Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp	(revision 22161)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp	(revision 22162)
@@ -50,4 +50,6 @@
 	iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
 	iomodel->FetchDataToInput(elements,"md.geometry.bed",BedEnum);
+	iomodel->FetchDataToInput(elements,"md.slr.Ngia",SealevelNGiaRateEnum);
+	iomodel->FetchDataToInput(elements,"md.slr.Ugia",SealevelUGiaRateEnum);
 
 	/*Initialize cumdeltalthickness and sealevel rise rate input: unfortunately, we don't have femmodel, so we 
@@ -59,7 +61,5 @@
 			element->InputUpdateFromConstant(0.0,SealevelriseCumDeltathicknessEnum);
 			element->InputUpdateFromConstant(0.0,SealevelNEsaRateEnum);
-			element->InputUpdateFromConstant(0.0,SealevelNGiaRateEnum);
 			element->InputUpdateFromConstant(0.0,SealevelUEsaRateEnum);
-			element->InputUpdateFromConstant(0.0,SealevelUGiaRateEnum);
 			counter++;
 		}
Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp	(revision 22161)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp	(revision 22162)
@@ -527,19 +527,23 @@
 void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia, Vector<IssmDouble>** pN_gia,FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/
 
+	/*variables:*/
 	Vector<IssmDouble> *U_gia  = NULL; 
 	Vector<IssmDouble> *N_gia  = NULL; 
-		
-	int configuration_type;
-	int  gsize;
+	
+	/*parameters:*/
+	int					frequency;
+	IssmDouble          dt;
 
 	/*retrieve some parameters:*/
-	femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-
-	/*find size of vectors:*/
-	gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
-
-	/*intialize vectors:*/
-	U_gia = new Vector<IssmDouble>(gsize); U_gia->Set(0);
-	N_gia = new Vector<IssmDouble>(gsize); N_gia->Set(0);
+	femmodel->parameters->FindParam(&frequency,SealevelriseGeodeticRunFrequencyEnum);
+	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
+	/*recover GIA rates:*/
+	GetVectorFromInputsx(&U_gia,femmodel,SealevelUGiaRateEnum,VertexSIdEnum);
+	GetVectorFromInputsx(&N_gia,femmodel,SealevelNGiaRateEnum,VertexSIdEnum);
+
+	/*we just loaded rates, that's not what's being asked, scale by time:*/ 
+	U_gia->Scale(frequency*dt);
+	N_gia->Scale(frequency*dt);
 
 	/*Assign output pointers:*/
Index: /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/slr.m
===================================================================
--- /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/slr.m	(revision 22161)
+++ /issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/slr.m	(revision 22162)
@@ -31,4 +31,6 @@
 		loop_increment = 0;
 		horiz = 0;
+		Ngia  = 0;
+		Ugia  = 0;
 		requested_outputs      = {};
 		transitions    = {};
@@ -118,4 +120,6 @@
 			md = checkfield(md,'fieldname','slr.loop_increment','NaN',1,'Inf',1,'>=',1);
 			md = checkfield(md,'fieldname','slr.horiz','NaN',1,'Inf',1,'values',[0 1]);
+			md = checkfield(md,'fieldname','slr.Ngia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
+			md = checkfield(md,'fieldname','slr.Ugia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
 
 			%check that love numbers are provided at the same level of accuracy: 
@@ -171,4 +175,6 @@
 			fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); 
 			fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)'); 
+			fielddisplay(self,'Ngia','rate of viscous (GIA) geoid expansion (in mm/yr)'); 
+			fielddisplay(self,'Ugia','rate of viscous (GIA) bedrock uplift (in mm/yr)'); 
 			fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)');
 			fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
@@ -203,4 +209,6 @@
 			WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
+			WriteData(fid,prefix,'object',self,'fieldname','Ngia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
+			WriteData(fid,prefix,'object',self,'fieldname','Ugia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
 			WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
 			WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');