Index: issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 26114)
+++ issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 26165)
@@ -4748,408 +4748,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){/*{{{*/
-
-	/*serialized vectors:*/
-	IssmDouble* RSLg_old=NULL;
-
-	IssmDouble* RSLgo  = NULL;
-	int* indices = NULL;
-	int         gsize;
-
-	bool computerigid = true;
-	bool computeelastic= true;
-
-	/*recover computational flags: */
-	this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
-
-	/*Initialize temporary vector that will be used to sum barystatic components on all local elements, prior
-	 * to assembly:*/
-	gsize = this->nodes->NumberOfDofs(GsetEnum);
-	RSLgo=xNewZeroInit<IssmDouble>(gsize);
-	indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
-
-	/*Serialize vectors from previous iteration:*/
-	RSLg_old=pRSLg_old->ToMPISerial();
-
-	/*Call the sal sea level change core only if required: */
-	if(computerigid){
-		for(Object* & object : this->elements->objects){
-			Element* element = xDynamicCast<Element*>(object);
-			element->SealevelchangeSal(RSLgo,RSLg_old,masks);
-		}
-	}
-	pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL);
-	pRSLgo->Assemble();
-
-	/*Free ressources:*/
-	xDelete<int>(indices);
-	xDelete<IssmDouble>(RSLgo);
-	xDelete<IssmDouble>(RSLg_old);
-
-}
-/*}}}*/
-void FemModel::SealevelchangeRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz,  SealevelMasks* masks){/*{{{*/
-
-	/*serialized vectors:*/
-	bool spherical=true;
-	IssmDouble* RSLg_old=NULL;
-	IssmDouble*	tide_love_h  = NULL;
-	IssmDouble*	tide_love_k  = NULL;
-	IssmDouble*	load_love_k  = NULL;
-	IssmDouble  tide_love_k2secular;
-	IssmDouble  moi_e, moi_p, omega, g;
-	IssmDouble	m1, m2, m3;
-	IssmDouble	lati, longi, radi, value;
-	IssmDouble          *latitude    = NULL;
-	IssmDouble          *longitude    = NULL;
-	IssmDouble          *radius    = NULL;
-
-	/*Serialize vectors from previous iteration:*/
-	RSLg_old=pRSLg_old->ToMPISerial();
-
-	IssmDouble moi_list[3]={0,0,0};
-	IssmDouble moi_list_cpu[3]={0,0,0};
-	for(Object* & object : this->elements->objects){
-		Element* element = xDynamicCast<Element*>(object);
-		element->SealevelchangeMomentOfInertia(&moi_list[0],RSLg_old,masks );
-		moi_list_cpu[0] += moi_list[0];
-		moi_list_cpu[1] += moi_list[1];
-		moi_list_cpu[2] += moi_list[2];
-	}
-	ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	//
-	ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-	//
-	ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
-	/*pull out some useful parameters: */
-	parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
-	parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
-	parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
-	parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
-	parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
-	parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
-	parameters->FindParam(&omega,RotationalAngularVelocityEnum);
-
-	/*compute perturbation terms for angular velocity vector: */
-	m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0];
-	m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1];
-	m3 = -(1+load_love_k[2])/moi_p * moi_list[2];	// term associated with fluid number (3-order-of-magnitude smaller) is negelected
-
-	/*recover lat,long and radius vectors from vertices: */
-	VertexCoordinatesx(&latitude,&longitude,&radius,this->vertices,spherical);
-
-	/* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11)
-	 * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3)
-	 * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered)
-	 * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180]
-	 */
-	for(Object* & object : vertices->objects){
-		Vertex* vertex=xDynamicCast<Vertex*>(object);
-		int sid=vertex->Sid();
-
-		lati=latitude[sid]/180*PI;
-		longi=longitude[sid]/180*PI;
-		radi=radius[sid];
-
-		/*only first order terms are considered now: */
-		value=((1.0+tide_love_k[2]-tide_love_h[2])/9.81)*pow(omega*radi,2.0)*
-						(-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi)));
-
-		pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
-	}
-
-	/*Assemble mesh velocity*/
-	pRSLgo_rot->Assemble();
-
-	/*Assign output pointers:*/
-	if(pIxz)*pIxz=moi_list[0];
-	if(pIyz)*pIyz=moi_list[1];
-	if(pIzz)*pIzz=moi_list[2];
-	xDelete<IssmDouble>(latitude);
-	xDelete<IssmDouble>(longitude);
-	xDelete<IssmDouble>(tide_love_h);
-	xDelete<IssmDouble>(tide_love_k);
-	xDelete<IssmDouble>(load_love_k);
-
-	xDelete<IssmDouble>(radius);
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(RSLg_old);
-
-}
-/*}}}*/
-void FemModel::SealevelchangeDeformation(Vector<IssmDouble>* pgeoid, Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, SealevelMasks* masks){/*{{{*/
-
-	/*serialized vectors:*/
-	IssmDouble* RSLg=NULL;
-
-	IssmDouble* Up  = NULL;
-	IssmDouble* North  = NULL;
-	IssmDouble* East  = NULL;
-	int* indices = NULL;
-	int  gsize;
-	int  horiz;
-
-	/*retrieve parameters:*/
-	this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
-
-	/*Serialize vectors from previous iteration:*/
-	RSLg=pRSLg->ToMPISerial();
-
-	/*Initialize temporary vector that will be used to sum barystatic components on all local elements, prior
-	 * to assembly:*/
-	gsize = this->nodes->NumberOfDofs(GsetEnum);
-	Up=xNewZeroInit<IssmDouble>(gsize);
-	if(horiz){
-		North=xNewZeroInit<IssmDouble>(gsize);
-		East=xNewZeroInit<IssmDouble>(gsize);
-	}
-	indices=xNew<int>(gsize); for (int i=0;i<gsize;i++)indices[i]=i;
-
-	/*Call the deformation from loading routines:*/
-	for(Object* & object : this->elements->objects){
-		Element* element = xDynamicCast<Element*>(object);
-		element->DeformationFromSurfaceLoads(Up,North,East,RSLg,masks);
-	}
-
-	pUp->SetValues(gsize,indices,Up,ADD_VAL);
-	pUp->Assemble();
-
-
-	/*Add RSL to Up to find the geoid: */
-	pUp->Copy(pgeoid); pgeoid->AXPY(pRSLg,1);
-
-	if (horiz){
-		pNorth->SetValues(gsize,indices,North,ADD_VAL);
-		pNorth->Assemble();
-		pEast->SetValues(gsize,indices,East,ADD_VAL);
-		pEast->Assemble();
-	}
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(Up);
-	if(horiz){
-		xDelete<IssmDouble>(North);
-		xDelete<IssmDouble>(East);
-	}
-	xDelete<int>(indices);
-	xDelete<IssmDouble>(RSLg);
-}
-/*}}}*/
-IssmDouble FemModel::SealevelchangeOceanAverage(Vector<IssmDouble>* RSLg,SealevelMasks* masks, IssmDouble oceanarea) { /*{{{*/
-
-	IssmDouble* RSLg_serial=NULL;
-	IssmDouble  oceanvalue,oceanvalue_cpu;
-
-	/*Serialize vectors from previous iteration:*/
-	RSLg_serial=RSLg->ToMPISerial();
-
-	/*Initialize:*/
-	oceanvalue_cpu=0;
-
-	/*Go through elements, and add contribution from each element and divide by overall ocean area:*/
-	for(Object* & object : this->elements->objects){
-		Element* element = xDynamicCast<Element*>(object);
-		oceanvalue_cpu += element->OceanAverage(RSLg_serial,masks);
-	}
-
-	ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-	ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(RSLg_serial);
-
-	return oceanvalue/oceanarea;
-}
-/*}}}*/
-void FemModel::IvinsDeformation(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, IssmDouble* x, IssmDouble* y){ /*{{{*/
-
-	/*Find the litho material to be used by all the elements:*/
-	Matlitho* matlitho=NULL;
-	for (Object* & object: this->materials->objects){
-		Material* material=xDynamicCast<Material*>(object);
-		if(material->ObjectEnum()==MatlithoEnum){
-			matlitho=xDynamicCast<Matlitho*>(material);
-			break;
-		}
-	}
-	_assert_(matlitho);
-
-
-	/*Go through elements, and add contribution from each element to the deflection vector wg:*/
-	for(Object* & object : this->elements->objects){
-		Element* element = xDynamicCast<Element*>(object);
-		element->GiaDeflection(wg,dwgdt, matlitho, x,y);
-	}
-
-	/*Assemble parallel vector:*/
-	dwgdt->Assemble();
-	wg->Assemble();
-}
-/*}}}*/
-#endif
 void FemModel::HydrologyEPLupdateDomainx(IssmDouble* pEplcount){ /*{{{*/
 
