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

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif


#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"

/*support routines local definitions:{{{*/
void TransferForcing(FemModel* femmodel,int forcingenum);
void TransferSealevel(FemModel* femmodel,int forcingenum);
bool slcconvergence(Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
IssmDouble  SealevelloadsOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble oceanarea);
void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom);
void ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom);
void ivins_deformation_core(FemModel* femmodel);
IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel);
/*}}}*/

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

	SealevelGeometry* slgeom=NULL;

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

	/*Parameters, variables:*/
	bool save_results;

	/*Retrieve parameters:*/
	femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
	
	/*Verbose: */
	if(VerboseSolution()) _printf0_("   computing sea level change\n");

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

	/*run geometry core: */
	slgeom=sealevelchange_geometry(femmodel);
	
	/*any external forcings?:*/
	solidearthexternal_core(femmodel);

	/*Run coupler input transfer:*/
	couplerinput_core(femmodel);

	/*Run geodetic:*/
	grd_core(femmodel,slgeom);

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

	/*Run coupler output transfer: */
	coupleroutput_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,SealevelchangeRequestedOutputsEnum);
		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);}
	}

	/*End profiler*/
	femmodel->profiler->Stop(SLRCORE);
}
/*}}}*/
void              solidearthexternal_core(FemModel* femmodel){ /*{{{*/

	/*variables:*/
	Vector<IssmDouble> *bedrock  = NULL; 
	Vector<IssmDouble> *bedrock_rate = NULL;
	Vector<IssmDouble> *bedrockeast  = NULL; 
	Vector<IssmDouble> *bedrockeast_rate = NULL;
	Vector<IssmDouble> *bedrocknorth  = NULL; 
	Vector<IssmDouble> *bedrocknorth_rate = NULL;
	Vector<IssmDouble> *geoid= NULL; 
	Vector<IssmDouble> *geoid_rate= NULL; 
	int horiz=0;
	int modelid=-1;
	int  isexternal=0;
 
	/*parameters: */
	IssmDouble          dt;

	/*Retrieve parameters:*/
	femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
	femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	femmodel->parameters->FindParam(&isexternal,SolidearthIsExternalEnum); 
	
	/*Early return:*/
	if (!isexternal)return;

	/*Verbose: */
	if(VerboseSolution()) _printf0_("	  computing external solid earth contributions\n");

	/*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/
	GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
	GetVectorFromInputsx(&geoid,femmodel,SealevelEnum,VertexSIdEnum); //In ISSM, Sealevel is absolute.
	if(horiz){
		GetVectorFromInputsx(&bedrockeast,femmodel,BedEastEnum,VertexSIdEnum);
		GetVectorFromInputsx(&bedrocknorth,femmodel,BedNorthEnum,VertexSIdEnum);
	}
	
	GetVectorFromInputsx(&geoid_rate,femmodel,SolidearthExternalGeoidRateEnum,VertexSIdEnum);
	GetVectorFromInputsx(&bedrock_rate,femmodel,SolidearthExternalDisplacementUpRateEnum,VertexSIdEnum);
	if(horiz){
		GetVectorFromInputsx(&bedrockeast_rate,femmodel,SolidearthExternalDisplacementEastRateEnum,VertexSIdEnum);
		GetVectorFromInputsx(&bedrocknorth_rate,femmodel,SolidearthExternalDisplacementNorthRateEnum,VertexSIdEnum);
	}

	/*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate)  * dt + steric_rate + dynamic_rate dt*/
	geoid->AXPY(geoid_rate,dt);
	bedrock->AXPY(bedrock_rate,dt);
	if(horiz){
		bedrockeast->AXPY(bedrockeast_rate,dt);
		bedrocknorth->AXPY(bedrocknorth_rate,dt);
	}

	/*update element inputs:*/
	InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum);	
	InputUpdateFromVectorx(femmodel,geoid,SealevelEnum,VertexSIdEnum);	
	if(horiz){
		InputUpdateFromVectorx(femmodel,bedrockeast,BedEastEnum,VertexSIdEnum);	
		InputUpdateFromVectorx(femmodel,bedrocknorth,BedNorthEnum,VertexSIdEnum);	
	}

	/*Free ressources:*/	
	delete bedrock; delete bedrock_rate;
	delete geoid; delete geoid_rate;
	if(horiz){
		delete bedrockeast; delete bedrockeast_rate;
		delete bedrocknorth; delete bedrocknorth_rate;
	}
}
/*}}}*/
void              couplerinput_core(FemModel* femmodel){  /*{{{*/

	/*Be very careful here, everything is well thought through, do not remove 
	 * without taking big risks:*/
	
	/*parameters:*/
	int  iscoupling;
	int  modelid,earthid;
	int  count,frequency;
	int  horiz;

	/*retrieve more parameters:*/
	femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum);
	femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
	femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);

	if(iscoupling){
		femmodel->parameters->FindParam(&modelid,ModelIdEnum);
		femmodel->parameters->FindParam(&earthid,EarthIdEnum);
	}
	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; 
	}

	/*if we are carrying loads but are not yet computing grd core, accumulate them and skip 
	 * the rest: */
	if (count<frequency){
		count++; 
		femmodel->parameters->SetParam(count,SealevelchangeRunCountEnum); 
		return;
	}

	/*Basins are supposed to accumulate loads and hand them over to the Earth
	  for slr computations every "frequency" time steps. If we are here, we
	  have reached the "frequency"'th time step, and we are going to pick up
	  the old loads, the current loads, and send them to the Earth for slr
	  computations.  So the Earth is never supposed to compute loads. Except,
	  when we are running the Earth as a giant basin (ex: running a
	  Peltier-style GIA model, or running a "Snow Ball" Earth) which only
	  happens when we are not coupled, at which point we are sending these
	  loads to ourselves (hence the convoluted condition).
	  */

	/*transer loads from basins to Earth for grd core computations. :*/
	if(iscoupling){
	
		/*transfer ice thickness change load from basins to earth: */
		TransferForcing(femmodel,DeltaIceThicknessEnum);
		TransferForcing(femmodel,DeltaBottomPressureEnum);
		TransferForcing(femmodel,DeltaTwsEnum);

		/*transfer external forcings back to Earth:*/
		TransferSealevel(femmodel,BedEnum);
		TransferSealevel(femmodel,SealevelEnum);
		if(horiz){
			TransferSealevel(femmodel,BedEastEnum);
			TransferSealevel(femmodel,BedNorthEnum);
		}
	}
	
}; /*}}}*/
void              grd_core(FemModel* femmodel, SealevelGeometry* slgeom) { /*{{{*/

	/*variables:{{{*/
	int nel;
	BarystaticContributions* barycontrib=NULL;
	GenericParam<BarystaticContributions*>* barycontribparam=NULL;
	IssmDouble rotationaxismotionvector[3]={0};
	
	GrdLoads*              loads=NULL;
	Vector<IssmDouble>*    oldsealevelloads=NULL;
	Vector<IssmDouble>*    oceanareas=NULL;
	IssmDouble             oceanarea;
	Vector<IssmDouble>*    subelementoceanareas=NULL;
	IssmDouble             oceanaverage;
	bool                   scaleoceanarea=false;
	IssmDouble             rho_water;

	IssmDouble           eps_rel;
	IssmDouble           eps_abs;
	int                  max_nonlinear_iterations;
	int                  iterations=0;
	int                  step;
	IssmDouble           time; 

	int  modelid,earthid;
	int  horiz;
	int  count,frequency,iscoupling;
	int  grd=0;
	int  grdmodel; 
	int  computesealevel=0;
	bool viscous=false;
	IssmDouble*           sealevelpercpu=NULL;

	/*}}}*/

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

	/*retrieve parameters:{{{*/
	femmodel->parameters->FindParam(&grd,SolidearthSettingsGRDEnum); 
	femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
	femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
	femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum);
	femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
	femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
	femmodel->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
	/*}}}*/

	/*only run if grd was requested, if we are the earth, and we have reached
	 * the necessary number of time steps dictated by :*/
	if(!grd)            return;
	if(count!=frequency)return;
	femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum);
	if(iscoupling){
		femmodel->parameters->FindParam(&modelid,ModelIdEnum);
		femmodel->parameters->FindParam(&earthid,EarthIdEnum);
		if(modelid!=earthid)return;
	}
	/*branch directly to Ivins deformation core if requested:*/
	if(grdmodel==IvinsEnum){
		ivins_deformation_core(femmodel);
		return;
	}

	/*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);
	femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
	/*}}}*/

	/*initialize loads and sea level loads:*/
	femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);

	loads=new GrdLoads(nel,slgeom);
	subelementoceanareas=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
	oceanareas=new Vector<IssmDouble>(nel);
	sealevelpercpu=xNewZeroInit<IssmDouble>(femmodel->vertices->Size());

	if(VerboseSolution()) _printf0_("	  starting  GRD convolutions\n");
	
	/*update viscous RSL:*/
	for(Object* & object : femmodel->elements->objects){
		Element* element = xDynamicCast<Element*>(object);
		element->SealevelchangeUpdateViscousFields();
	}

	/*buildup loads: */
	for(Object* & object : femmodel->elements->objects){
		Element* element = xDynamicCast<Element*>(object);
		element->SealevelchangeBarystaticLoads(loads, barycontrib,slgeom); 
	}

	//broadcast loads 
	loads->BroadcastLoads();

	//compute rotation axis motion:
	RotationAxisMotion(&rotationaxismotionvector[0],femmodel,loads,slgeom);

	/*skip computation of sea level if requested, which means sea level loads should be zeroed */
	if(!computesealevel){
		loads->sealevelloads=xNewZeroInit<IssmDouble>(nel);
		loads->subsealevelloads=xNewZeroInit<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
		goto deformation;
	}

	if(VerboseSolution()) _printf0_("	  converging GRD convolutions\n");
	for(;;){

		oldsealevelloads=loads->vsealevelloads->Duplicate(); loads->vsealevelloads->Copy(oldsealevelloads);

		/*convolve load and sealevel loads on oceans:*/
		for(Object* & object : femmodel->elements->objects){
			Element* element = xDynamicCast<Element*>(object);
			element->SealevelchangeConvolution(sealevelpercpu, loads , rotationaxismotionvector,slgeom);
		}
		
		/*retrieve sea level average  and ocean area:*/
		for(Object* & object : femmodel->elements->objects){
			Element* element = xDynamicCast<Element*>(object);
			element->SealevelchangeOceanAverage(loads, oceanareas, subelementoceanareas, sealevelpercpu, slgeom);
		}

		loads->AssembleSealevelLoads();
	
		/*compute ocean areas:*/
		if(!loads->sealevelloads){ //first time in the loop
			oceanareas->Assemble(); 
			subelementoceanareas->Assemble();
			oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
			if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
		}
	
		//Conserve ocean mass: 
		oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, oceanarea);
		ConserveOceanMass(femmodel,loads,barycontrib->Total()/oceanarea - oceanaverage,slgeom);

		//broadcast sea level loads 
		loads->BroadcastSealevelLoads();

		//compute rotation axis motion:
		RotationAxisMotion(&rotationaxismotionvector[0],femmodel,loads, slgeom);

		//convergence?
		if(slcconvergence(loads->vsealevelloads,oldsealevelloads,eps_rel,eps_abs))break;

		//early return?
		if(iterations>=max_nonlinear_iterations)break;
		iterations++;
	}

	deformation:
	
	if(VerboseSolution()) _printf0_("	  deformation GRD convolutions\n");

	/*convolve loads and sea level loads to get the deformation:*/
	for(Object* & object : femmodel->elements->objects){
		Element* element = xDynamicCast<Element*>(object);
		element->SealevelchangeDeformationConvolution(loads, rotationaxismotionvector,slgeom);
	}

	if(VerboseSolution()) _printf0_("	  updating GRD fields\n");

	/*Update bedrock motion and geoid:*/
	if(computesealevel){
		femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water); //given that we converged, no need to recompute ocean average

		//cumulate barystatic contributions and save to results: 
		barycontrib->Cumulate(femmodel->parameters);
		barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);

		//Avoid double counting barystatic contributions in future time steps: reset all non-cumulative bscl to 0
		for(Object* & object : femmodel->elements->objects){
			Element* element = xDynamicCast<Element*>(object);
			barycontrib->Reset(element->Sid());
		}
	}

	femmodel->inputs->AXPY(1,SealevelGRDEnum,SealevelEnum);
	femmodel->inputs->AXPY(1,BedGRDEnum,BedEnum);
	if(horiz){
		femmodel->inputs->AXPY(1,BedEastGRDEnum,BedEastEnum);
		femmodel->inputs->AXPY(1,BedNorthGRDEnum, BedNorthEnum);
	}

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

	/*variables:*/
	Vector<IssmDouble> *sealevel  = NULL; 
	Vector<IssmDouble> *deltadsl  = NULL; 
	Vector<IssmDouble> *deltastr = NULL;

	/*parameters: */
	int  step;
	int computesealevel=0;
	bool isocean=false;
	IssmDouble time;

	IssmDouble cumgmtslc=0;
	IssmDouble cumbslc=0;
	IssmDouble cumgmslc=0;
	IssmDouble gmtslc=0;

	/*early return if we are not computing sea level, but rather deformation: */
	femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum);
	if (!computesealevel)return;

	/*early return if we have no ocean transport:*/
	femmodel->parameters->FindParam(&isocean,TransientIsoceantransportEnum);
	if(!isocean)return;
	
	/*Verbose: */
	if(VerboseSolution()) _printf0_("	  computing steric and dynamic sea level change\n");

	/*Retrieve sealevel and add steric + dynamic rates:*/
	GetVectorFromInputsx(&sealevel,femmodel,SealevelEnum,VertexSIdEnum);
	GetVectorFromInputsx(&deltadsl,femmodel,DeltaDslEnum,VertexSIdEnum);
	GetVectorFromInputsx(&deltastr,femmodel,DeltaStrEnum,VertexSIdEnum);

	/*compute: sea level change = initial sea level + steric + dynamic*/
	sealevel->AXPY(deltadsl,1);
	sealevel->AXPY(deltastr,1);

	/*cumulate thermal steric rate:*/
	femmodel->parameters->FindParam(&cumgmtslc,CumGmtslcEnum); 
	femmodel->parameters->FindParam(&cumbslc,CumBslcEnum); 

	gmtslc=deltastr->Norm(NORM_TWO);
	cumgmtslc+=gmtslc;
	cumgmslc=cumbslc+cumgmtslc;

	femmodel->parameters->SetParam(cumgmtslc,CumGmtslcEnum);
	femmodel->parameters->SetParam(cumgmslc,CumGmslcEnum);
	
	/*Outputs some metrics:*/
	femmodel->parameters->FindParam(&step,StepEnum);
	femmodel->parameters->FindParam(&time,TimeEnum);

	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,GmtslcEnum,gmtslc,step,time));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,CumGmtslcEnum,cumgmtslc,step,time));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,CumGmslcEnum,cumgmslc,step,time));

	/*update element inputs:*/
	InputUpdateFromVectorx(femmodel,sealevel,SealevelEnum,VertexSIdEnum);	

	/*Free ressources:*/	
	delete sealevel;
	delete deltadsl;
	delete deltastr;
}
/*}}}*/
void              coupleroutput_core(FemModel* femmodel){  /*{{{*/
	
	/*parameters:*/
	int iscoupling;
	int horiz=0;

	/*retrieve more parameters:*/
	femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum);
	femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
		
	if(iscoupling){
		/*transfer sea level back to ice caps:*/
		TransferSealevel(femmodel,SealevelEnum);
		TransferSealevel(femmodel,BedEnum);
		if(horiz){
			TransferSealevel(femmodel,BedNorthEnum);
			TransferSealevel(femmodel,BedEastEnum);
		}
	}
}; /*}}}*/
void              ivins_deformation_core(FemModel* femmodel){ /*{{{*/

	int  gsize;
	Vector<IssmDouble> *bedup  = NULL; 
	Vector<IssmDouble> *beduprate= NULL; 
	IssmDouble          *xx     = NULL;
	IssmDouble          *yy     = NULL;
	
	if(VerboseSolution()) _printf0_("	  computing vertical deformation using Ivins model. \n");

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

	/*Find the litho material to be used by all the elements:*/
	Matlitho* matlitho=NULL;
	for (Object* & object: femmodel->materials->objects){
		Material* material=xDynamicCast<Material*>(object);
		if(material->ObjectEnum()==MatlithoEnum){
			matlitho=xDynamicCast<Matlitho*>(material);
			break;
		}
	}

	/*initialize vectors:*/
	bedup = new Vector<IssmDouble>(gsize);
	beduprate = new Vector<IssmDouble>(gsize);
	
	/*retrieve geometric information: */
	VertexCoordinatesx(&xx,&yy,NULL,femmodel->vertices); 

	/*Go through elements, and add contribution from each element to the deflection vector wg:*/
	for(Object* & object : femmodel->elements->objects){
		Element* element = xDynamicCast<Element*>(object);
		element->GiaDeflection(bedup,beduprate, matlitho, xx,yy);
	}

	/*Assemble parallel vector:*/
	beduprate->Assemble();
	bedup->Assemble();

	/*Save results:*/
	InputUpdateFromVectorx(femmodel,bedup,BedGRDEnum,VertexSIdEnum);
	femmodel->inputs->AXPY(1,BedGRDEnum,BedEnum);

	/*Free ressources: */
	xDelete<IssmDouble>(xx);
	xDelete<IssmDouble>(yy);
	delete beduprate;
	delete bedup;
}
/*}}}*/
void              sealevelchange_initialgeometry(FemModel* femmodel) {  /*{{{*/

	/*Geometry core where we compute geometrical kernels and weights:*/

	/*parameters: */
	IssmDouble *xxe    = NULL;
	IssmDouble *yye    = NULL;
	IssmDouble *zze    = NULL;
	IssmDouble* areae  = NULL;
	int  nel;
	int  grdmodel=0;

	/*retrieve parameters:*/
	femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
	nel=femmodel->elements->NumberOfElements();
		
	/*early return?:*/
	if(grdmodel==IvinsEnum) return;

	/*Verbose: */
	if(VerboseSolution()) _printf0_("	  computing initial sea level geometrical kernels and weights.\n");

	/*recover x,y,z and areas from elements: */
	ElementCoordinatesx(&xxe,&yye,&zze,&areae,femmodel->elements);

	/*Run sealevel geometry routine in elements:*/
	for(Object* & object : femmodel->elements->objects){
		Element*   element=xDynamicCast<Element*>(object);
		element->SealevelchangeGeometryInitial(xxe,yye,zze,areae);
	}

	femmodel->parameters->AddObject(new DoubleVecParam(XxeEnum,xxe,nel));
	femmodel->parameters->AddObject(new DoubleVecParam(YyeEnum,yye,nel));
	femmodel->parameters->AddObject(new DoubleVecParam(ZzeEnum,zze,nel));
	femmodel->parameters->AddObject(new DoubleVecParam(AreaeEnum,areae,nel));


	#ifdef _ISSM_DEBUG_
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,XxeEnum,xxe,nel,1,1,1));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,YyeEnum,yye,nel,1,1,1));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,ZzeEnum,zze,nel,1,1,1));
	femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,AreaeEnum,areae,nel,1,1,1));
	#endif

	return;


}/*}}}*/
SealevelGeometry* sealevelchange_geometry(FemModel* femmodel) {  /*{{{*/

	/*Geometry core where we compute updates to the Green function kernels and weights, dependent 
	 * on the evolution of levelsets: */

	/*parameters: */
	IssmDouble *xxe    = NULL;
	IssmDouble *yye    = NULL;
	IssmDouble *zze    = NULL;
	IssmDouble* areae  = NULL;

	int nel;
	int  grdmodel=0;
	SealevelGeometry* slgeom=NULL;

	/*early return?:*/
	femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
	if(grdmodel==IvinsEnum) return NULL;

	/*retrieve parameters:*/
	femmodel->parameters->FindParam(&xxe,&nel,XxeEnum);
	femmodel->parameters->FindParam(&yye,&nel,YyeEnum);
	femmodel->parameters->FindParam(&zze,&nel,ZzeEnum);
	femmodel->parameters->FindParam(&areae,&nel,AreaeEnum);
	
	/*initialize SealevelMasks structure: */
	slgeom=new SealevelGeometry(femmodel->elements->Size(),femmodel->vertices->Size());
	
	/*Verbose: */
	if(VerboseSolution()) _printf0_("	  computing sea level geometrical kernel and weight updates.\n");
	
	
	/*Run sealevel geometry routine for elements with full loading:*/
	for(Object* & object : femmodel->elements->objects){
		Element*   element=xDynamicCast<Element*>(object);
		element->SealevelchangeGeometryCentroidLoads(slgeom,xxe,yye,zze,areae);
	}
	
	/*Initialize fractional loading mapping: */
	slgeom->InitializeMappingsAndBarycentres();

	
	/*Run sealevel geometry routine for elements with fractional loading:*/
	for(Object* & object : femmodel->elements->objects){
		Element*   element=xDynamicCast<Element*>(object);
		element->SealevelchangeGeometrySubElementLoads(slgeom,areae);
	}

	/*Assemble barycentres of fraction loading elements:*/
	slgeom->Assemble();

	/*Create fractional green function kernels: */
	for(Object* & object : femmodel->elements->objects){
		Element*   element=xDynamicCast<Element*>(object);
		element->SealevelchangeGeometryFractionKernel(slgeom);
	}


	femmodel->parameters->AddObject(new DoubleVecParam(XxeEnum,xxe,nel));
	femmodel->parameters->AddObject(new DoubleVecParam(YyeEnum,yye,nel));
	femmodel->parameters->AddObject(new DoubleVecParam(ZzeEnum,zze,nel));
	femmodel->parameters->AddObject(new DoubleVecParam(AreaeEnum,areae,nel));

	return slgeom;

}/*}}}*/

/*subroutines:*/
bool slcconvergence(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 (RSL_old=" << RSLg_old->Norm(NORM_TWO) << " RSL=" << RSLg->Norm(NORM_TWO) << ")");
	}

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

	//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*/
	return converged;

} /*}}}*/
IssmDouble  SealevelloadsOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* suboceanareas, IssmDouble oceanarea){ /*{{{*/

	IssmDouble sealevelloadsaverage;	
	IssmDouble subsealevelloadsaverage;	

	Vector<IssmDouble>* vsealevelloadsvolume=loads->vsealevelloads->Duplicate();
	Vector<IssmDouble>* vsubsealevelloadsvolume=loads->vsubsealevelloads->Duplicate();

	vsealevelloadsvolume->PointwiseMult(loads->vsealevelloads,oceanareas);
	vsubsealevelloadsvolume->PointwiseMult(loads->vsubsealevelloads,suboceanareas);
	
	vsealevelloadsvolume->Sum(&sealevelloadsaverage);
	vsubsealevelloadsvolume->Sum(&subsealevelloadsaverage);
	delete vsealevelloadsvolume; 
	delete vsubsealevelloadsvolume; 
	
	return (sealevelloadsaverage+subsealevelloadsaverage)/oceanarea;
} /*}}}*/
void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom){ /*{{{*/

	IssmDouble  moi_list[3]={0,0,0};
	IssmDouble  moi_list_sub[3]={0,0,0};
	IssmDouble  moi_list_cpu[3]={0,0,0};
	IssmDouble*	tide_love_h  = NULL;
	IssmDouble*	tide_love_k  = NULL;
	IssmDouble*	load_love_k  = NULL;
	IssmDouble  tide_love_k2secular;
	IssmDouble  moi_e, moi_p;
	IssmDouble	m1, m2, m3;
	bool rotation=false;
	
	/*early return?:*/
	femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
	if(!rotation)return;

	/*retrieve parameters: */
	femmodel->parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
	femmodel->parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
	femmodel->parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
	femmodel->parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
	femmodel->parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
	femmodel->parameters->FindParam(&moi_p,RotationalPolarMoiEnum);

	for(Object* & object : femmodel->elements->objects){
		Element* element = xDynamicCast<Element*>(object);

		element->SealevelchangeMomentOfInertiaCentroid(&moi_list[0],loads,slgeom);

		element->SealevelchangeMomentOfInertiaSubElement(&moi_list_sub[0],loads, slgeom);

		moi_list_cpu[0] += moi_list[0]+moi_list_sub[0];
		moi_list_cpu[1] += moi_list[1]+moi_list_sub[1];
		moi_list_cpu[2] += moi_list[2]+moi_list_sub[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());

	/*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

	/*Assign output pointers:*/
	m[0]=m1;
	m[1]=m2;
	m[2]=m3;
} /*}}}*/
void        ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/

	/*Shift sealevel loads by ocean average, only on ocean! :*/
	for(Object* & object : femmodel->elements->objects){
		Element* element = xDynamicCast<Element*>(object);
		element->SealevelchangeShift(loads, offset,slgeom);
	}
	loads->AssembleSealevelLoads();

} /*}}}*/
IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel){ /*{{{*/

	int* indices=xNew<int>(nel);
	for(int i=0;i<nel;i++)indices[i]=i;
	
	Vector<IssmDouble>* vloadcopy=new Vector<IssmDouble>(nel);
	IssmDouble* loadcopy=xNew<IssmDouble>(nel);
	
	vloadcopy->SetValues(nel,indices,load,INS_VAL);
	vloadcopy->Assemble();


	if(subload){
		for (int i=0;i<femmodel->elements->Size();i++){
			if (slgeom->issubelement[loadtype][i]){
				int se= slgeom->subelementmapping[loadtype][i];
				IssmDouble subloadi=subload[se];
				Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
				vloadcopy->SetValue(element->Sid(),subloadi,ADD_VAL);
			}
		}
	}
	vloadcopy->Assemble();
	loadcopy=vloadcopy->ToMPISerial();

	return loadcopy;

} /*}}}*/

/*Coupling 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;
	int          existforcing=0;         

	/*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();
		existforcing=reCast<int>(femmodel->inputs->Exist(forcingenum));
		if(existforcing)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(&existforcing, 1, ISSM_MPI_INT, 0,i, fromcomms[i], &status);
				if(existforcing){
					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{
					forcings[i]=NULL;
				}
			}

		}
		else{
			ISSM_MPI_Send(&existforcing, 1, ISSM_MPI_INT, 0, modelid, tocomm);
			if(existforcing){
				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,SealevelchangeTransitionsEnum);

		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!
				if(forcingfromcap){
					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]; 
			if(temp)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,SealevelchangeTransitionsEnum);

			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);
	}
	/*}}}*/

} /*}}}*/
