Index: /issm/trunk-jpl/src/c/cores/cores.h
===================================================================
--- /issm/trunk-jpl/src/c/cores/cores.h	(revision 26118)
+++ /issm/trunk-jpl/src/c/cores/cores.h	(revision 26119)
@@ -63,7 +63,7 @@
 #endif
 void grd_core(FemModel* femmodel);
+void grd_core_optim(FemModel* femmodel);
 void solidearthexternal_core(FemModel* femmodel);
 void dynstr_core(FemModel* femmodel);
-SealevelMasks* sealevel_masks(FemModel* femmodel);
 Vector<IssmDouble>* sealevelchange_core_barystatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea);
 Vector<IssmDouble>* sealevelchange_core_sal(FemModel* femmodel,SealevelMasks* masks, Vector<IssmDouble>* RSLg_barystatic,IssmDouble oceanarea);
Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26118)
+++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26119)
@@ -17,4 +17,5 @@
 void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
 IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea);
+SealevelMasks* sealevel_masks(FemModel* femmodel);
 /*}}}*/
 
@@ -390,11 +391,11 @@
 }; /*}}}*/
 
-void grd_core_optim(FemModel* femmodel,SealevelMasks* masks) { /*{{{*/
+void grd_core_optim(FemModel* femmodel) { /*{{{*/
 
 	/*variables:{{{*/
 	int nel;
 	BarystaticContributions* barycontrib=NULL;
+	SealevelMasks* masks=NULL;
 	GenericParam<BarystaticContributions*>* barycontribparam=NULL;
-	IssmDouble               barystatic;
 	
 	Vector<IssmDouble>*    loads=NULL;
@@ -402,5 +403,4 @@
 	Vector<IssmDouble>*    sealevelloads=NULL;
 	Vector<IssmDouble>*    oldsealevelloads=NULL;
-	IssmDouble             sealevelloadsaverage;
 	IssmDouble*            allsealevelloads=NULL;
 	Vector<IssmDouble>*    oceanareas=NULL;
@@ -413,11 +413,31 @@
 	int                  step;
 	IssmDouble           time; 
-	
-	IssmDouble cumbslc;
-	IssmDouble cumbslcice;
-	IssmDouble cumbslchydro;
+
+	int  modelid,earthid;
+	int  horiz;
+	int  count,frequency,iscoupling;
+	int  grd=0;
 	/*}}}*/
 
-	/*retrieve parameters: */
+	/*Verbose: */
+	if(VerboseSolution()) _printf0_("	  computing GRD patterns\n");
+
+	/*retrieve parameters:{{{*/
+	femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); 
+	femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
+	femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
+	/*}}}*/
+
+	/*only run if grd was requested, if we are the earth, and we have reached
+	 * the necessary number of time steps dictated by :*/
+	if(!grd)            return;
+	if(count!=frequency)return;
+	femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum);
+	if(iscoupling){
+		femmodel->parameters->FindParam(&modelid,ModelIdEnum);
+		femmodel->parameters->FindParam(&earthid,EarthIdEnum);
+		if(modelid!=earthid)return;
+	}
+	/*retrieve parameters: {{{*/ 
 	femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
 	barycontribparam = xDynamicCast<GenericParam<BarystaticContributions*>*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum));
@@ -426,5 +446,7 @@
 	femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
 	femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
-	
+	femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
+	/*}}}*/
+
 	/*initialize matrices and vectors:*/
 	femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
@@ -432,4 +454,7 @@
 	sealevelloads=new Vector<IssmDouble>(nel);
 	oceanareas=new Vector<IssmDouble>(nel);
+
+	/*call masks core: */
+	masks=sealevel_masks(femmodel);
 	
 	/*buildup loads: */
@@ -439,5 +464,5 @@
 	}
 
-	//Communicate loads from local to global: 
+	//broadcast loads 
 	loads->Assemble(); allloads=loads->ToMPISerial();
 
@@ -447,15 +472,14 @@
 		element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,masks);
 	}
-
+	sealevelloads->Assemble();
+	
 	//Get ocean area: 
 	oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
 	if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
 
-	//Get sea level loads ocean average: 
-	sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
-	
 	//substract ocean average and barystatic contributionfrom sea level loads: 
-	barystatic=barycontrib->Total()/oceanarea/rho_water;
-	sealevelloads->Shift(-sealevelloadsaverage+barystatic);
+	sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
+	
+	//broadcast sea level loads 
 	allsealevelloads=sealevelloads->ToMPISerial();
 
@@ -471,10 +495,10 @@
 		}
 		sealevelloads->Assemble();
+	
 		
 		//substract ocean average and barystatic contribution
-		sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
-		sealevelloads->Shift(-sealevelloadsaverage+barystatic);
-
-		//broadcast loads 
+		sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
+		
+		//broadcast sea level loads 
 		allsealevelloads=sealevelloads->ToMPISerial();
 
@@ -484,7 +508,22 @@
 	}
 
+	/*convolve loads and sea level loads to get the deformation:*/
+	for(Object* & object : femmodel->elements->objects){
+		Element* element = xDynamicCast<Element*>(object);
+		//element->SealevelchangeDeformationConvolution(allsealevelloads, allloads,masks);
+	}
+
+	/*Update bedrock motion and geoid:*/
+	/*inputs->AXPY(SealevelEnum,SealevelGRD,1);
+	inputs->AXPY(BedEnum,BedGRD,1);
+	if(horiz){
+		inputs->AXPY(BedrockeastEnum,BedEastGRD,1);
+		inputs->AXPY(BedrocknorthEnum,BedNorthGRD,1);
+	}*/
+
 	/*cumulate barystatic contributions and save to results: */
 	barycontrib->Cumulate(femmodel->parameters);
 	barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
+
 
 }