Index: ../trunk-jpl/src/c/cores/cores.h
===================================================================
--- ../trunk-jpl/src/c/cores/cores.h	(revision 26118)
+++ ../trunk-jpl/src/c/cores/cores.h	(revision 26119)
@@ -62,9 +62,9 @@
 void sealevelchange_geometry(FemModel* femmodel);
 #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);
 void sealevelchange_core_deformation(Vector<IssmDouble>** pN_radial, Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks);
Index: ../trunk-jpl/src/c/cores/sealevelchange_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26118)
+++ ../trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26119)
@@ -16,6 +16,7 @@
 void TransferSealevel(FemModel* femmodel,int forcingenum);
 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);
 /*}}}*/
 
 /*main cores:*/
@@ -389,19 +390,18 @@
 	}
 }; /*}}}*/
 
-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;
 	IssmDouble*            allloads=NULL; 
 	Vector<IssmDouble>*    sealevelloads=NULL;
 	Vector<IssmDouble>*    oldsealevelloads=NULL;
-	IssmDouble             sealevelloadsaverage;
 	IssmDouble*            allsealevelloads=NULL;
 	Vector<IssmDouble>*    oceanareas=NULL;
 	IssmDouble             oceanarea;
@@ -412,13 +412,33 @@
 	IssmDouble           eps_abs;
 	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));
 	barycontrib=barycontribparam->GetParameterValue();
@@ -425,12 +445,17 @@
 	femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
 	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);
 	loads=new Vector<IssmDouble>(nel);
 	sealevelloads=new Vector<IssmDouble>(nel);
 	oceanareas=new Vector<IssmDouble>(nel);
+
+	/*call masks core: */
+	masks=sealevel_masks(femmodel);
 	
 	/*buildup loads: */
 	for(Object* & object : femmodel->elements->objects){
@@ -438,7 +463,7 @@
 		element->SealevelchangeBarystaticLoads(loads, barycontrib,masks); 
 	}
 
-	//Communicate loads from local to global: 
+	//broadcast loads 
 	loads->Assemble(); allloads=loads->ToMPISerial();
 
 	/*convolve loads:*/
@@ -446,17 +471,16 @@
 		Element* element = xDynamicCast<Element*>(object);
 		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: 
+	sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
 	
-	//substract ocean average and barystatic contributionfrom sea level loads: 
-	barystatic=barycontrib->Total()/oceanarea/rho_water;
-	sealevelloads->Shift(-sealevelloadsaverage+barystatic);
+	//broadcast sea level loads 
 	allsealevelloads=sealevelloads->ToMPISerial();
 
 	bool converged=false;
@@ -470,12 +494,12 @@
 			element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,masks);
 		}
 		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();
 
 		//convergence?
@@ -483,10 +507,25 @@
 		if (converged)break;
 	}
 
+	/*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);
 
+
 }
 /*}}}*/