Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 26113)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 26114)
@@ -4749,4 +4749,162 @@
 #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/FemModel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 26113)
+++ /issm/trunk-jpl/src/c/classes/FemModel.h	(revision 26114)
@@ -164,4 +164,5 @@
 		#endif
 		#ifdef _HAVE_SEALEVELCHANGE_
+		void SealevelchangeBarystatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition, SealevelMasks* masks); 
 		void SealevelchangeSal(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old,  SealevelMasks* masks,bool verboseconvolution);
 		void SealevelchangeRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks);