/*!\file: sealevelchange_core.cpp
 * \brief: core of the sea-level change solution 
 */ 

#include "./cores.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../classes/Inputs/TriaInput.h"
#include "../classes/Inputs/TransientInput.h"
#include "../classes/Inputs/DatasetInput.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"


/*main cores:*/
void sealevelchange_core(FemModel* femmodel){ /*{{{*/

	/*Start profiler*/
	femmodel->profiler->Start(SLRCORE);

	/*Parameters, variables:*/
	bool save_results;
	bool isslr=0;
	bool isgia=0;
	int  grd=0;
	int solution_type;

	/*Retrieve parameters:*/
	femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
	femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); 

	/*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
	if(solution_type==SealevelriseSolutionEnum){
		isslr=1;
		isgia=1;
	}

	/*Should we be here?:*/
	if(!isslr)return;

	/*Verbose: */
	if(VerboseSolution()) _printf0_("   computing sea level rise\n");

	/*Run gia core: */
	if(isgia){
		#ifdef _HAVE_GIA_
		gia_core(femmodel);
		#else
		_error_("ISSM was not compiled with gia capabilities. Exiting");
		#endif
	}

	/*set SLR configuration: */
	femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);

	/*run geometry core: */
	sealevelrise_core_geometry(femmodel);

	/*Run geodetic:*/
	if(grd)grd_core(femmodel);

	/*Run steric core for sure:*/
	dynstr_core(femmodel);

	/*Save results: */
	if(save_results){
		int     numoutputs;
		char **requested_outputs = NULL;
		if(VerboseSolution()) _printf0_("   saving results\n");
		femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
		femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
		if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
	}

	/*requested dependents: */
	if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
	
	/*End profiler*/
	femmodel->profiler->Stop(SLRCORE);
}
/*}}}*/
void grd_core(FemModel* femmodel){ /*{{{*/

	/*Gravity rotation deformation core GRD: */

	/*variables:*/
	Vector<IssmDouble> *RSLg    = NULL;
	Vector<IssmDouble> *BPg    = NULL;
	Vector<IssmDouble> *RSLg_rate    = NULL;
	Vector<IssmDouble> *RSLg_eustatic  = NULL; 
	Vector<IssmDouble> *U_esa  = NULL; 
	Vector<IssmDouble> *U_esa_rate  = NULL; 
	Vector<IssmDouble> *N_esa  = NULL; 
	Vector<IssmDouble> *N_esa_rate  = NULL; 
	Vector<IssmDouble> *U_north_esa   = NULL; 
	Vector<IssmDouble> *U_east_esa    = NULL; 
	Vector<IssmDouble> *N_gia= NULL; 
	Vector<IssmDouble> *U_gia= NULL; 
	Vector<IssmDouble> *N_gia_rate= NULL; 
	Vector<IssmDouble> *U_gia_rate= NULL; 
	SealevelMasks* masks=NULL;

	/*parameters:*/
	bool iscoupler;
	int  solution_type;
	int  modelid,earthid;
	bool istransientmasstransport;
	int  frequency,count;
	int  horiz;
	IssmDouble dt;
	IssmDouble oceanarea;
	int        bp_compute_fingerprints=0;

	
	/*Verbose: */
	if(VerboseSolution()) _printf0_("	  computing GRD sea level patterns\n");

	/*retrieve more parameters:*/
	femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
	femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
	femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum);
	femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);

	if(iscoupler){
		femmodel->parameters->FindParam(&modelid,ModelIdEnum);
		femmodel->parameters->FindParam(&earthid,EarthIdEnum);
		femmodel->parameters->FindParam(&istransientmasstransport,TransientIsmasstransportEnum);
	}
	else{
		/* we are here, we are not running in a coupler, so we will indeed compute SLR,
		 * so make sure we are identified as being the Earth.:*/
		modelid=1; earthid=1; 
		/* in addition, if we are running solution_type SealevelriseSolutionEnum, make sure we 
		 * run, irresepective of the time settings:*/
		count=frequency;
	}

	/*If we are running in coupled mode, the Earth model needs to run its own mass transport (if 
	 * not already done by the mass trasnport module. For ice caps, they rely on the transient mass 
	 * transport module exclusively:*/
	if(iscoupler) if(modelid==earthid) if(!istransientmasstransport) EarthMassTransport(femmodel);

	/*increment counter, or call solution core if count==frequency:*/
	if (count<frequency){
		count++; femmodel->parameters->SetParam(count,SealevelriseRunCountEnum); 
		return;
	}

	/*call sea-level rise sub cores:*/
	if(iscoupler){
		/*transfer cumulated deltathickness forcing from ice caps to earth model: */
		TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum);

		/*we have accumulated thicknesses, dump them in deltathcikness: */
		if(modelid==earthid)InputDuplicatex(femmodel,SealevelriseCumDeltathicknessEnum,SurfaceloadIceThicknessChangeEnum);
	}

	/*run cores:*/
	if(modelid==earthid){

		/*call masks core: */
		masks=sealevelrise_core_masks(femmodel);

		/*set bottom pressures:*/
		SetBottomPressure(femmodel);


		/*call eustatic core  (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
		RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&oceanarea); 

		/*call non-eustatic core (ocean loading tems  - 2nd and 5th terms on the RHS of Farrel and Clark) */
		RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,oceanarea); 

		/*compute other elastic signatures, such as components of 3-D crustal motion: */
		sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,masks);

		/*compute viscosus (GIA) signatures:*/
		sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg);

		/*compute sea-level rise (low-order spherical harmonics coefficients) diagnostics:*/
		sealevelrise_diagnostics(femmodel,RSLg);

		/*recover N_esa  = U_esa + RSLg:*/
		N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1);

		/*transform these values into rates (as we only run this once each frequency turn:*/
		N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency));
		U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency));
		N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency));
		U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency));
		RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency));

		/*get some results into elements:{{{*/
		InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum); 
		InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum); 
		InputUpdateFromVectorx(femmodel,U_gia_rate,UGiaRateEnum,VertexSIdEnum);
		InputUpdateFromVectorx(femmodel,N_gia_rate,NGiaRateEnum,VertexSIdEnum); 
		InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum); 
		InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum); 
		InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum); 
		InputUpdateFromVectorx(femmodel,U_gia,UGiaEnum,VertexSIdEnum); 
		InputUpdateFromVectorx(femmodel,N_gia,NGiaEnum,VertexSIdEnum); 
		InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum); 
		InputUpdateFromVectorx(femmodel,RSLg_eustatic,SealevelRSLEustaticEnum,VertexSIdEnum); 

		if (horiz){
			InputUpdateFromVectorx(femmodel,U_north_esa,SealevelUNorthEsaEnum,VertexSIdEnum);	// north motion 
			InputUpdateFromVectorx(femmodel,U_east_esa,SealevelUEastEsaEnum,VertexSIdEnum);		// east motion 
		} /*}}}*/
	}

	if(iscoupler){
		/*transfer sea level back to ice caps:*/
		TransferSealevel(femmodel,SealevelNEsaRateEnum);
		TransferSealevel(femmodel,NGiaRateEnum);
		TransferSealevel(femmodel,SealevelUEsaRateEnum);
		TransferSealevel(femmodel,UGiaRateEnum);

		//reset cumdeltathickness  to 0: 
		InputUpdateFromConstantx(femmodel->inputs,femmodel->elements,0.,SealevelriseCumDeltathicknessEnum);
	}

	/*reset counter to 1:*/
	femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter.

	/*free ressources:{{{*/
	delete RSLg;
	delete RSLg_rate;
	delete RSLg_eustatic;
	delete U_esa;
	delete U_esa_rate;
	delete N_esa;
	delete N_esa_rate;
	delete BPg;

	if(horiz){
		delete U_north_esa;
		delete U_east_esa;
	}
	delete N_gia;
	delete U_gia;
	delete N_gia_rate;
	delete U_gia_rate;
	//delete masks;
	/*}}}*/

} 
/*}}}*/
void dynstr_core(FemModel* femmodel){ /*{{{*/

	/*variables:*/
	Vector<IssmDouble> *bedrock  = NULL; 
	Vector<IssmDouble> *SL  = NULL; 
	Vector<IssmDouble> *steric_rate_g  = NULL; 
	Vector<IssmDouble> *dynamic_rate_g = NULL;
	Vector<IssmDouble> *U_esa_rate= NULL;
	Vector<IssmDouble> *N_esa_rate= NULL;
	Vector<IssmDouble> *U_gia_rate= NULL;
	Vector<IssmDouble> *N_gia_rate= NULL;

	/*parameters: */
	bool isslr=0;
	int  solution_type;
	IssmDouble          dt;
	int  grd=0;
	int  step;
	IssmDouble time;
	IssmDouble cumgmtslr=0;
	IssmDouble cumbslr=0;
	IssmDouble cumgmslr=0;
	IssmDouble gmtslr=0;

	/*Retrieve parameters:*/
	femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
	femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum);
	femmodel->parameters->FindParam(&cumgmtslr,CumGmtslrEnum); 
	femmodel->parameters->FindParam(&cumbslr,CumBslrEnum); 
	femmodel->parameters->FindParam(&step,StepEnum);
	femmodel->parameters->FindParam(&time,TimeEnum);

	/*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
	if(solution_type==SealevelriseSolutionEnum)isslr=1;

	/*Should we be here?:*/
	if(!isslr)return;

	/*Verbose: */
	if(VerboseSolution()) _printf0_("	  computing steric sea level rise\n");

	/*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/
	GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
	GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum);
	GetStericRate(&steric_rate_g,femmodel);
	GetDynamicRate(&dynamic_rate_g,femmodel);
	if(grd){
		GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum);
		GetVectorFromInputsx(&U_gia_rate,femmodel,UGiaRateEnum,VertexSIdEnum);
		GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum);
		GetVectorFromInputsx(&N_gia_rate,femmodel,NGiaRateEnum,VertexSIdEnum);
	}

	/*Save steric rate, and cumulate it through time, cumulate with barystatic sea level rise, and cumulate through time:*/
	gmtslr=steric_rate_g->Max()*dt;
	cumgmtslr+=gmtslr;
	cumgmslr=cumbslr+cumgmtslr;
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,GmtslrEnum,gmtslr,step,time));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,CumGmtslrEnum,cumgmtslr,step,time));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,CumGmslrEnum,cumgmslr,step,time));
	femmodel->parameters->SetParam(cumgmtslr,CumGmtslrEnum);

	/*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate)  * dt + steric_rate + dynamic_rate dt*/
	if(grd){
		SL->AXPY(N_gia_rate,dt);
		SL->AXPY(N_esa_rate,dt);
	}
	SL->AXPY(steric_rate_g,dt);
	SL->AXPY(dynamic_rate_g,dt);

	/*compute new bedrock position: */
	if(grd){
		bedrock->AXPY(U_esa_rate,dt);
		bedrock->AXPY(U_gia_rate,dt);
	}

	/*update element inputs:*/
	InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum);	
	InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum);	

	/*Free ressources:*/	
	delete bedrock;
	delete SL;
	delete steric_rate_g;
	delete dynamic_rate_g;
	if(grd){
		delete U_esa_rate;
		delete U_gia_rate;
		delete N_esa_rate;
		delete N_gia_rate;
	}
}
/*}}}*/

SealevelMasks* sealevelrise_core_masks(FemModel* femmodel) {  /*{{{*/

	if(VerboseSolution()) _printf0_("	  computing sea level masks\n");
	
	/*initialize SealevelMasks structure: */
	SealevelMasks* masks=new SealevelMasks(femmodel->elements->Size());

	/*go through elements and fill the masks: */
	for(Object* & object : femmodel->elements->objects){
		Element*   element=xDynamicCast<Element*>(object);
		element->SetSealevelMasks(masks);
	}

	return masks;
}/*}}}*/
void sealevelrise_core_geometry(FemModel* femmodel) {  /*{{{*/

	/*Geometry core where we compute indices into tables pre computed in the SealevelRiseAnalysis: */

	/*parameters: */
	bool spherical=true;
	IssmDouble *latitude  = NULL;
	IssmDouble *longitude = NULL;
	IssmDouble *radius    = NULL;
	IssmDouble *xx    = NULL;
	IssmDouble *yy    = NULL;
	IssmDouble *zz    = NULL;
	int  horiz;
	bool geometrydone = false;

		
	/*retrieve parameters:*/
	femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
	femmodel->parameters->FindParam(&geometrydone,SealevelriseGeometryDoneEnum);

	if(geometrydone){
		if(VerboseSolution()) _printf0_("	  geometrical offsets have already been computed, skipping \n");
		return; //don't need to run this again.
	}

	/*Verbose: */
	if(VerboseSolution()) _printf0_("	  computing geometrical offsets into precomputed Green tables \n");

	/*first, recover lat,long and radius vectors from vertices: */
	VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); 
	if(horiz) VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); 


	/*Run sealevelrise geometry routine in elements:*/
	for(Object* & object : femmodel->elements->objects){
		Element*   element=xDynamicCast<Element*>(object);
		element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz);
	}

	/*Free ressources:*/
	if(horiz){
		xDelete<IssmDouble>(xx);
		xDelete<IssmDouble>(yy);
		xDelete<IssmDouble>(zz);
	}
	xDelete<IssmDouble>(latitude);
	xDelete<IssmDouble>(longitude);
	xDelete<IssmDouble>(radius);

	/*Record the fact that we ran this module already: */
	femmodel->parameters->SetParam(true,SealevelriseGeometryDoneEnum); 


}/*}}}*/
Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/

	/*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/

	Vector<IssmDouble> *RSLgi    = NULL;
	IssmDouble          RSLgi_oceanaverage   = 0;

	/*parameters: */
	int  gsize;
	IssmDouble oceanarea;
	int        step;
	IssmDouble time;

	/*barystatic contribution:*/
	IssmDouble bslr;
	IssmDouble bslrice;
	IssmDouble* bslrice_partition=NULL;
	IssmDouble bslrhydro;
	IssmDouble* bslrhydro_partition=NULL;
	IssmDouble cumbslr;
	IssmDouble cumbslrice;
	IssmDouble cumbslrhydro;
	IssmDouble* cumbslrice_partition=NULL;
	int npartice;
	IssmDouble* cumbslrhydro_partition=NULL;
	int nparthydro;
	
	if(VerboseSolution()) _printf0_("	  computing bslr components on ice\n");

	/*Figure out size of g-set deflection vector and allocate solution vector: */
	gsize = femmodel->nodes->NumberOfDofs(GsetEnum);

	/*some parameters:*/
	femmodel->parameters->FindParam(&step,StepEnum);
	femmodel->parameters->FindParam(&time,TimeEnum);
	femmodel->parameters->FindParam(&cumbslr,CumBslrEnum);
	femmodel->parameters->FindParam(&cumbslrice,CumBslrIceEnum);
	femmodel->parameters->FindParam(&cumbslrhydro,CumBslrHydroEnum);
	femmodel->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
	femmodel->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
	if(npartice) femmodel->parameters->FindParam(&cumbslrice_partition,&npartice,NULL,CumBslrIcePartitionEnum);
	if(nparthydro) femmodel->parameters->FindParam(&cumbslrhydro_partition,&nparthydro,NULL,CumBslrHydroPartitionEnum);

	/*Initialize:*/
	RSLgi = new Vector<IssmDouble>(gsize);

	/*call the bslr main module: */
	femmodel->SealevelriseEustatic(RSLgi,&oceanarea,&bslr, &bslrice, &bslrhydro, &bslrice_partition, &bslrhydro_partition,masks); //this computes 

	/*we need to average RSLgi over the ocean: RHS term  4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
	RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,masks, oceanarea);

	/*RSLg is the sum of the pure bslr component (term 3) and the contribution from the perturbation to the graviation potential due to the 
	 * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/
	RSLgi->Shift(-bslr-RSLgi_oceanaverage);

	/*save bslr and cumulated bslr value for results:{{{ */
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,BslrEnum,-bslr,step,time));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,BslrIceEnum,-bslrice,step,time));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,BslrHydroEnum,-bslrhydro,step,time));

	//cumulative barystatic contribution: 
	cumbslr=cumbslr-bslr;
	cumbslrice=cumbslrice-bslrice;
	cumbslrhydro=cumbslrhydro-bslrhydro;

	femmodel->parameters->SetParam(cumbslr,CumBslrEnum);
	femmodel->parameters->SetParam(cumbslrice,CumBslrIceEnum);
	femmodel->parameters->SetParam(cumbslrhydro,CumBslrHydroEnum);

	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,CumBslrEnum,cumbslr,step,time));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,CumBslrIceEnum,cumbslrice,step,time));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,CumBslrHydroEnum,cumbslrhydro,step,time));

	//cumulative barystatic contributions by partition:
	if(npartice){
		for(int i=0;i<npartice;i++) cumbslrice_partition[i] -= bslrice_partition[i];
		femmodel->parameters->SetParam(cumbslrice_partition,npartice,1,CumBslrIcePartitionEnum);
		femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,CumBslrIcePartitionEnum,cumbslrice_partition,npartice,1,step,time));
	}

	if(nparthydro){
		for(int i=0;i<nparthydro;i++) cumbslrhydro_partition[i] -= bslrhydro_partition[i];
		femmodel->parameters->SetParam(cumbslrhydro_partition,nparthydro,1,CumBslrHydroPartitionEnum);
		femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,CumBslrHydroPartitionEnum,cumbslrhydro_partition,nparthydro,1,step,time));
	}
	/*}}}*/
	
	/*Assign output pointers and return: */
	*poceanarea=oceanarea;
	return RSLgi;
}/*}}}*/
Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble oceanarea){ /*{{{*/

	/*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
	  non eustatic core of the SLR solution */

	Vector<IssmDouble> *RSLg    = NULL;
	Vector<IssmDouble> *RSLg_old    = NULL;

	Vector<IssmDouble> *RSLgo    = NULL; //ocean convolution of the perturbation to gravity potential.
	Vector<IssmDouble> *RSLgo_rot= NULL; // rotational feedback 
	IssmDouble          RSLgo_oceanaverage = 0;  //average of RSLgo over the ocean.

	/*parameters: */
	int count;
	bool save_results;
	int  gsize;
	bool converged=true;
	bool rotation=true;
	bool verboseconvolution=true;
	int max_nonlinear_iterations;
	IssmDouble           eps_rel;
	IssmDouble           eps_abs;
	IssmDouble			Ixz, Iyz, Izz; 
	
	if(VerboseSolution()) _printf0_("	  converging on ocean components\n");

	/*Recover some parameters: */
	femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
	femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
	femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);

	/*computational flag: */
	femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);

	/*Figure out size of g-set deflection vector and allocate solution vector: */
	gsize = femmodel->nodes->NumberOfDofs(GsetEnum);

	/*Initialize:*/
	RSLg = new Vector<IssmDouble>(gsize);
	RSLg->Assemble();
	RSLg_eustatic->Copy(RSLg);  //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic.

	RSLg_old = new Vector<IssmDouble>(gsize);
	RSLg_old->Assemble();

	count=1;
	converged=false;

	/*Start loop: */
	for(;;){

		//save pointer to old sea level rise
		delete RSLg_old; RSLg_old=RSLg; 

		/*Initialize solution vector: */
		RSLg  = new Vector<IssmDouble>(gsize); RSLg->Assemble();
		RSLgo = new Vector<IssmDouble>(gsize); RSLgo->Assemble();

		/*call the non eustatic module: */
		femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old,  masks, verboseconvolution);

		/*assemble solution vector: */
		RSLgo->Assemble(); 

		if(rotation){

			/*call rotational feedback  module: */
			RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
			femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks); 
			RSLgo_rot->Assemble(); 

			/*save changes in inertia tensor as results: */
			femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz));
			femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz));
			femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,Izz));

			RSLgo->AXPY(RSLgo_rot,1); 
		}

		/*we need to average RSLgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
		RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo,masks, oceanarea);

		/*RSLg is the sum of the eustatic term, and the ocean terms: */
		RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1); 
		RSLg->Shift(-RSLgo_oceanaverage);

		/*convergence criterion:*/
		slrconvergence(&converged,RSLg,RSLg_old,eps_rel,eps_abs);
		
	
		/*free ressources: */
		delete RSLgo;
		delete RSLgo_rot;

		/*Increase count: */
		count++;
		if(converged==true){
			break;
		}
		if(count>=max_nonlinear_iterations){
			_printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); 
			converged=true;
			break;
		}	

		/*some minor verbosing adjustment:*/
		if(count>1)verboseconvolution=false;

	}
	if(VerboseConvergence()) _printf0_("\n              total number of iterations: " << count-1 << "\n");
	
	
	delete RSLg_old;

	return RSLg;
} /*}}}*/
void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks){ /*{{{*/

	Vector<IssmDouble> *U_esa  = NULL; 
	Vector<IssmDouble> *U_north_esa   = NULL; 
	Vector<IssmDouble> *U_east_esa    = NULL; 

	/*parameters: */
	int  gsize;
	bool spherical=true;

	IssmDouble          *latitude   = NULL;
	IssmDouble          *longitude  = NULL;
	IssmDouble          *radius     = NULL;
	IssmDouble          *xx     = NULL;
	IssmDouble          *yy     = NULL;
	IssmDouble          *zz     = NULL;
	int  horiz;
	
	if(VerboseSolution()) _printf0_("	  computing vertical and horizontal geodetic signatures\n");

	/*retrieve some parameters:*/
	femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);

	/*find size of vectors:*/
	gsize      = femmodel->nodes->NumberOfDofs(GsetEnum);

	/*intialize vectors:*/
	U_esa = new Vector<IssmDouble>(gsize);
	if (horiz){
		U_north_esa = new Vector<IssmDouble>(gsize);
		U_east_esa = new Vector<IssmDouble>(gsize);
	}

	/*retrieve geometric information: */
	VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); 
	VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); 

	/*call the elastic main modlule:*/ 
	femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg, masks);

	/*Assign output pointers:*/
	*pU_esa=U_esa;
	if(horiz){
		*pU_east_esa=U_east_esa;
		*pU_north_esa=U_north_esa;
	}

	/*Free ressources: */
	xDelete<IssmDouble>(longitude);
	xDelete<IssmDouble>(latitude);
	xDelete<IssmDouble>(xx);
	xDelete<IssmDouble>(yy);
	xDelete<IssmDouble>(zz);
	xDelete<IssmDouble>(radius);
}
/*}}}*/
void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia, Vector<IssmDouble>** pN_gia,FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/

	/*variables:*/
	Vector<IssmDouble> *U_gia  = NULL; 
	Vector<IssmDouble> *N_gia  = NULL; 

	/*parameters:*/
	int					frequency;
	IssmDouble          dt;
	
	if(VerboseSolution()) _printf0_("	  computing viscous components\n");

	/*retrieve some parameters:*/
	femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);

	/*recover GIA rates:*/
	GetVectorFromInputsx(&U_gia,femmodel,UGiaRateEnum,VertexSIdEnum);
	GetVectorFromInputsx(&N_gia,femmodel,NGiaRateEnum,VertexSIdEnum);

	/*we just loaded rates, that's not what's being asked, scale by time:*/ 
	U_gia->Scale(frequency*dt);
	N_gia->Scale(frequency*dt);

	/*Assign output pointers:*/
	*pU_gia=U_gia;
	*pN_gia=N_gia;

}
/*}}}*/
void sealevelrise_diagnostics(FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/

	/*compute spherical harmonics deg 1 and deg 2 coefficeints:*/

}
/*}}}*/
void GetDynamicRate(Vector<IssmDouble> ** pdynamic_rate_g, FemModel* femmodel){ /*{{{*/

	int dslmodel=-1;
	IssmDouble time;

	/*variables:*/
	Vector<IssmDouble> *dynamic_rate_g  = NULL; 

	/*Update steric rates before retrieving them on Vertex SID set:*/
	femmodel->parameters->FindParam(&dslmodel,DslModelEnum);
	femmodel->parameters->FindParam(&time,TimeEnum);
	if(dslmodel==1){
		TransientInput* transient_input  = femmodel->inputs->GetTransientInput(DslSeaSurfaceHeightChangeAboveGeoidEnum);
		TriaInput* tria_input=transient_input->GetTriaInput(time);
		Input* tria_input_copy=tria_input->copy();
		tria_input_copy->ChangeEnum(DslDynamicRateEnum);
		femmodel->inputs->AddInput(tria_input_copy);
	}
	else if(dslmodel==2){
	
		IssmDouble modelid;
		
		/*Recover modelid:*/
		femmodel->parameters->FindParam(&modelid,DslModelidEnum);
		modelid--; //from matlab. 
		
		/*find the DslSeaSurfaceHeightChangeAboveGeoidEnum dataset of transient inputs:*/
		DatasetInput* dataset_input=femmodel->inputs->GetDatasetInput(DslSeaSurfaceHeightChangeAboveGeoidEnum);
		
		/*Go find the modelid'th transient input:*/
		TriaInput* tria_input=dataset_input->GetTriaInputByOffset(reCast<int, IssmDouble>(modelid));
		
		/*Plug into DslDynamicRate input: */
		Input* tria_input_copy=tria_input->copy();
		tria_input_copy->ChangeEnum(DslDynamicRateEnum);
		femmodel->inputs->AddInput(tria_input_copy);
	}
	else _error_("not implemented yet");

	GetVectorFromInputsx(&dynamic_rate_g,femmodel,DslDynamicRateEnum,VertexSIdEnum);
	*pdynamic_rate_g=dynamic_rate_g;
}
/*}}}*/
void GetStericRate(Vector<IssmDouble> ** psteric_rate_g, FemModel* femmodel){ /*{{{*/

	int dslmodel=-1;
	IssmDouble time;

	/*variables:*/
	Vector<IssmDouble> *steric_rate_g  = NULL; 

	/*Update steric rates before retrieving them on Vertex SID set:*/
	femmodel->parameters->FindParam(&dslmodel,DslModelEnum);
	femmodel->parameters->FindParam(&time,TimeEnum);
	if(dslmodel==1){
		TransientInput* transient_input  = femmodel->inputs->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
		TriaInput* tria_input=transient_input->GetTriaInput(time);
		Input* tria_input_copy=tria_input->copy();
		tria_input_copy->ChangeEnum(DslStericRateEnum);
		femmodel->inputs->AddInput(tria_input_copy);
	}
	else if (dslmodel==2){
		IssmDouble modelid;
		
		/*Recover modelid:*/
		femmodel->parameters->FindParam(&modelid,DslModelidEnum);
		
		modelid--; //from matlab. 
		
		/*find the DslGlobalAverageThermostericSeaLevelChangeEnum dataset of transient inputs:*/
		DatasetInput* dataset_input=femmodel->inputs->GetDatasetInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
		
		/*Go find the modelid'th transient input:*/
		TriaInput* tria_input=dataset_input->GetTriaInputByOffset(reCast<int, IssmDouble>(modelid));
		
		/*Plug into DslStericRate input: */
		Input* tria_input_copy=tria_input->copy();
		tria_input_copy->ChangeEnum(DslStericRateEnum);
		femmodel->inputs->AddInput(tria_input_copy);
	}
	else _error_("not implemented yet");

	GetVectorFromInputsx(&steric_rate_g,femmodel,DslStericRateEnum,VertexSIdEnum);
	*psteric_rate_g=steric_rate_g;
}
/*}}}*/
void SetBottomPressure(FemModel* femmodel){ /*{{{*/

	int dslmodel=-1;
	int type;

	/*If we are running a dsl mme model, we need to update the bottom pressure first: */
	femmodel->parameters->FindParam(&dslmodel,DslModelEnum);
	if(dslmodel==2){
	
		IssmDouble modelid;
		/*figure out the type of DslSeaWaterPressureChangeAtSeaFloor input:*/
		type=femmodel->inputs->GetInputObjectEnum(DslSeaWaterPressureChangeAtSeaFloorEnum);

		if(type==DatasetInputEnum){
		
			/*find the DslSeaWaterPressureChangeAtSeaFloor dataset of transient inputs:*/
			DatasetInput* dataset_input=femmodel->inputs->GetDatasetInput(DslSeaWaterPressureChangeAtSeaFloorEnum);
		
			/*Recover modelid:*/
			femmodel->parameters->FindParam(&modelid,DslModelidEnum);
			modelid--; //from matlab. 
		
			/*Go find the modelid'th transient input:*/
			TransientInput* transient_input=dataset_input->GetTransientInputByOffset(reCast<int, IssmDouble>(modelid));

			/*Plug into DslSeaWaterPressureChangeAtSeaFloor input: */
			Input* transient_input_copy=transient_input->copy();
			transient_input_copy->ChangeEnum(DslSeaWaterPressureChangeAtSeaFloorEnum);
			femmodel->inputs->AddInput(transient_input_copy);
		}
	}
}
/*}}}*/

/*support routines:*/
void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/

	/*forcing being transferred from models to earth: */
	IssmDouble** forcings=NULL;
	IssmDouble*  forcing=NULL; 
	Vector<IssmDouble>* forcingglobal=NULL; 
	int*         nvs=NULL;

	/*transition vectors:*/
	IssmDouble** transitions=NULL;
	int          ntransitions; 
	int*         transitions_m=NULL;
	int*         transitions_n=NULL;
	int          nv;

	/*communicators:*/
	ISSM_MPI_Comm tocomm;
	ISSM_MPI_Comm* fromcomms=NULL;
	ISSM_MPI_Status status;
	int         my_rank;
	int         modelid,earthid;
	int         nummodels;

	/*Recover some parameters: */
	femmodel->parameters->FindParam(&modelid,ModelIdEnum);
	femmodel->parameters->FindParam(&earthid,EarthIdEnum);
	femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
	my_rank=IssmComm::GetRank();

	/*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */
	if(modelid==earthid){
		GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
		if(!parcoms)_error_("TransferForcing error message: could not find IcecapToEarthComm communicator");
		fromcomms=parcoms->GetParameterValue();
	}
	else {
		GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
		if(!parcom)_error_("TransferForcing error message: could not find IcecapToEarthComm communicator");
		tocomm=parcom->GetParameterValue();
	}

	/*For each icecap, retrieve the forcing vector that will be sent to the earth model: */
	if(modelid!=earthid){
		nv=femmodel->vertices->NumberOfVertices();
		GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum);
	}

	/*Send the forcing to the earth model:{{{*/
	if(my_rank==0){
		if(modelid==earthid){
			forcings=xNew<IssmDouble*>(nummodels-1);
			nvs=xNew<int>(nummodels-1);
			for(int i=0;i<earthid;i++){
				ISSM_MPI_Recv(nvs+i, 1, ISSM_MPI_INT, 0,i, fromcomms[i], &status);
				forcings[i]=xNew<IssmDouble>(nvs[i]);
				ISSM_MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status);
			}

		}
		else{
			ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, modelid, tocomm);
			ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, tocomm);
		}
	}
	/*}}}*/

	/*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/
	if(modelid==earthid){

		/*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions. 
		 *First, build the global delta thickness vector in the earth model: */
		nv=femmodel->vertices->NumberOfVertices();
		GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);

		/*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/
		femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);

		if(ntransitions!=earthid)_error_("TransferForcing error message: number of transition vectors is not equal to the number of icecaps!");

		/*Go through all the delta thicknesses coming from each ice cap: */
		if(my_rank==0){
			for(int i=0;i<earthid;i++){

				IssmDouble* forcingfromcap= forcings[i]; //careful, this only exists on rank 0 of the earth model!
				IssmDouble* transition=transitions[i];
				int         M=transitions_m[i];

				/*build index to plug values: */
				int*        index=xNew<int>(M); for(int i=0;i<M;i++)index[i]=reCast<int>(transition[i])-1; //matlab indexing!

				/*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular 
				 * ice cap: */
				forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL);
				xDelete<int>(index);
			}
		}

		/*Assemble vector:*/
		forcingglobal->Assemble();

		/*Plug into elements:*/
		InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum);
	} 
	/*}}}*/

	/*Free ressources:{{{*/
	if(forcings){
		for(int i=0;i<nummodels-1;i++){
			IssmDouble* temp=forcings[i]; xDelete<IssmDouble>(temp);
		}
		xDelete<IssmDouble*>(forcings);
	}
	if(forcing)xDelete<IssmDouble>(forcing);
	if(forcingglobal)delete forcingglobal;
	if(transitions){
		for(int i=0;i<earthid;i++){
			IssmDouble* temp=transitions[i];
			xDelete<IssmDouble>(temp);
		}
		xDelete<IssmDouble*>(transitions);
		xDelete<int>(transitions_m);
		xDelete<int>(transitions_n);
	}
	if(nvs)xDelete<int>(nvs);
	/*}}}*/

} /*}}}*/
void TransferSealevel(FemModel* femmodel,int forcingenum){ /*{{{*/

	/*forcing being transferred from earth to ice caps: */
	IssmDouble*  forcing=NULL; 
	IssmDouble*  forcingglobal=NULL; 

	/*transition vectors:*/
	IssmDouble** transitions=NULL;
	int          ntransitions; 
	int*         transitions_m=NULL;
	int*         transitions_n=NULL;
	int          nv;

	/*communicators:*/
	ISSM_MPI_Comm fromcomm;
	ISSM_MPI_Comm* tocomms=NULL;
	ISSM_MPI_Status status;
	int         my_rank;
	int         modelid,earthid;
	int         nummodels;
	int         numcoms;

	/*Recover some parameters: */
	femmodel->parameters->FindParam(&modelid,ModelIdEnum);
	femmodel->parameters->FindParam(&earthid,EarthIdEnum);
	femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
	my_rank=IssmComm::GetRank();

	/*retrieve the inter communicators that will be used to send data from earth to ice caps:*/
	if(modelid==earthid){
		GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
		if(!parcoms)_error_("TransferSealevel error message: could not find IcecapToEarthComm communicator");
		tocomms=parcoms->GetParameterValue();
		//femmodel->parameters->FindParam((int**)(&tocomms),&numcoms,IcecapToEarthCommEnum);
	}
	else{
		GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
		if(!parcom)_error_("TransferSealevel error message: could not find IcecapToEarthComm communicator");
		fromcomm=parcom->GetParameterValue();
		//femmodel->parameters->FindParam((int*)(&fromcomm), IcecapToEarthCommEnum);
	}

	/*Retrieve sea-level on earth model: */
	if(modelid==earthid){
		nv=femmodel->vertices->NumberOfVertices();
		GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
	}

	/*Send the forcing to the ice caps:{{{*/
	if(my_rank==0){

		if(modelid==earthid){

			/*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */
			femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);

			if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!");

			for(int i=0;i<earthid;i++){
				nv=transitions_m[i];
				forcing=xNew<IssmDouble>(nv);
				IssmDouble* transition=transitions[i];
				for(int j=0;j<nv;j++){
					forcing[j]=forcingglobal[reCast<int>(transition[j])-1];
				}
				ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, i, tocomms[i]);
				ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, i, tocomms[i]);
			}
		}
		else{
			ISSM_MPI_Recv(&nv, 1, ISSM_MPI_INT, 0, modelid, fromcomm, &status);
			forcing=xNew<IssmDouble>(nv);
			ISSM_MPI_Recv(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, fromcomm, &status);
		}
	}
	/*}}}*/

	/*On each ice cap, spread the forcing across cpus, and update the elements dataset accordingly: {{{*/
	if(modelid!=earthid){

		ISSM_MPI_Bcast(&nv,1,ISSM_MPI_INT,0,IssmComm::GetComm());
		if(my_rank!=0)forcing=xNew<IssmDouble>(nv);
		ISSM_MPI_Bcast(forcing,nv,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());

		/*Plug into elements:*/
		InputUpdateFromVectorx(femmodel,forcing,forcingenum,VertexSIdEnum);
	} 
	/*}}}*/

	/*Free ressources:{{{*/
	if(forcingglobal)xDelete<IssmDouble>(forcingglobal);
	if(forcing)xDelete<IssmDouble>(forcing);
	if(transitions){
		for(int i=0;i<ntransitions;i++){
			IssmDouble* temp=transitions[i];
			xDelete<IssmDouble>(temp);
		}
		xDelete<IssmDouble*>(transitions);
		xDelete<int>(transitions_m);
		xDelete<int>(transitions_n);
	}
	/*}}}*/

} /*}}}*/
void EarthMassTransport(FemModel* femmodel){ /*{{{*/

	IssmDouble time,dt;
	Vector<IssmDouble> *oldthickness    = NULL;
	Vector<IssmDouble> *newthickness    = NULL;
	Vector<IssmDouble> *deltathickness    = NULL;
	Vector<IssmDouble> *cumdeltathickness    = NULL;
	int nv;

	if(VerboseSolution()) _printf0_("              computing earth mass transport\n");

	/*This module has to be recoded! We do not grab spc thicknesses from an earth model 
	 * anymore! Thicknesses come from a mass transport module applied to each basin! 
	 * Commeting out for now:*/
	_error_("EarthMassTransport error message: not supported anymore!");

	/*This mass transport module for the Earth is because we might have thickness variations as spcs 
	 * specified in the md.slr class, outside of what we will get from the icecaps. That's why we get t
	 * the thickness variations from SealevelriseSpcthicknessEnum.*/

	/*No mass transport module was called, so we are just going to retrieve the geometry thickness 
	 * at this time step, at prior time step, and plug the difference as deltathickness: */
	/*femmodel->parameters->FindParam(&time,TimeEnum);
	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	nv=femmodel->vertices->NumberOfVertices();

	GetVectorFromInputsx(&newthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum);
	GetVectorFromInputsx(&oldthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum,time-dt);*/

	/*compute deltathickness: */
	/*deltathickness = new Vector<IssmDouble>(nv); 
	newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1); */

	/*plug into elements:*/
//	InputUpdateFromVectorx(femmodel,deltathickness,SurfaceloadIceThicknessChangeEnum,VertexSIdEnum);

	/*add to cumulated delta thickness: */
	/*GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
	cumdeltathickness->AXPY(deltathickness,1); 
	InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);*/

	/*free ressources:*/
	/*delete oldthickness;
	delete newthickness;
	delete deltathickness;
	delete cumdeltathickness;*/

} /*}}}*/
void slrconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/

	bool converged=true;
	IssmDouble ndS,nS; 
	Vector<IssmDouble> *dRSLg    = NULL;

	//compute norm(du) and norm(u) if requested
	dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
	ndS=dRSLg->Norm(NORM_TWO); 

	if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");

	if(!xIsNan<IssmDouble>(eps_rel)){
		nS=RSLg_old->Norm(NORM_TWO);
		if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
	}

	//clean up
	delete dRSLg;

	//print
	if(!xIsNan<IssmDouble>(eps_rel)){
		if((ndS/nS)<eps_rel){
			if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " < " << eps_rel*100 << " %\n");
		}
		else{ 
			if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " > " << eps_rel*100 << " %\n");
			converged=false;
		}
	}
	if(!xIsNan<IssmDouble>(eps_abs)){
		if(ndS<eps_abs){
			if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)" << ndS << " < " << eps_abs << " \n");
		}
		else{ 
			if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)" << ndS << " > " << eps_abs << " \n");
			converged=false;
		}
	}

	/*assign output*/
	*pconverged=converged;

} /*}}}*/
