Index: /issm/trunk-jpl/src/c/Makefile.am
===================================================================
--- /issm/trunk-jpl/src/c/Makefile.am	(revision 26109)
+++ /issm/trunk-jpl/src/c/Makefile.am	(revision 26110)
@@ -97,4 +97,5 @@
 	./classes/Hook.cpp \
 	./classes/Radar.cpp \
+	./classes/BarystaticContributions.cpp \
 	./classes/ExternalResults/Results.cpp \
 	./classes/Elements/Element.cpp \
Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp	(revision 26109)
+++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp	(revision 26110)
@@ -183,5 +183,7 @@
 		xDelete<IssmDouble>(partitionocean);
 	}
-
+	/*New optimized code:*/
+	BarystaticContributions* barystaticcontributions=new BarystaticContributions(iomodel);
+	parameters->AddObject(new GenericParam<BarystaticContributions*>(barystaticcontributions,BarystaticContributionsEnum));
 	
 	/*Deal with external multi-model ensembles: {{{*/
Index: /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp	(revision 26110)
+++ /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp	(revision 26110)
@@ -0,0 +1,149 @@
+/*!\file BarystaticContributions.c
+ * \brief: implementation of the BarystaticContributions object
+ */
+
+/*Include files: {{{*/
+#ifdef HAVE_CONFIG_H
+	#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./BarystaticContributions.h" 
+#include "../toolkits/toolkits.h"
+#include "./classes.h"
+/*}}}*/
+
+/*Constructors and destructors:*/
+BarystaticContributions::BarystaticContributions(IoModel* iomodel ){ /*{{{*/
+
+	int nel;
+
+	iomodel->FetchData(&nice,"md.solidearth.npartice");
+	if(nice){
+		iomodel->FetchData(&pice,&nel,NULL,"md.solidearth.partitionice");
+		ice=new Vector<IssmDouble>(nice);
+		cumice=new Vector<IssmDouble>(nice); cumice->Set(0); cumice->Assemble();
+	}
+
+	iomodel->FetchData(&nhydro,"md.solidearth.nparthydro");
+	if(nhydro){
+		iomodel->FetchData(&phydro,&nel,NULL,"md.solidearth.partitionhydro");
+		hydro=new Vector<IssmDouble>(nhydro);
+		cumhydro=new Vector<IssmDouble>(nhydro); cumhydro->Set(0); cumhydro->Assemble();
+	}
+	iomodel->FetchData(&nocean,"md.solidearth.npartocean");
+	if(nocean){
+		iomodel->FetchData(&pocean,&nel,NULL,"md.solidearth.partitionocean");
+		ocean=new Vector<IssmDouble>(nocean);
+		cumocean=new Vector<IssmDouble>(nocean); cumocean->Set(0); cumocean->Assemble();
+	}
+
+} /*}}}*/
+BarystaticContributions::~BarystaticContributions(){ /*{{{*/
+	delete ice;   delete cumice;
+	delete hydro; delete cumhydro;
+	delete ocean; delete cumocean;
+	if(nice)xDelete<IssmDouble>(pice);
+	if(nhydro)xDelete<IssmDouble>(phydro);
+	if(nocean)xDelete<IssmDouble>(pocean);
+}; /*}}}*/
+
+/*Support routines:*/
+IssmDouble BarystaticContributions::Total(){ /*{{{*/
+
+	IssmDouble  sumice,sumhydro,sumocean;
+	
+	ice->Assemble();
+	hydro->Assemble();
+	ocean->Assemble();
+
+	ice->Sum(&sumice);
+	hydro->Sum(&sumhydro);
+	ocean->Sum(&sumocean);
+
+	return sumice+sumhydro+sumocean;
+
+} /*}}}*/
+IssmDouble BarystaticContributions::CumTotal(){ /*{{{*/
+
+	IssmDouble sumice,sumhydro,sumocean;
+
+	cumice->Assemble();
+	cumhydro->Assemble();
+	cumocean->Assemble();
+
+	cumice->Sum(&sumice);
+	cumhydro->Sum(&sumhydro);
+	cumocean->Sum(&sumocean);
+
+
+	return sumice+sumhydro+sumocean;
+
+} /*}}}*/
+void BarystaticContributions::Cumulate(Parameters* parameters){ /*{{{*/
+
+	cumice->AXPY(ice,1);
+	cumocean->AXPY(ocean,1);
+	cumhydro->AXPY(hydro,1);
+
+
+} /*}}}*/
+void BarystaticContributions::Save(Results* results, Parameters* parameters, IssmDouble oceanarea){ /*{{{*/
+
+	int        step;
+	IssmDouble time;
+	IssmDouble rho_water;
+
+	IssmDouble* cumice_serial=NULL;
+	IssmDouble* cumhydro_serial=NULL;
+	IssmDouble* cumocean_serial=NULL;
+
+	IssmDouble sumice,sumhydro,sumocean;
+
+	parameters->FindParam(&step,StepEnum);
+	parameters->FindParam(&time,TimeEnum);
+	parameters->FindParam(&rho_water,TimeEnum);
+
+	ice->Sum(&sumice); hydro->Sum(&sumhydro); ocean->Sum(&sumocean);
+	results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcEnum,-this->Total()/oceanarea/rho_water,step,time));
+	results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcIceEnum,-sumice/oceanarea/rho_water,step,time));
+	results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcHydroEnum,-sumice/oceanarea/rho_water,step,time));
+	results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcOceanEnum,-sumocean/oceanarea/rho_water,step,time));
+
+	cumice->Sum(&sumice); cumhydro->Sum(&sumhydro); cumocean->Sum(&sumocean);
+	results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcEnum,this->CumTotal()/oceanarea/rho_water,step,time));
+	results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcIceEnum,sumice/oceanarea/rho_water,step,time));
+	results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcHydroEnum,sumhydro/oceanarea/rho_water,step,time));
+	results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcOceanEnum,sumocean/oceanarea/rho_water,step,time));
+
+	cumice_serial=this->cumice->ToMPISerial0(); for (int i=0;i<nice;i++)cumice_serial[i]=cumice_serial[i]/oceanarea/rho_water;
+	cumhydro_serial=this->cumhydro->ToMPISerial0(); for (int i=0;i<nhydro;i++)cumhydro_serial[i]=cumhydro_serial[i]/oceanarea/rho_water;
+	cumocean_serial=this->cumocean->ToMPISerial0(); for (int i=0;i<nocean;i++)cumocean_serial[i]=cumocean_serial[i]/oceanarea/rho_water;
+	
+	results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcIcePartitionEnum,cumice_serial,nice,1,step,time));
+	results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcHydroPartitionEnum,cumhydro_serial,nhydro,1,step,time));
+	results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcOceanPartitionEnum,cumocean_serial,nocean,1,step,time));
+
+	if(IssmComm::GetRank()==0){
+		xDelete<IssmDouble>(cumice_serial);
+		xDelete<IssmDouble>(cumhydro_serial);
+		xDelete<IssmDouble>(cumocean_serial);
+	}
+	return;
+
+} /*}}}*/
+void BarystaticContributions::Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue){ /*{{{*/
+	
+	int id;
+		
+	id=reCast<int>(pice[eid]);
+	ice->SetValue(id,icevalue,ADD_VAL);
+
+	id=reCast<int>(phydro[eid]);
+	hydro->SetValue(id,hydrovalue,ADD_VAL);
+
+	id=reCast<int>(pocean[eid]);
+	ocean->SetValue(id,oceanvalue,ADD_VAL);
+
+} /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/BarystaticContributions.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/BarystaticContributions.h	(revision 26110)
+++ /issm/trunk-jpl/src/c/classes/BarystaticContributions.h	(revision 26110)
@@ -0,0 +1,46 @@
+/*!\file BarystaticContributions.h
+ * \brief: header file for barystatic contribution object
+ */
+
+#ifndef _BARYSTATICCONTRIBUTIONS_H_
+#define _BARYSTATICCONTRIBUTIONS_H_
+
+/*Headers:*/
+class IoModel;
+class Parameters;
+class Results;
+template <class doubletype> class Vector;
+#include "../shared/shared.h"
+
+class BarystaticContributions {
+
+	public: 
+
+		Vector<IssmDouble>* ice;  //contributions to every ice partition
+		Vector<IssmDouble>* cumice;  //cumulated contributions to every ice partition
+		int                 nice; //number of ice partitions 
+		IssmDouble*         pice; //ice partition
+
+		Vector<IssmDouble>* hydro;  //contributions to every hydro partition
+		Vector<IssmDouble>* cumhydro;  //cumulated contributions to every hydro partition
+		int                 nhydro; //number of hydro partitions 
+		IssmDouble*         phydro; //hydro partition
+
+		Vector<IssmDouble>* ocean;  //contributions to every ocean partition
+		Vector<IssmDouble>* cumocean;  //cumulated contributions to every ocean partition
+		int                 nocean; //number of ocean partitions 
+		IssmDouble*         pocean; //ocean partition
+
+		/*BarystaticContributions constructors, destructors :*/
+		BarystaticContributions(IoModel* iomodel );
+		~BarystaticContributions();
+
+		/*routines:*/
+		IssmDouble Total();
+		IssmDouble CumTotal();
+		void Cumulate(Parameters* parameters);
+		void Save(Results* results, Parameters* parameters, IssmDouble oceanarea);
+		void Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue);
+
+};
+#endif  /* _BARYSTATICCONTRIBUTIONS_H_ */
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26110)
@@ -38,4 +38,5 @@
 class ElementMatrix;
 class ElementVector;
+class BarystaticContributions;
 /*}}}*/
 
@@ -389,4 +390,10 @@
 		virtual void          DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0;
 		virtual void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0;
+
+		virtual void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0;
+		virtual void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0;
+		virtual void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks)=0;
+		virtual void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks)=0;
+		virtual void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0;
 		#endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 26110)
@@ -227,4 +227,10 @@
 		void    DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
 		void           GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
+
+		void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
+		void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
+		void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+		void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+		void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
 		#endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 26110)
@@ -182,4 +182,11 @@
 		IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
 		void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
+
+		void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
+		void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
+		void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+		void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+		void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
+
 #endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 26110)
@@ -189,4 +189,10 @@
 		IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
 		void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
+		void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
+		void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
+		void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+		void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+		void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
+
 #endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26109)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26110)
@@ -5432,4 +5432,5 @@
 #endif
 #ifdef _HAVE_SEALEVELCHANGE_
+//old code
 IssmDouble Tria::OceanAverage(IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
 
@@ -5451,5 +5452,5 @@
 }
 /*}}}*/
-void	Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
+void	      Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
 	/*early return if we are not on an ice cap OR ocean:*/
 	if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){
@@ -5545,5 +5546,5 @@
 	return;
 }/*}}}*/
-void    Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
+void       Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
 
 	masks->isiceonly[this->lid]=this->IsIceOnlyInElement();
@@ -5558,8 +5559,7 @@
 	if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true;
 	else masks->notfullygrounded[this->lid]=false;
-
-}
-/*}}}*/
-void    Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
+}
+/*}}}*/
+void       Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
 
 	/*diverse:*/
@@ -5951,5 +5951,5 @@
 }
 /*}}}*/
-void    Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/
+void       Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/
 
 	/*diverse:*/
@@ -6001,5 +6001,5 @@
 }
 /*}}}*/
-void    Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/
+void       Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/
 
 	/*diverse:*/
@@ -6043,5 +6043,5 @@
 }
 /*}}}*/
-void    Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
+void       Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
 
 	/*diverse:*/
@@ -6132,5 +6132,5 @@
 }
 /*}}}*/
-void    Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/
+void       Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/
 
 	IssmDouble xyz_list[NUMVERTICES][3];
@@ -6238,6 +6238,6 @@
 /*}}}*/
 
-//new logic 
-void    Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
+//new code 
+void       Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
 
 	/*diverse:*/
@@ -6428,5 +6428,5 @@
 }
 /*}}}*/
-void    Tria::SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix<IssmDouble>* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/
+void       Tria::SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){ /*{{{*/
 
 	/*diverse:*/
@@ -6436,5 +6436,4 @@
 	IssmDouble I,W,BP;  //change in ice thickness, water column or bottom pressure (Farrel and Clarke, Equ. 4)
 	bool notfullygrounded=false;
-	bool scaleoceanarea= false;
 	bool computerigid= false;
 	int  glfraction=1;
@@ -6469,7 +6468,4 @@
 			constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
 		#endif
-		/*barystatic_contribution[0]+=0;
-		barystatic_contribution[1]+=0;
-		barystatic_contribution[2]+=0;*/
 		return;
 		}
@@ -6483,7 +6479,4 @@
 			this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
 			#endif
-			/*barystatic_contribution[0]+=0;
-			barystatic_contribution[1]+=0;
-			barystatic_contribution[2]+=0;*/
 			return;
 		}
@@ -6494,7 +6487,4 @@
 	if(!isice  && !ishydro){
 		if(!masks->isoceanin[this->lid]){
-			/*barystatic_contribution[0]+=0;
-			barystatic_contribution[1]+=0;
-			barystatic_contribution[2]+=0;*/
 			return;
 		}
@@ -6509,5 +6499,4 @@
 	rho_water=FindParam(MaterialsRhoSeawaterEnum);
 	this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
-	this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
 	this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
 	this->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
@@ -6518,7 +6507,7 @@
 
 	/*Deal with ice loads if we are on grounded ice:*/
-	if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){ /*{{{*/
-
-		/*Compute fraction of the element that is grounded: */
+	if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){ 
+
+		/*Compute fraction of the element that is grounded: {{{*/
 		if(notfullygrounded){
 			IssmDouble xyz_list[NUMVERTICES][3];
@@ -6532,4 +6521,5 @@
 		}
 		else phi_ice=1.0;
+		/*}}}*/
 
 		/*Inform mask: */
@@ -6569,13 +6559,11 @@
 		/*}}}*/
 
-		/*Compute barystatic contribution:*/
-		_assert_(oceanarea>0.);
-		if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
-		bslcice = rho_ice*area*I/(oceanarea*rho_water);
+		/*Compute barystatic contribution in kg:*/
+		bslcice = rho_ice*area*I;
 		_assert_(!xIsNan<IssmDouble>(bslcice));
 
 		/*Transfer thickness change into kg/m^2:*/
 		I=I*rho_ice*phi_ice; 
-	} /*}}}*/
+	} 
 
 	/*Deal with water loads if we are on ground:*/
@@ -6587,8 +6575,6 @@
 		deltathickness_input->GetInputAverage(&W);
 
-		/*Compute barystatic component:*/
-		_assert_(oceanarea>0.);
-		if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
-		bslchydro = rho_freshwater*area*phi_water*W/(oceanarea*rho_water);
+		/*Compute barystatic component in kg:*/
+		bslchydro = rho_freshwater*area*phi_water*W;
 		_assert_(!xIsNan<IssmDouble>(bslchydro));
 
@@ -6605,6 +6591,6 @@
 		bottompressure_change_input->GetInputAverage(&BP);
 		
-		/*Compute barystatic component:*/
-		bslcbp = rho_water*area*BP/(oceanarea*rho_water);
+		/*Compute barystatic component in kg:*/
+		bslcbp = rho_water*area*BP;
 
 		/*convert from m to kg/m^2:*/
@@ -6613,26 +6599,17 @@
 
 	/*Plug all loads into total load vector:*/
-	localloads[this->lid]+=I+W+BP;
+	loads->SetValue(this->sid,I+W+BP,INS_VAL);
 	
-	/*Plug bslcice into barystatic contribution vector:*/
-	if(barystatic_contribution_onpartition){
-		int idi=reCast<int>(partition[this->Sid()])+1;
-		int idj=0;
-		idj=0;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcice,ADD_VAL);
-		idj=1;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslchydro,ADD_VAL);
-		idj=2;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcbp,ADD_VAL);
-	}
-
-	barystatic_contribution[0]+=bslcice;
-	barystatic_contribution[1]+=bslchydro;
-	barystatic_contribution[2]+=bslcbp;
-
-}
-/*}}}*/
-IssmDouble Tria::SealevelchangeConvolution(IssmDouble* loads){ /*{{{*/
+	/*Keep track of barystatic contributions:*/
+	barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
+
+}
+/*}}}*/
+void       Tria::SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){ /*{{{*/
 
 	/*sal green function:*/
 	IssmDouble* G=NULL;
-	IssmDouble Sealevel[NUMVERTICES]={0,0,0};
+	IssmDouble SealevelGRD[NUMVERTICES]={0,0,0};
+	IssmDouble oceanaverage,oceanarea=0;
 
 	bool sal = false;
@@ -6647,13 +6624,113 @@
 		for(int i=0;i<NUMVERTICES;i++) {
 			for (int e=0;e<nel;e++){
-				Sealevel[i]+=G[i*nel+e]*loads[e]; 
-			}
-		}
-
-		this->AddInput(SealevelRSLEnum,Sealevel,P1Enum);
-	}
-
+				SealevelGRD[i]+=G[i*nel+e]*loads[e]; 
+			}
+		}
+	}
+
+	/*compute ocean average over element:*/
+	OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
+	
+	/*add ocean average in the global sealevelloads vector:*/
+	sealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
+	oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
+	
+	return;
 } /*}}}*/
-
+void       Tria::SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){ /*{{{*/
+
+	bool converged=false;
+	IssmDouble SealevelGRD[3]={0,0,0};
+	IssmDouble oceanaverage,oceanarea=0;
+	int nel;
+	bool sal = false;
+	IssmDouble* G=NULL;
+	int size;
+	
+	this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
+	this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
+		
+	if(sal){
+		this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
+
+		for(int i=0;i<NUMVERTICES;i++) {
+			for (int e=0;e<nel;e++){
+				SealevelGRD[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]);
+			}
+		}
+	}
+	OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
+	newsealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
+
+} /*}}}*/
+void       Tria::OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
+
+	IssmDouble phi=1.0;
+	bool iscoastline=false;
+	IssmDouble area;
+	IssmDouble Sg_avg=0; //output
+	
+	/*Do we have an ocean?:*/
+	if(!masks->isoceanin[this->lid]){
+		*poceanarea=0;
+		*poceanaverage=0;
+	}
+
+	/*Do we have a coastline?:*/
+	if(!masks->isfullyfloating[this->lid])iscoastline=true;
+
+	if(iscoastline){
+		IssmDouble xyz_list[NUMVERTICES][3];
+		::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+		phi=1.0-this->GetGroundedPortion(&xyz_list[0][0]); 
+	}
+
+	/*Get area of element:*/
+	this->Element::GetInputValue(&area,AreaEnum);
+
+	/*Average over ocean if there is no coastline, area of the ocean 
+	 *is the areaa of the element:*/
+	if(!iscoastline){
+
+		/*Average Sg over vertices:*/
+		for(int i=0;i<NUMVERTICES;i++) Sg_avg+=Sg[i]/NUMVERTICES;
+
+		*poceanaverage=Sg_avg;
+		*poceanarea=area;
+		return;
+	}
+
+	/*Average over  the ocean only if there is a coastline. Area of the 
+	 * ocean will be the fraction times the area of the element:*/
+	area=phi*area;
+   
+	IssmDouble total_weight=0;
+	bool mainlyfloating = true;
+	int         point1;
+	IssmDouble  fraction1,fraction2;
+
+	/*Recover portion of element that is grounded*/
+	this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+	//!mainlyfloating so that the integration is done on the ocean (and not the grounded) part
+	Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,!mainlyfloating,2); 
+
+	/* Start  looping on the number of gaussian points and average over these gaussian points: */
+	total_weight=0;
+	Sg_avg=0;
+	while(gauss->next()){
+		IssmDouble Sg_gauss=0;
+		TriaRef::GetInputValue(&Sg_gauss, Sg, gauss,P1Enum);
+		Sg_avg+=Sg_gauss*gauss->weight;
+		total_weight+=gauss->weight;
+	}
+	Sg_avg=Sg_avg/total_weight;
+	delete gauss;
+
+	*poceanaverage=Sg_avg;
+	*poceanarea=area;
+	return;
+
+}
+/*}}}*/
 #endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 26110)
@@ -162,18 +162,20 @@
 		#endif
 		#ifdef _HAVE_SEALEVELCHANGE_
+		void       SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
+		void       SealevelchangeGeometry(IssmDouble* lat, IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
+		IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea);
+		IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea);
+		void       SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks);
+		void       SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks);
+		void       DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks);
+		void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
 		IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks);
-		void    SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
-		void    SetSealevelMasks(SealevelMasks* masks);
-		void    SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
-		void    SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
-
-		IssmDouble    SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea);
-		IssmDouble    SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea);
-		void    SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks);
-		void    SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks);
-		void    DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks);
-		void    GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
-		void    SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix<IssmDouble>* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea);
-		IssmDouble SealevelchangeConvolution(IssmDouble* loads);
+		void       SetSealevelMasks(SealevelMasks* masks);
+		
+		void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
+		void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks);
+		void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks);
+		void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks);
+		void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks);
 		#endif
 		/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 26109)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 26110)
@@ -4749,162 +4749,4 @@
 #endif
 #ifdef _HAVE_SEALEVELCHANGE_
-void FemModel::SealevelchangeBarystatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition,SealevelMasks* masks) { /*{{{*/
-
-	/*serialized vectors:*/
-	IssmDouble  bslcice       = 0.;
-	IssmDouble  bslcice_cpu   = 0.;
-	IssmDouble  bslchydro       = 0.;
-	IssmDouble  bslchydro_cpu   = 0.;
-	IssmDouble  area      = 0.;
-	IssmDouble  oceanarea      = 0.;
-	IssmDouble  oceanarea_cpu  = 0.;
-	int bp_compute_fingerprints= 0;
-	bool isoceantransport=false;
-
-	Vector<IssmDouble>* bslcice_partition=NULL;
-	IssmDouble* bslcice_partition_serial=NULL;
-	IssmDouble* partitionice=NULL;
-	int npartice,nel;
-
-	Vector<IssmDouble>* bslchydro_partition=NULL;
-	IssmDouble* bslchydro_partition_serial=NULL;
-	IssmDouble* partitionhydro=NULL;
-	bool istws=0;
-	int nparthydro;
-		
-	int npartocean;
-	Vector<IssmDouble>* bslcocean_partition=NULL;
-	IssmDouble* bslcocean_partition_serial=NULL;
-	IssmDouble* partitionocean=NULL;
-
-   /*Initialize temporary vector that will be used to sum barystatic components
-    * on all local elements, prior to assembly:*/
-	int gsize = this->nodes->NumberOfDofs(GsetEnum);
-	IssmDouble* RSLgi=xNewZeroInit<IssmDouble>(gsize);
-	int* indices=xNew<int>(gsize);
-   for(int i=0;i<gsize;i++) indices[i]=i;
-
-	/*First, figure out the area of the ocean, which is needed to compute the barystatic component: */
-	int i = -1;
-	for(Object* & object : this->elements->objects){
-		i +=1;
-		Element* element = xDynamicCast<Element*>(object);
-		element->GetInputValue(&area,AreaEnum);
-		if (masks->isoceanin[i]) oceanarea_cpu += area;
-	}
-	ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	_assert_(oceanarea>0.);
-
-	/*Initialize partition vectors to retrieve barystatic contributions: */
-	this->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
-	if(npartice){
-		this->parameters->FindParam(&partitionice,&nel,NULL,SolidearthPartitionIceEnum);
-		bslcice_partition= new Vector<IssmDouble>(npartice);
-	}
-
-	this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
-	if(nparthydro){
-		this->parameters->FindParam(&partitionhydro,&nel,NULL,SolidearthPartitionHydroEnum);
-		bslchydro_partition= new Vector<IssmDouble>(nparthydro);
-	}
-
-	this->parameters->FindParam(&npartocean,SolidearthNpartOceanEnum);
-	if(npartocean){
-		this->parameters->FindParam(&partitionocean,&nel,NULL,SolidearthPartitionOceanEnum);
-		bslchydro_partition= new Vector<IssmDouble>(npartocean);
-	}
-	/*For later:
-	npartbarystatic=npartice;
-	if(nparthydro>npartbarystatic)npartbarystatic=nparthydro;
-	if(npartocean>npartbarystatic)npartbarystatic=npartocean;
-	bslc_partition=new Matrix(IssmDouble>(npartbarystatic,3);
-
-	bslc_cpu[0]=0; bslc_cpu[1]=0; bslc_cpu[2]=0;
-	for(Object* & object : this->elements->objects){
-		Element* element = xDynamicCast<Element*>(object);
-		element->SealevelchangeBarystaticLoads(&bslc_cpu[0], localloads,masks, bslcice_partition,partitionice,oceanarea);
-	}
-	MPI Bcast localloads -> loads
-
-	for(Object* & object : this->elements->objects){
-		Element* element = xDynamicCast<Element*>(object);
-		element->SealevelchangeConvolution(loads);
-	}
-	*/
-
-
-
-
-
-	/*Call the barystatic sea level change core for ice : */
-	bslcice_cpu=0;
-	for(Object* & object : this->elements->objects){
-		Element* element = xDynamicCast<Element*>(object);
-		bslcice_cpu+=element->SealevelchangeBarystaticIce(RSLgi,masks, bslcice_partition,partitionice,oceanarea);
-	}
-
-	/*Call the barystatic sea level change core for hydro: */
-	bslchydro_cpu=0; //make sure to initialize this, so we have a total barystatic contribution computed at 0.
-	this->parameters->FindParam(&istws,TransientIshydrologyEnum);
-	if(istws){
-		for(int i=0;i<elements->Size();i++){
-			Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-			bslchydro_cpu+=element->SealevelchangeBarystaticHydro(RSLgi,masks, bslchydro_partition,partitionhydro,oceanarea);
-		}
-	}
-
-	/*Call the barystatic sea level change core for bottom pressures: */
-	this->parameters->FindParam(&bp_compute_fingerprints,SolidearthSettingsComputeBpGrdEnum);
-	this->parameters->FindParam(&isoceantransport,TransientIsoceantransportEnum);
-	if(bp_compute_fingerprints && isoceantransport){
-		for(int i=0;i<elements->Size();i++){
-			Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-			element->SealevelchangeBarystaticBottomPressure(RSLgi,masks);
-		}
-	}
-
-	/*Plug values once and assemble: */
-	pRSLgi->SetValues(gsize,indices,RSLgi,ADD_VAL);
-	pRSLgi->Assemble();
-
-	/*Sum all barystatic components from all cpus:*/
-	ISSM_MPI_Reduce (&bslcice_cpu,&bslcice,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&bslcice,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	_assert_(!xIsNan<IssmDouble>(bslcice));
-
-	ISSM_MPI_Reduce (&bslchydro_cpu,&bslchydro,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&bslchydro,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	_assert_(!xIsNan<IssmDouble>(bslchydro));
-
-	/*Take care of partition vectors:*/
-	if(bslcice_partition){
-		bslcice_partition->Assemble();
-		bslcice_partition_serial=bslcice_partition->ToMPISerial();
-	}
-	if(bslchydro_partition){
-		bslchydro_partition->Assemble();
-		bslchydro_partition_serial=bslchydro_partition->ToMPISerial();
-	}
-
-
-	/*Free ressources:*/
-	xDelete<int>(indices);
-	xDelete<IssmDouble>(RSLgi);
-	if(bslchydro_partition)delete bslchydro_partition;
-	if(bslcice_partition)delete bslcice_partition;
-	if(partitionhydro)xDelete<IssmDouble>(partitionhydro);
-	if(partitionice)xDelete<IssmDouble>(partitionice);
-
-	/*Assign output pointers:*/
-	*poceanarea = oceanarea;
-	*pbslcice  = bslcice;
-	*pbslchydro  = bslchydro;
-	*pbslc=bslchydro+bslcice;
-	*pbslcice_partition=bslcice_partition_serial;
-	*pbslchydro_partition=bslchydro_partition_serial;
-
-}
-/*}}}*/
 void FemModel::SealevelchangeSal(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old,  SealevelMasks* masks, bool verboseconvolution){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/classes/classes.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/classes.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/classes/classes.h	(revision 26110)
@@ -19,4 +19,5 @@
 #include "./Misfit.h"
 #include "./SealevelMasks.h"
+#include "./BarystaticContributions.h"
 #include "./Nodalvalue.h"
 #include "./Numberedcostfunction.h"
Index: /issm/trunk-jpl/src/c/cores/cores.h
===================================================================
--- /issm/trunk-jpl/src/c/cores/cores.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/cores/cores.h	(revision 26110)
@@ -81,8 +81,5 @@
 void ResetBoundaryConditions(FemModel* femmodel, int analysis_type);
 void PrintBanner(void);
-void TransferForcing(FemModel* femmodel,int forcingenum);
-void TransferSealevel(FemModel* femmodel,int forcingenum);
 void EarthMassTransport(FemModel* femmodel);
-void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
 
 //solution configuration
Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26109)
+++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26110)
@@ -12,4 +12,10 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+/*support routines local definitions:{{{*/
+void TransferForcing(FemModel* femmodel,int forcingenum);
+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);
+/*}}}*/
 
 /*main cores:*/
@@ -383,4 +389,105 @@
 	}
 }; /*}}}*/
+
+void grd_core_optim(FemModel* femmodel,SealevelMasks* masks) { /*{{{*/
+
+	/*variables:{{{*/
+	int nel;
+	BarystaticContributions* barycontrib=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;
+	bool                   scaleoceanarea=false;
+	IssmDouble             rho_water;
+
+	IssmDouble           eps_rel;
+	IssmDouble           eps_abs;
+	int                  step;
+	IssmDouble           time; 
+	
+	IssmDouble cumbslc;
+	IssmDouble cumbslcice;
+	IssmDouble cumbslchydro;
+	/*}}}*/
+
+	/*retrieve parameters: */
+	femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
+	barycontribparam = xDynamicCast<GenericParam<BarystaticContributions*>*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum));
+	barycontrib=barycontribparam->GetParameterValue();
+	femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
+	femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
+	femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
+	
+	/*initialize matrices and vectors:*/
+	femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
+	loads=new Vector<IssmDouble>(nel);
+	sealevelloads=new Vector<IssmDouble>(nel);
+	oceanareas=new Vector<IssmDouble>(nel);
+	
+	/*buildup loads: */
+	for(Object* & object : femmodel->elements->objects){
+		Element* element = xDynamicCast<Element*>(object);
+		element->SealevelchangeBarystaticLoads(loads, barycontrib,masks); 
+	}
+
+	//Communicate loads from local to global: 
+	loads->Assemble(); allloads=loads->ToMPISerial();
+
+	/*convolve loads:*/
+	for(Object* & object : femmodel->elements->objects){
+		Element* element = xDynamicCast<Element*>(object);
+		element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,masks);
+	}
+
+	//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);
+	allsealevelloads=sealevelloads->ToMPISerial();
+
+	bool converged=false;
+	for(;;){
+			
+		oldsealevelloads=sealevelloads->Duplicate();
+
+		/*convolve load and sealevel loads on oceans:*/
+		for(Object* & object : femmodel->elements->objects){
+			Element* element = xDynamicCast<Element*>(object);
+			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 
+		allsealevelloads=sealevelloads->ToMPISerial();
+
+		//convergence?
+		slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs);
+		if (converged)break;
+	}
+
+	/*cumulate barystatic contributions and save to results: */
+	barycontrib->Cumulate(femmodel->parameters);
+	barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
+
+}
+/*}}}*/
 
 //Geometry:
@@ -1062,2 +1169,10 @@
 
 } /*}}}*/
+IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea){ /*{{{*/
+	IssmDouble sealevelloadsaverage;	
+	Vector<IssmDouble>* sealevelloadsvolume=sealevelloads->Duplicate();
+	sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas);
+	sealevelloadsvolume->Sum(&sealevelloadsaverage);
+	delete sealevelloadsvolume; 
+	return sealevelloadsaverage/oceanarea;
+} /*}}}*/
Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26109)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26110)
@@ -747,4 +747,5 @@
 syn keyword cConstant BslcIceEnum
 syn keyword cConstant BslcHydroEnum
+syn keyword cConstant BslcOceanEnum
 syn keyword cConstant BslcRateEnum
 syn keyword cConstant GmtslcEnum
@@ -1071,4 +1072,5 @@
 syn keyword cConstant BasalforcingsIsmip6Enum
 syn keyword cConstant BasalforcingsPicoEnum
+syn keyword cConstant BarystaticContributionsEnum
 syn keyword cConstant BeckmannGoosseFloatingMeltRateEnum
 syn keyword cConstant BedSlopeSolutionEnum
@@ -1433,4 +1435,5 @@
 syn keyword cType AmrNeopz
 syn keyword cType ArrayInput
+syn keyword cType BarystaticContributions
 syn keyword cType BoolInput
 syn keyword cType BoolParam
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26110)
@@ -743,4 +743,5 @@
 	BslcIceEnum,
 	BslcHydroEnum,
+	BslcOceanEnum,
 	BslcRateEnum,
 	GmtslcEnum,
@@ -1070,4 +1071,5 @@
 	BasalforcingsIsmip6Enum,
 	BasalforcingsPicoEnum,
+	BarystaticContributionsEnum,
 	BeckmannGoosseFloatingMeltRateEnum,
 	BedSlopeSolutionEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26109)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26110)
@@ -749,4 +749,5 @@
 		case BslcIceEnum : return "BslcIce";
 		case BslcHydroEnum : return "BslcHydro";
+		case BslcOceanEnum : return "BslcOcean";
 		case BslcRateEnum : return "BslcRate";
 		case GmtslcEnum : return "Gmtslc";
@@ -1073,4 +1074,5 @@
 		case BasalforcingsIsmip6Enum : return "BasalforcingsIsmip6";
 		case BasalforcingsPicoEnum : return "BasalforcingsPico";
+		case BarystaticContributionsEnum : return "BarystaticContributions";
 		case BeckmannGoosseFloatingMeltRateEnum : return "BeckmannGoosseFloatingMeltRate";
 		case BedSlopeSolutionEnum : return "BedSlopeSolution";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26109)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26110)
@@ -767,4 +767,5 @@
 	      else if (strcmp(name,"BslcIce")==0) return BslcIceEnum;
 	      else if (strcmp(name,"BslcHydro")==0) return BslcHydroEnum;
+	      else if (strcmp(name,"BslcOcean")==0) return BslcOceanEnum;
 	      else if (strcmp(name,"BslcRate")==0) return BslcRateEnum;
 	      else if (strcmp(name,"Gmtslc")==0) return GmtslcEnum;
@@ -874,9 +875,9 @@
 	      else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
 	      else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
-	      else if (strcmp(name,"SmbT")==0) return SmbTEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
+	      if (strcmp(name,"SmbT")==0) return SmbTEnum;
+	      else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
 	      else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
 	      else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
@@ -997,9 +998,9 @@
 	      else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
 	      else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
-	      else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
+	      if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
+	      else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
 	      else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
 	      else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
@@ -1097,4 +1098,5 @@
 	      else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum;
 	      else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum;
+	      else if (strcmp(name,"BarystaticContributions")==0) return BarystaticContributionsEnum;
 	      else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
 	      else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
@@ -1119,10 +1121,10 @@
 	      else if (strcmp(name,"Closed")==0) return ClosedEnum;
 	      else if (strcmp(name,"Colinear")==0) return ColinearEnum;
-	      else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
-	      else if (strcmp(name,"Contact")==0) return ContactEnum;
          else stage=10;
    }
    if(stage==10){
-	      if (strcmp(name,"Contour")==0) return ContourEnum;
+	      if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
+	      else if (strcmp(name,"Contact")==0) return ContactEnum;
+	      else if (strcmp(name,"Contour")==0) return ContourEnum;
 	      else if (strcmp(name,"Contours")==0) return ContoursEnum;
 	      else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
@@ -1242,10 +1244,10 @@
 	      else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
 	      else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
-	      else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
-	      else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
          else stage=11;
    }
    if(stage==11){
-	      if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
+	      if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
+	      else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
+	      else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
 	      else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
 	      else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum;
@@ -1365,10 +1367,10 @@
 	      else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum;
 	      else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;
-	      else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
-	      else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
          else stage=12;
    }
    if(stage==12){
-	      if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
+	      if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
+	      else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
+	      else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
 	      else if (strcmp(name,"Scaled")==0) return ScaledEnum;
 	      else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h	(revision 26110)
@@ -52,4 +52,5 @@
 		virtual void PointwiseMult(IssmAbsVec* x,IssmAbsVec* y)=0;
 		virtual void Pow(doubletype scale_factor)=0;
+		virtual void Sum(doubletype* pvalue)=0;
 };
 
Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h	(revision 26110)
@@ -594,4 +594,8 @@
 		}
 		/*}}}*/
+		void Sum(doubletype* pvalue){/*{{{*/
+			_error_("not support yet!");
+		}
+		/*}}}*/
 		void BucketsBuildScatterBuffers(int** pnumvalues_forcpu,int** prow_indices_forcpu,doubletype** pvalues_forcpu,int** pmodes_forcpu,DataSet** bucketsforcpu,int num_procs){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h	(revision 26110)
@@ -324,4 +324,12 @@
 		}
 		/*}}}*/
+		void Sum(doubletype* pvalue){/*{{{*/
+
+			doubletype value=0;
+			int i;
+			for(i=0;i<this->M;i++)value+=this->vector[i];
+			*pvalue=value;
+		}
+		/*}}}*/
 		
 };
Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h	(revision 26110)
@@ -220,4 +220,8 @@
 		}
 		/*}}}*/
+		void Sum(doubletype*  pvalue){/*{{{*/
+			vector->Sum(pvalue);
+		}
+		/*}}}*/
 };
 
Index: /issm/trunk-jpl/src/c/toolkits/objects/Matrix.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/objects/Matrix.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/toolkits/objects/Matrix.h	(revision 26110)
@@ -282,4 +282,20 @@
 		}
 		/*}}}*/
+		doubletype* ToMPISerial(void){/*{{{*/
+
+			doubletype* output=NULL;
+
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				output=this->pmatrix->ToMPISerial();
+				#endif
+			}
+			else{
+				_error_("not implemented yet!");
+			}
+
+			return output;
+		}
+		/*}}}*/
 		void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){/*{{{*/
 
@@ -307,8 +323,6 @@
 		}
 		/*}}}*/
-		/*
-		* sets all values to 0 but keeps the structure of a sparse matrix
-		*/
 		void SetZero(void) {/*{{{*/
+			// sets all values to 0 but keeps the structure of a sparse matrix
 			if(type==PetscMatType){
 				#ifdef _HAVE_PETSC_
Index: /issm/trunk-jpl/src/c/toolkits/objects/Vector.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/objects/Vector.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/toolkits/objects/Vector.h	(revision 26110)
@@ -406,4 +406,15 @@
 		}
 		/*}}}*/
-};
+void Sum(doubletype* pvalue){ /*{{{*/
+	_assert_(this);/*{{{*/
+
+	if(type==PetscVecType){
+		#ifdef _HAVE_PETSC_
+		this->pvector->Sum(pvalue);
+		#endif
+	}
+	else this->ivector->Sum(pvalue);
+}
+/*}}}*/
+}; /*}}}*/
 #endif //#ifndef _VECTOR_H_
Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp	(revision 26109)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp	(revision 26110)
@@ -249,4 +249,11 @@
 }
 /*}}}*/
+void PetscVec::Sum(IssmDouble* pvalue){/*{{{*/
+
+	_assert_(this->vector);
+	VecSum(this->vector,pvalue);
+
+}
+/*}}}*/
 IssmDouble PetscVec::Dot(PetscVec* input){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h	(revision 26110)
@@ -58,4 +58,5 @@
 		void        Scale(IssmDouble scale_factor);
 		void        Pow(IssmDouble scale_factor);
+		void        Sum(IssmDouble* pvalue);
 		void        PointwiseDivide(PetscVec* x,PetscVec* y);
 		void        PointwiseMult(PetscVec* x,PetscVec* y);
Index: /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h	(revision 26109)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h	(revision 26110)
@@ -32,4 +32,5 @@
 void MatMultPatch(Mat A,Vec X, Vec AX,ISSM_MPI_Comm comm);
 void MatToSerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
+void MatToMPISerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
 Vec  SerialToVec(double* vector,int vector_size);
 InsertMode ISSMToPetscInsertMode(InsMode mode);
Index: /issm/trunk-jpl/test/NightlyRun/test2004.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2004.m	(revision 26109)
+++ /issm/trunk-jpl/test/NightlyRun/test2004.m	(revision 26110)
@@ -134,4 +134,8 @@
 		disp('      reading bedrock');
 		md.geometry.bed=-ones(md.mesh.numberofvertices,1);
+		md.geometry.base=md.geometry.bed;
+		md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1);
+		md.geometry.surface=md.geometry.bed+md.geometry.thickness;
+
 	end % }}}
 	%Slc: {{{
@@ -167,9 +171,9 @@
 		end
 
-		md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
-
-		md.dsl.global_average_thermosteric_sea_level_change=[0;0];
-		md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
-		md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
+		md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
+
+		md.dsl.global_average_thermosteric_sea_level=[0;0];
+		md.dsl.sea_surface_height_above_geoid=zeros(md.mesh.numberofvertices+1,1);
+		md.dsl.sea_water_pressure_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
 
 	end %}}}
@@ -313,4 +317,7 @@
 	di=md.materials.rho_ice/md.materials.rho_water;
 	md.geometry.bed=-ones(md.mesh.numberofvertices,1);
+	md.geometry.base=md.geometry.bed;
+	md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1);
+	md.geometry.surface=md.geometry.bed+md.geometry.thickness;
 	% }}}
 	%materials:  {{{
@@ -346,4 +353,7 @@
 sl.transfer('mask.ocean_levelset');
 sl.transfer('geometry.bed');
+sl.transfer('geometry.surface');
+sl.transfer('geometry.thickness');
+sl.transfer('geometry.base');
 sl.transfer('mesh.lat');
 sl.transfer('mesh.long');
@@ -400,4 +410,14 @@
 md.transient.isslc=1;
 
+%Initializations: 
+md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
+md.initialization.bottompressure=zeros(md.mesh.numberofvertices,1);
+md.initialization.dsl=zeros(md.mesh.numberofvertices,1);
+md.initialization.str=0;
+md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
 
 %max number of iterations reverted back to 10 (i.e. the original default value)