/*!\file Element.cpp
 * \brief: implementation of the Element object
 */
/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "../classes.h"
#include "../../shared/shared.h"
#include "../../modules/SurfaceMassBalancex/SurfaceMassBalancex.h"
#include "../Inputs/BoolInput.h"
#include "../Inputs/TransientInput.h"
#include "../Inputs/ElementInput.h"
#include "../Inputs/PentaInput.h"
#include "../Inputs/DatasetInput.h"
#include "../Inputs/ControlInput.h"
#include "../Inputs/ArrayInput.h"
/*}}}*/
#define MAXVERTICES 6 /*Maximum number of vertices per element, currently Penta, to avoid dynamic mem allocation*/

#ifdef _HAVE_SEMIC_
/* SEMIC prototype {{{*/
extern "C" void run_semic_(IssmDouble *sf_in, IssmDouble *rf_in, IssmDouble *swd_in, IssmDouble *lwd_in, IssmDouble *wind_in, IssmDouble *sp_in, IssmDouble *rhoa_in,
			IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *saccu_out, IssmDouble *smelt_out);
#endif
// _HAVE_SEMIC_
/*}}}*/
/*Constructors/destructor/copy*/
Element::Element(){/*{{{*/
	this->id  = -1;
	this->sid = -1;
	this->lid = -1;
	this->inputs    = NULL;
	this->nodes      = NULL;
	this->vertices   = NULL;
	this->material   = NULL;
	this->parameters = NULL;
	this->element_type_list=NULL;
}/*}}}*/
Element::~Element(){/*{{{*/
	xDelete<int>(element_type_list);
}
/*}}}*/

/*Other*/
bool       Element::AnyFSet(){/*{{{*/

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = this->GetNumberOfNodes();

	for(int i=0;i<numnodes;i++){
		if(nodes[i]->fsize) return true;
	}
	return false;
}/*}}}*/
void       Element::ComputeLambdaS(){/*{{{*/

	/*Intermediaries*/
	IssmDouble vx,vy,vz,vmag;
	IssmDouble dvx[3],dvy[3],dvz[3],dvmag[3];
	IssmDouble eps[3][3],epseff,epsprime;
	IssmDouble lambdas[MAXVERTICES];
	int         dim;
	IssmDouble *xyz_list = NULL;

	/*Retrieve all inputs we will be needing: */
	this->GetVerticesCoordinates(&xyz_list);
	parameters->FindParam(&dim,DomainDimensionEnum);
	Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input);
	Input* vz_input=NULL;
	if(dim==3){vz_input=this->GetInput(VzEnum); _assert_(vz_input);}

	/*Allocate arrays*/
	const int NUM_VERTICES = this->GetNumberOfVertices();

	/* Start looping on the number of vertices: */
	Gauss* gauss=this->NewGauss();
	for (int iv=0;iv<NUM_VERTICES;iv++){
		gauss->GaussVertex(iv);

		/*Get velocity derivatives in all directions*/
		_assert_(dim>1);
		_assert_(vx_input);
		vx_input->GetInputValue(&vx,gauss);
		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
		_assert_(vy_input);
		vy_input->GetInputValue(&vy,gauss);
		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
		if(dim==3){
			_assert_(vz_input);
			vz_input->GetInputValue(&vz,gauss);
			vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
		}
		else{
			vz = 0.;
			dvz[0] = 0.; dvz[1] = 0.; dvz[2] = 0.;
			dvx[2]= 0.;
			dvy[2]= 0.;
		}
		/*Calculate velocity magnitude and its derivative*/
		vmag = sqrt(vx*vx+vy*vy+vz*vz);
		if(vmag<1e-12){
			vmag=1e-12;
			dvmag[0]=0;
			dvmag[1]=0;
			dvmag[2]=0;
		}
		else{
			dvmag[0]=1./(2*sqrt(vmag))*(2*vx*dvx[0]+2*vy*dvy[0]+2*vz*dvz[0]);
			dvmag[1]=1./(2*sqrt(vmag))*(2*vx*dvx[1]+2*vy*dvy[1]+2*vz*dvz[1]);
			dvmag[2]=1./(2*sqrt(vmag))*(2*vx*dvx[2]+2*vy*dvy[2]+2*vz*dvz[2]);
		}
		/*Build strain rate tensor*/
		eps[0][0] = dvx[0];             eps[0][1] = .5*(dvx[1]+dvy[0]); eps[0][2] = .5*(dvx[2]+dvz[0]);
		eps[1][0] = .5*(dvx[1]+dvy[0]); eps[1][1] = dvy[1];             eps[1][2] = .5*(dvy[2]+dvz[1]);
		eps[2][0] = .5*(dvx[2]+dvz[0]); eps[2][1] = .5*(dvy[2]+dvz[1]); eps[2][2] = dvz[2];

		/*effective strain rate*/
		epseff = 0.;
		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
		epseff = sqrt(eps[0][0]*eps[0][0] + eps[1][1]*eps[1][1] + eps[0][1]*eps[0][1] +  eps[0][2]*eps[0][2] + eps[1][2]*eps[1][2] + eps[0][0]*eps[1][1]);

		EstarStrainrateQuantities(&epsprime,vx,vy,vz,vmag,&dvx[0],&dvy[0],&dvz[0],&dvmag[0]);
		lambdas[iv]=EstarLambdaS(epseff,epsprime);
	}

	/*Add Stress tensor components into inputs*/
	this->AddInput(LambdaSEnum,&lambdas[0],P1Enum);

	/*Clean up and return*/
	delete gauss;
	xDelete<IssmDouble>(xyz_list);
}/*}}}*/
void       Element::ComputeNewDamage(){/*{{{*/

	IssmDouble *xyz_list=NULL;
	IssmDouble  eps_xx,eps_xy,eps_yy,eps_xz,eps_yz,eps_zz,eps_eff;
	IssmDouble  epsmin=1.e-27;
	IssmDouble  eps_0,kappa,sigma_0,B,D,n,envelopeD;
	int         dim,counter=0;
	IssmDouble  k1,k2,threshold=1.e-12;

	/* Retrieve parameters */
	this->GetVerticesCoordinates(&xyz_list);
	this->ComputeStrainRate();
	parameters->FindParam(&dim,DomainDimensionEnum);
	parameters->FindParam(&kappa,DamageKappaEnum);
	parameters->FindParam(&sigma_0,DamageStressThresholdEnum);

	/* Retrieve inputs */
	Input* eps_xx_input=this->GetInput(StrainRatexxEnum); _assert_(eps_xx_input);
	Input* eps_yy_input=this->GetInput(StrainRateyyEnum); _assert_(eps_yy_input);
	Input* eps_xy_input=this->GetInput(StrainRatexyEnum); _assert_(eps_xy_input);
	Input* eps_xz_input=NULL;
	Input* eps_yz_input=NULL;
	Input* eps_zz_input=NULL;
	if(dim==3){
		eps_xz_input=this->GetInput(StrainRatexzEnum); _assert_(eps_xz_input);
		eps_yz_input=this->GetInput(StrainRateyzEnum); _assert_(eps_yz_input);
		eps_zz_input=this->GetInput(StrainRatezzEnum); _assert_(eps_zz_input);
	}

	/* Fetch number of nodes and allocate output*/
	int numnodes = this->GetNumberOfNodes();
	IssmDouble* newD = xNew<IssmDouble>(numnodes);

	/* Retrieve domain-dependent inputs */
	Input* n_input=this->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
	Input* damage_input = NULL;
	Input* B_input = NULL;
	int domaintype;
	parameters->FindParam(&domaintype,DomainTypeEnum);
	if(domaintype==Domain2DhorizontalEnum){
		damage_input = this->GetInput(DamageDbarOldEnum);  _assert_(damage_input);
		B_input=this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input);
	}
	else{
		damage_input = this->GetInput(DamageDOldEnum);   _assert_(damage_input);
		B_input=this->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
	}

	/* Start looping on the number of nodes: */
	Gauss* gauss=this->NewGauss();
	for (int i=0;i<numnodes;i++){
		gauss->GaussNode(this->GetElementType(),i);

		eps_xx_input->GetInputValue(&eps_xx,gauss);
		eps_yy_input->GetInputValue(&eps_yy,gauss);
		eps_xy_input->GetInputValue(&eps_xy,gauss);
		if(dim==3){
			eps_xz_input->GetInputValue(&eps_xz,gauss);
			eps_yz_input->GetInputValue(&eps_yz,gauss);
			eps_zz_input->GetInputValue(&eps_zz,gauss);
		}
		else{eps_xz=0; eps_yz=0; eps_zz=0;}

		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
		eps_eff=sqrt(eps_xx*eps_xx+eps_yy*eps_yy+eps_xy*eps_xy+eps_xz*eps_xz+eps_yz*eps_yz+eps_xx*eps_yy+epsmin*epsmin);

		B_input->GetInputValue(&B,gauss);
		n_input->GetInputValue(&n,gauss);
		damage_input->GetInputValue(&D,gauss);

		/* Compute threshold strain rate from threshold stress */
		eps_0=pow(sigma_0/B,n);

		if(eps_eff>eps_0){
			/* Compute damage on envelope curve for existing level of effective strain rate */
			envelopeD=1.-pow(eps_0/eps_eff,1./n)*exp(-(eps_eff-eps_0)/(eps_0*(kappa-1.)));

			if(envelopeD>D){
				newD[i]=envelopeD;
			}
			else newD[i]=D;
		}
		else newD[i]=D;
	}

	/* Add new damage input to DamageEnum and NewDamageEnum */
	this->AddInput(NewDamageEnum,newD,P1DGEnum);
	if(domaintype==Domain2DhorizontalEnum){
		this->AddInput(DamageDbarEnum,newD,this->GetElementType());
	}
	else{
		this->AddInput(DamageDEnum,newD,this->GetElementType());
	}

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(newD);
	delete gauss;

}/*}}}*/
void       Element::ComputeStrainRate(){/*{{{*/

	int         dim;
	IssmDouble *xyz_list = NULL;
	IssmDouble  epsilon[6];

	/*Retrieve all inputs we will be needing: */
	this->GetVerticesCoordinates(&xyz_list);
	parameters->FindParam(&dim,DomainDimensionEnum);
	Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input);
	Input* vz_input=NULL;
	if(dim==3){vz_input=this->GetInput(VzEnum); _assert_(vz_input);}

	/*Allocate arrays*/
	const int NUM_VERTICES = this->GetNumberOfVertices();
	_assert_(NUM_VERTICES<=MAXVERTICES);

	IssmDouble eps_xx[MAXVERTICES];
	IssmDouble eps_yy[MAXVERTICES];
	IssmDouble eps_zz[MAXVERTICES];
	IssmDouble eps_xy[MAXVERTICES];
	IssmDouble eps_xz[MAXVERTICES];
	IssmDouble eps_yz[MAXVERTICES];
	IssmDouble eps_ef[MAXVERTICES];

	/* Start looping on the number of vertices: */
	Gauss* gauss=this->NewGauss();
	for (int iv=0;iv<NUM_VERTICES;iv++){
		gauss->GaussVertex(iv);

		/*Compute strain rate viscosity and pressure: */
		if(dim==2)
		 this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
		else
		 this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);

		if(dim==2){
			/* epsilon=[exx,eyy,exy];*/
			eps_xx[iv]=epsilon[0];
			eps_yy[iv]=epsilon[1];
			eps_xy[iv]=epsilon[2];
			/* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
			eps_ef[iv] = 1./sqrt(2.)*sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + 2.*epsilon[2]*epsilon[2]);
		}
		else{
			/*epsilon=[exx eyy ezz exy exz eyz]*/
			eps_xx[iv]=epsilon[0];
			eps_yy[iv]=epsilon[1];
			eps_zz[iv]=epsilon[2];
			eps_xy[iv]=epsilon[3];
			eps_xz[iv]=epsilon[4];
			eps_yz[iv]=epsilon[5];
			/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
			eps_ef[iv] = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[3]*epsilon[3] +  epsilon[4]*epsilon[4] + epsilon[5]*epsilon[5] + epsilon[0]*epsilon[1]);
		}
	}

	/*Add Stress tensor components into inputs*/
	this->AddInput(StrainRatexxEnum,&eps_xx[0],P1Enum);
	this->AddInput(StrainRatexyEnum,&eps_xy[0],P1Enum);
	this->AddInput(StrainRatexzEnum,&eps_xz[0],P1Enum);
	this->AddInput(StrainRateyyEnum,&eps_yy[0],P1Enum);
	this->AddInput(StrainRateyzEnum,&eps_yz[0],P1Enum);
	this->AddInput(StrainRatezzEnum,&eps_zz[0],P1Enum);
	this->AddInput(StrainRateeffectiveEnum,&eps_ef[0],P1Enum);

	/*Clean up and return*/
	delete gauss;
	xDelete<IssmDouble>(xyz_list);
}
/*}}}*/
void       Element::CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/

	int         i,counter;
	int         numdofs   = 0;
	IssmDouble  norm;
	IssmDouble *transform = NULL;
	IssmDouble  coord_system[3][3];

	/*Some checks in debugging mode*/
	_assert_(numnodes && nodes_list);

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Allocate and initialize transform matrix*/
	transform=xNew<IssmDouble>(numdofs*numdofs);
	for(i=0;i<numdofs*numdofs;i++) transform[i]=0.0;

	/*Create transform matrix for all nodes (x,y for 2d and x,y,z for 3d). It is a block matrix
	 *for 3 nodes:

	 *     | T1 0  0 |
	 * Q = | 0  T2 0 |
	 *     | 0  0  T3|
	 *
	 * Where T1 is the transform matrix for node 1. It is a simple copy of the coordinate system
	 * associated to this node*/
	counter=0;
	for(i=0;i<numnodes;i++){
		nodes_list[i]->GetCoordinateSystem(&coord_system[0][0]);
		switch(cs_array[i]){
			case PressureEnum:
				/*DO NOT change anything*/
				transform[(numdofs)*(counter) + counter] = 1.;
				counter+=1;
				break;
			case XYEnum:
				/*We remove the z component, we need to renormalize x and y: x=[x1 x2 0] y=[-x2 x1 0]*/
				norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]); _assert_(norm>1.e-4);
				transform[(numdofs)*(counter+0) + counter+0] =   coord_system[0][0]/norm;
				transform[(numdofs)*(counter+0) + counter+1] = - coord_system[1][0]/norm;
				transform[(numdofs)*(counter+1) + counter+0] =   coord_system[1][0]/norm;
				transform[(numdofs)*(counter+1) + counter+1] =   coord_system[0][0]/norm;
				counter+=2;
				break;
			case XYZEnum:
				/*The 3 coordinates are changed (x,y,z)*/
				transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0];
				transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1];
				transform[(numdofs)*(counter+0) + counter+2] = coord_system[0][2];
				transform[(numdofs)*(counter+1) + counter+0] = coord_system[1][0];
				transform[(numdofs)*(counter+1) + counter+1] = coord_system[1][1];
				transform[(numdofs)*(counter+1) + counter+2] = coord_system[1][2];
				transform[(numdofs)*(counter+2) + counter+0] = coord_system[2][0];
				transform[(numdofs)*(counter+2) + counter+1] = coord_system[2][1];
				transform[(numdofs)*(counter+2) + counter+2] = coord_system[2][2];
				counter+=3;
				break;
			default:
				_error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Assign output pointer*/
	*ptransform=transform;
}
/*}}}*/
void       Element::DeepEcho(void){/*{{{*/

	_printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
	_printf_("   id : "<<this->id <<"\n");
	_printf_("   sid: "<<this->sid<<"\n");
	_printf_("   lid: "<<this->lid<<"\n");
	if(vertices){
		const int NUM_VERTICES = this->GetNumberOfVertices();
		for(int i=0;i<NUM_VERTICES;i++) vertices[i]->Echo();
	}
	else _printf_("vertices = NULL\n");

	if(nodes){
		int numnodes = this->GetNumberOfNodes();
		for(int i=0;i<numnodes;i++) nodes[i]->DeepEcho();
	}
	else _printf_("nodes = NULL\n");

	if (material) material->DeepEcho();
	else _printf_("material = NULL\n");

	_printf_("   parameters\n");
	if (parameters) parameters->DeepEcho();
	else _printf_("parameters = NULL\n");

	_printf_("   inputs\n");
	if(inputs) inputs->DeepEcho();
	else _printf_("inputs=NULL\n");

	return;
}
/*}}}*/
void       Element::DeleteMaterials(void){/*{{{*/
	delete this->material;
}/*}}}*/
void       Element::Delta18oParameterization(void){/*{{{*/

	/*Are we on the base? If not, return*/
	if(!IsOnBase()) return;

	const int NUM_VERTICES = this->GetNumberOfVertices();
	const int NUM_VERTICES_MONTHS_PER_YEAR	= NUM_VERTICES * 12;

	int        i;
	int*        vertexlids=xNew<int>(NUM_VERTICES);
	IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* TemperaturesLgm=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);

	/*Recover parameters*/
	IssmDouble time,yts,finaltime;
	this->parameters->FindParam(&time,TimeEnum);
	this->parameters->FindParam(&yts,ConstantsYtsEnum);
	this->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
	this->GetVerticesLidList(vertexlids);
	IssmDouble time_yr=floor(time/yts)*yts;

	/*Recover present day temperature and precipitation*/
	DatasetInput* dinput1=this->GetDatasetInput(SmbTemperaturesPresentdayEnum);   _assert_(dinput1);
	DatasetInput* dinput2=this->GetDatasetInput(SmbTemperaturesLgmEnum);          _assert_(dinput2);
	DatasetInput* dinput3=this->GetDatasetInput(SmbPrecipitationsPresentdayEnum); _assert_(dinput3);

	/*loop over vertices: */
	Gauss* gauss=this->NewGauss();
	for(int month=0;month<12;month++){
		for(int iv=0;iv<NUM_VERTICES;iv++){
			gauss->GaussVertex(iv);
			dinput1->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month);
			dinput2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss,month);
			dinput3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month);

			PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts;
		}
	}

	/*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/
	IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime;
	IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime;
	this->parameters->FindParam(&Delta18oPresent,SmbDelta18oEnum,finaltime);
	this->parameters->FindParam(&Delta18oLgm,SmbDelta18oEnum,(finaltime-(21000*yts)));
	this->parameters->FindParam(&Delta18oTime,SmbDelta18oEnum,time);
	this->parameters->FindParam(&Delta18oSurfacePresent,SmbDelta18oSurfaceEnum,finaltime);
	this->parameters->FindParam(&Delta18oSurfaceLgm,SmbDelta18oSurfaceEnum,(finaltime-(21000*yts)));
	this->parameters->FindParam(&Delta18oSurfaceTime,SmbDelta18oSurfaceEnum,time);

	/*Compute the temperature and precipitation*/
	for(int iv=0;iv<NUM_VERTICES;iv++){
		ComputeDelta18oTemperaturePrecipitation(Delta18oSurfacePresent, Delta18oSurfaceLgm, Delta18oSurfaceTime,
					Delta18oPresent, Delta18oLgm, Delta18oTime,
					&PrecipitationsPresentday[iv*12],
					&TemperaturesLgm[iv*12], &TemperaturesPresentday[iv*12],
					&monthlytemperatures[iv*12], &monthlyprec[iv*12]);
	}

	/*Update inputs*/
	for (int imonth=0;imonth<12;imonth++) {
		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
		switch(this->ObjectEnum()){
			case TriaEnum: this->inputs->SetTriaDatasetInput(SmbMonthlytemperaturesEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			case PentaEnum: this->inputs->SetPentaDatasetInput(SmbMonthlytemperaturesEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			default: _error_("Not implemented yet");
		}
		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
		switch(this->ObjectEnum()){
			case TriaEnum: this->inputs->SetTriaDatasetInput(SmbPrecipitationEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			case PentaEnum: this->inputs->SetPentaDatasetInput(SmbPrecipitationEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			default: _error_("Not implemented yet");
		}
	}

	switch(this->ObjectEnum()){
		case TriaEnum: break;
		case PentaEnum:
		case TetraEnum:
         this->DatasetInputExtrude(SmbMonthlytemperaturesEnum,-1);
         this->DatasetInputExtrude(SmbPrecipitationEnum,-1);
         break;
		default: _error_("Not implemented yet");
	}

	/*clean-up*/
	delete gauss;
	xDelete<IssmDouble>(monthlytemperatures);
	xDelete<IssmDouble>(monthlyprec);
	xDelete<IssmDouble>(TemperaturesPresentday);
	xDelete<IssmDouble>(TemperaturesLgm);
	xDelete<IssmDouble>(PrecipitationsPresentday);
	xDelete<IssmDouble>(tmp);
	xDelete<int>(vertexlids);

} /*}}}*/
void       Element::Delta18opdParameterization(void){/*{{{*/
	/*Are we on the base? If not, return*/
	if(!IsOnBase()) return;

	const int NUM_VERTICES 					= this->GetNumberOfVertices();
	const int NUM_VERTICES_MONTHS_PER_YEAR	= NUM_VERTICES * 12;

	int        	i,offset;
	int*        vertexlids=xNew<int>(NUM_VERTICES);
	IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* TemperaturesReconstructed=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* PrecipitationsReconstructed=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble Delta18oTime;
	IssmDouble f;
	IssmDouble time,yts,time_yr,month,time_climt,time_climp,del_clim;
	this->parameters->FindParam(&time,TimeEnum);
	this->parameters->FindParam(&yts,ConstantsYtsEnum);
	this->parameters->FindParam(&f,SmbFEnum);
	this->GetVerticesLidList(vertexlids);
	time_yr=floor(time/yts)*yts;
	time_climt=ceil(time/yts + 1e-10)*yts;
	time_climp=ceil(time/yts + 1e-10)*yts;

	/*Get some pdd parameters*/
	bool isTemperatureScaled,isPrecipScaled;
	IssmDouble dpermil=this->FindParam(SmbDpermilEnum);
	this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum);
	this->parameters->FindParam(&isPrecipScaled,SmbIsprecipscaledEnum);

	/*Recover present day temperature and precipitation*/
	DatasetInput *dinput3 = NULL;
	DatasetInput *dinput4 = NULL;
	int            offset_t,offset_p,N;
	if(!isTemperatureScaled){
		IssmDouble* time_temp_scaled = NULL;
		parameters->FindParam(&time_temp_scaled,&N,SmbTemperaturesReconstructedYearsEnum);
		if(!binary_search(&offset_t,time_climt,time_temp_scaled,N)) _error_("time not sorted?");
		if(offset_t<0) offset_t=0;
		xDelete<IssmDouble>(time_temp_scaled);
		dinput3=this->GetDatasetInput(SmbTemperaturesReconstructedEnum); _assert_(dinput3);
	}
	if(!isPrecipScaled){
		IssmDouble* time_precip_scaled = NULL;
		parameters->FindParam(&time_precip_scaled,&N,SmbPrecipitationsReconstructedYearsEnum);
		if(!binary_search(&offset_p,time_climt,time_precip_scaled,N)) _error_("time not sorted?");
		if(offset_p<0) offset_p=0;
		xDelete<IssmDouble>(time_precip_scaled);
		dinput4=this->GetDatasetInput(SmbPrecipitationsReconstructedEnum); _assert_(dinput4);
	}

	/*Get present day temp and precip (monthly)*/
	DatasetInput *dinput1 = this->GetDatasetInput(SmbTemperaturesPresentdayEnum);   _assert_(dinput1);
	DatasetInput *dinput2 = this->GetDatasetInput(SmbPrecipitationsPresentdayEnum); _assert_(dinput2);

	/*loop over vertices: */
	Gauss* gauss=this->NewGauss();
	for(int month=0;month<12;month++) {
		for(int iv=0;iv<NUM_VERTICES;iv++) {
			gauss->GaussVertex(iv);
			dinput1->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month);
			dinput2->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month);
			PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts;

			if(!isTemperatureScaled){
				dinput3->GetInputValue(&TemperaturesReconstructed[iv*12+month],gauss,offset_t*12+month);
			}
			if(!isPrecipScaled){
				dinput4->GetInputValue(&PrecipitationsReconstructed[iv*12+month],gauss,offset_p*12+month);
				PrecipitationsReconstructed[iv*12+month]=PrecipitationsReconstructed[iv*12+month]*yts;
			}
		}
	}

	/*Recover interpolation parameters at time t*/
	this->parameters->FindParam(&Delta18oTime,SmbDelta18oEnum,time);

	/*Compute the temperature and precipitation*/
	for(int iv=0;iv<NUM_VERTICES;iv++){
		ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,isTemperatureScaled,isPrecipScaled,
					f,&PrecipitationsPresentday[iv*12], &TemperaturesPresentday[iv*12],
					&PrecipitationsReconstructed[iv*12], &TemperaturesReconstructed[iv*12],
					&monthlytemperatures[iv*12], &monthlyprec[iv*12]);
	}

	/*Update inputs*/
	for (int imonth=0;imonth<12;imonth++) {
		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
		switch(this->ObjectEnum()){
			case TriaEnum:  this->inputs->SetTriaDatasetInput(SmbMonthlytemperaturesEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			case PentaEnum: this->inputs->SetPentaDatasetInput(SmbMonthlytemperaturesEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			default: _error_("Not implemented yet");
		}
		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
		switch(this->ObjectEnum()){
			case TriaEnum:  this->inputs->SetTriaDatasetInput(SmbPrecipitationEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			case PentaEnum: this->inputs->SetPentaDatasetInput(SmbPrecipitationEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			default: _error_("Not implemented yet");
		}
	}

	switch(this->ObjectEnum()){
		case TriaEnum: break;
		case PentaEnum:
		case TetraEnum:
         this->DatasetInputExtrude(SmbMonthlytemperaturesEnum,-1);
         this->DatasetInputExtrude(SmbPrecipitationEnum,-1);
         break;
		default: _error_("Not implemented yet");
	}

	/*clean-up*/
	delete gauss;
	xDelete<IssmDouble>(monthlytemperatures);
	xDelete<IssmDouble>(monthlyprec);
	xDelete<IssmDouble>(TemperaturesPresentday);
	xDelete<IssmDouble>(PrecipitationsPresentday);
	xDelete<IssmDouble>(TemperaturesReconstructed);
	xDelete<IssmDouble>(PrecipitationsReconstructed);
	xDelete<IssmDouble>(tmp);
	xDelete<int>(vertexlids);

} /*}}}*/
void       Element::SmbGradCompParameterization(void){/*{{{*/

	/*Are we on the base? If not, return*/
	if(!IsOnBase()) return;
	const int NUM_VERTICES = this->GetNumberOfVertices();

	int        i;
	IssmDouble accuref, runoffref; //reference values at given altitude
	IssmDouble accualti, runoffalti; //reference altitudes
	IssmDouble accugrad, runoffgrad; //gradients from reference altitude
	IssmDouble rho_water, rho_ice;
	IssmDouble time;

	IssmDouble*		smb	 = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble*		surf	 = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble*		accu	 = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble*		runoff = xNew<IssmDouble>(NUM_VERTICES);

	/*Get material parameters :*/
	rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
	rho_ice=this->FindParam(MaterialsRhoIceEnum);

	/*Recover parameters*/
	parameters->FindParam(&time,TimeEnum);
	parameters->FindParam(&accualti,SmbAccualtiEnum);
	parameters->FindParam(&runoffalti,SmbRunoffaltiEnum);

	/*Recover reference values at current time*/
	parameters->FindParam(&accugrad,SmbAccugradEnum,time);
	parameters->FindParam(&runoffgrad,SmbRunoffgradEnum,time);
	parameters->FindParam(&accuref,SmbAccurefEnum,time);
	parameters->FindParam(&runoffref,SmbRunoffrefEnum,time);

	/*Recover surface elevation*/
	GetInputListOnVertices(&surf[0],SurfaceEnum);

	/*Compute the temperature and precipitation*/
	for(int iv=0;iv<NUM_VERTICES;iv++){
		accu[iv]=max(0.,(accuref+(surf[iv]-accualti)*accugrad));
		runoff[iv]=max(0.,(runoffref+(surf[iv]-runoffalti)*runoffgrad));
		smb[iv]=(accu[iv]-runoff[iv])*rho_ice/rho_water;
	}

	switch(this->ObjectEnum()){
	case TriaEnum:
		this->AddInput(SmbMassBalanceSubstepEnum,&smb[0],P1Enum);
		this->AddInput(SmbRunoffSubstepEnum,&runoff[0],P1Enum);
		break;
	case PentaEnum:
		this->AddInput(SmbMassBalanceSubstepEnum,&smb[0],P1Enum);
		this->AddInput(SmbRunoffSubstepEnum,&runoff[0],P1Enum);
		this->InputExtrude(SmbMassBalanceSubstepEnum,-1);
		this->InputExtrude(SmbRunoffSubstepEnum,-1);
		break;
	default: _error_("Not implemented yet");
	}
	/*clean-up*/
	xDelete<IssmDouble>(surf);
	xDelete<IssmDouble>(accu);
	xDelete<IssmDouble>(runoff);
	xDelete<IssmDouble>(smb);
}
/*}}}*/
IssmDouble Element::Divergence(void){/*{{{*/
	/*Compute element divergence*/

	/*Intermediaries*/
	int        dim;
	IssmDouble Jdet;
	IssmDouble divergence=0.;
	IssmDouble dvx[3],dvy[3],dvz[3];
	IssmDouble *xyz_list = NULL;

	/*Get inputs and parameters*/
	this->FindParam(&dim,DomainDimensionEnum);
	Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input);
	Input* vz_input = NULL;
	if(dim==3){
		vz_input = this->GetInput(VzEnum); _assert_(vz_input);
	}
	this->GetVerticesCoordinates(&xyz_list);

	Gauss* gauss=this->NewGauss(5);
	while(gauss->next()){

		this->JacobianDeterminant(&Jdet,xyz_list,gauss);

		/*Get strain rate assuming that epsilon has been allocated*/
		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
		if(dim==2){
			divergence += (dvx[0]+dvy[1])*gauss->weight*Jdet;
		}
		else{
			vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
			divergence += (dvx[0]+dvy[1]+dvz[2])*gauss->weight*Jdet;
		}

	}

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	delete gauss;
	return divergence;
}/*}}}*/
void       Element::dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/

	/*Intermediaries*/
	int materialstype;
	IssmDouble dmudB;
	IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,exy,exz,eyz];    */
	IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
	IssmDouble eps_eff;
	IssmDouble eps0=1.e-27;

	if(dim==3){
		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
		this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
	}
	else{
		/* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
	}
	/*Get viscosity*/
	materialstype=this->material->ObjectEnum();
	switch(materialstype){
		case MaticeEnum:
			material->GetViscosity_B(&dmudB,eps_eff,gauss);
			break;
		case MatestarEnum:
			material->ViscosityBFS(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,vz_input,eps_eff);
			break;
		default: _error_("not supported");
	}

	/*Assign output pointer*/
	*pdmudB=dmudB;

}
/*}}}*/
void       Element::dViscositydBHO(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	int materialstype;
	IssmDouble dmudB;
	IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exy,exz,eyz];    */
	IssmDouble epsilon2d[2];/* epsilon=[exx,eyy,exy];    */
	IssmDouble eps_eff;
	IssmDouble eps0=1.e-27;

	if(dim==3){
		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
		this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
	}
	else{
		/* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
		this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1] + eps0*eps0);
	}
	/*Get viscosity*/
	materialstype=this->material->ObjectEnum();
	switch(materialstype){
		case MaticeEnum:
			material->GetViscosity_B(&dmudB,eps_eff,gauss);
			break;
		case MatestarEnum:
			material->ViscosityBHO(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,eps_eff);
			break;
		default: _error_("not supported");
	}

	/*Assign output pointer*/
	*pdmudB=dmudB;

}
/*}}}*/
void       Element::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	int materialstype;
	IssmDouble dmudB;
	IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
	IssmDouble epsilon1d;   /* epsilon=[exx];    */
	IssmDouble eps_eff;

	if(dim==2){
		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
		eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
	}
	else{
		/* eps_eff^2 = 1/2 exx^2*/
		this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
		eps_eff = sqrt(epsilon1d*epsilon1d/2.);
	}

	/*Get viscosity*/
	materialstype=this->material->ObjectEnum();
	switch(materialstype){
		case MaticeEnum:
			material->GetViscosity_B(&dmudB,eps_eff,gauss);
			break;
		case MatestarEnum:
			material->ViscosityBSSA(&dmudB,dim,xyz_list,gauss,vx_input,vy_input,eps_eff);
			break;
		default: _error_("not supported");
	}

	/*Assign output pointer*/
	*pdmudB=dmudB;

}
/*}}}*/
void       Element::dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble dmudB;
	IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
	IssmDouble epsilon1d;   /* epsilon=[exx];    */
	IssmDouble eps_eff;

	if(dim==2){
		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
		eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
	}
	else{
		/* eps_eff^2 = 1/2 exx^2*/
		this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
		eps_eff = sqrt(epsilon1d*epsilon1d/2.);
	}

	/*Get viscosity*/
	material->GetViscosity_D(&dmudB,eps_eff,gauss);

	/*Assign output pointer*/
	*pdmudB=dmudB;

}
/*}}}*/
void       Element::Echo(void){/*{{{*/
	_printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
	_printf_("   id : "<<this->id <<"\n");
	_printf_("   sid: "<<this->sid<<"\n");
	_printf_("   lid: "<<this->lid<<"\n");
	if(vertices){
		const int NUM_VERTICES = this->GetNumberOfVertices();
		for(int i=0;i<NUM_VERTICES;i++) vertices[i]->Echo();
	}
	else _printf_("vertices = NULL\n");

	if(nodes){
		int numnodes = this->GetNumberOfNodes();
		for(int i=0;i<numnodes;i++) {
			_printf_("nodes[" << i << "] = " << nodes[i]);
			nodes[i]->Echo();
		}
	}
	else _printf_("nodes = NULL\n");

	if (material) material->Echo();
	else _printf_("material = NULL\n");

	_printf_("   parameters\n");
	if (parameters) parameters->Echo();
	else _printf_("parameters = NULL\n");

	_printf_("   inputs\n");
	if (inputs) inputs->Echo();
	else _printf_("inputs=NULL\n");
}
/*}}}*/
void       Element::FindParam(bool* pvalue,int paramenum){/*{{{*/
	this->parameters->FindParam(pvalue,paramenum);
}/*}}}*/
void       Element::FindParam(int* pvalue,int paramenum){/*{{{*/
	this->parameters->FindParam(pvalue,paramenum);
}/*}}}*/
void       Element::FindParam(IssmDouble* pvalue,int paramenum){/*{{{*/
	this->parameters->FindParam(pvalue,paramenum);
}/*}}}*/
IssmDouble Element::FindParam(int paramenum){/*{{{*/
	return this->parameters->FindParam(paramenum);
}/*}}}*/
void       Element::FindParam(int** pvalues,int* psize,int paramenum){/*{{{*/
	this->parameters->FindParam(pvalues,psize,paramenum);
}/*}}}*/
IssmDouble Element::FloatingArea(IssmDouble* mask, bool scaled){/*{{{*/

	/*Retrieve values of the mask defining the element: */
	for(int i=0;i<this->GetNumberOfVertices();i++){
		if(mask[this->vertices[i]->Sid()]<=0.){
			return 0.;
		}
	}

	/*Return: */
	return this->FloatingArea(scaled);
}
/*}}}*/
void       Element::GetDofList(int** pdoflist,int approximation_enum,int setenum){/*{{{*/

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = this->GetNumberOfNodes();

	/*First, figure out size of doflist and create it: */
	int numberofdofs=0;
	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);

	/*Allocate output*/
	int* doflist=xNew<int>(numberofdofs);

	/*Populate: */
	int count=0;
	for(int i=0;i<numnodes;i++){
		nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
		count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
	}

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
void       Element::GetDofListLocal(int** pdoflist,int approximation_enum,int setenum){/*{{{*/

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = this->GetNumberOfNodes();

	/*First, figure out size of doflist and create it: */
	int numberofdofs=0;
	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);

	/*Allocate output*/
	int* doflist=xNew<int>(numberofdofs);

	/*Populate: */
	int count=0;
	for(int i=0;i<numnodes;i++){
		nodes[i]->GetDofListLocal(doflist+count,approximation_enum,setenum);
		count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
	}

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
void       Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/

	/*Fetch number of nodes and dof for this finite element*/
	int vnumnodes = this->NumberofNodesVelocity();
	int pnumnodes = this->NumberofNodesPressure();

	/*First, figure out size of doflist and create it: */
	int numberofdofs=0;
	for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);

	/*Allocate output*/
	int* doflist=xNew<int>(numberofdofs);

	/*Populate: */
	int count=0;
	for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++){
		nodes[i]->GetDofList(doflist+count,FSApproximationEnum,setenum);
		count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
	}

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
void       Element::GetDofListVelocity(int** pdoflist,int setenum){/*{{{*/

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = this->NumberofNodesVelocity();

	/*First, figure out size of doflist and create it: */
	int numberofdofs=0;
	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);

	/*Allocate output*/
	int* doflist=xNew<int>(numberofdofs);

	/*Populate: */
	int count=0;
	for(int i=0;i<numnodes;i++){
		nodes[i]->GetDofList(doflist+count,FSvelocityEnum,setenum);
		count+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
	}

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
void       Element::GetDofListLocalPressure(int** pdoflist,int setenum){/*{{{*/

	/*Fetch number of nodes and dof for this finite element*/
	int vnumnodes = this->NumberofNodesVelocity();
	int pnumnodes = this->NumberofNodesPressure();

	/*First, figure out size of doflist and create it: */
	int numberofdofs=0;
	for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);

	/*Allocate output*/
	int* doflist=xNew<int>(numberofdofs);

	/*Populate: */
	int count=0;
	for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++){
		nodes[i]->GetDofListLocal(doflist+count,FSApproximationEnum,setenum);
		count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
	}

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
void       Element::GetDofListLocalVelocity(int** pdoflist,int setenum){/*{{{*/

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = this->NumberofNodesVelocity();

	/*First, figure out size of doflist and create it: */
	int numberofdofs=0;
	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);

	/*Allocate output*/
	int* doflist=xNew<int>(numberofdofs);

	/*Populate: */
	int count=0;
	for(int i=0;i<numnodes;i++){
		nodes[i]->GetDofListLocal(doflist+count,FSvelocityEnum,setenum);
		count+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
	}

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
void       Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
	Input *input    = this->GetInput(enumtype);
	this->GetInputListOnNodes(pvalue,input,defaultvalue);
}
/*}}}*/
void       Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype){/*{{{*/

	Input *input    = this->GetInput(enumtype);
	if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
	this->GetInputListOnNodes(pvalue,input,0.);

}
/*}}}*/
void       Element::GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype){/*{{{*/

	_assert_(pvalue);

	int     numnodes = this->NumberofNodesVelocity();
	Input *input    = this->GetInput(enumtype);
	if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");

	/* Start looping on the number of vertices: */
	Gauss* gauss=this->NewGauss();
	for(int iv=0;iv<numnodes;iv++){
		gauss->GaussNode(this->VelocityInterpolation(),iv);
		input->GetInputValue(&pvalue[iv],gauss);
	}
	delete gauss;
}
/*}}}*/
void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){/*{{{*/

	/*Recover input*/
	Input* input2=this->GetInput(enumtype);
	if(!input2) _error_("input "<<EnumToStringx(enumtype)<<" not found in element");
	this->GetInputListOnVertices(pvalue,input2,0.);
}
/*}}}*/
void       Element::GetInputListOnVerticesAtTime(IssmDouble* pvalue, int enumtype, IssmDouble time){/*{{{*/

	/*Recover input*/
	Input* input=this->GetInput(enumtype,time);
	if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
	this->GetInputListOnVertices(pvalue,input,0.);
}
/*}}}*/
void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
	Input* input=this->GetInput(enumtype);
	this->GetInputListOnVertices(pvalue,input,defaultvalue);
}
/*}}}*/
void       Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug){/*{{{*/

	/*Get number of nodes for this element*/
	int numnodes = this->GetNumberOfNodes();

	/*Some checks to avoid segmentation faults*/
	_assert_(ug);
	_assert_(numnodes>0);
	_assert_(nodes);

	/*Get element minimum/maximum*/
	IssmDouble input_min = ug[nodes[0]->GetDof(0,GsetEnum)];
	IssmDouble input_max = input_min;
	for(int i=1;i<numnodes;i++){
		if(ug[nodes[i]->GetDof(0,GsetEnum)] < input_min) input_min = ug[nodes[i]->GetDof(0,GsetEnum)];
		if(ug[nodes[i]->GetDof(0,GsetEnum)] > input_max) input_max = ug[nodes[i]->GetDof(0,GsetEnum)];
	}

	/*Second loop to reassign min and max with local extrema*/
	for(int i=0;i<numnodes;i++){
		if(min[nodes[i]->GetDof(0,GsetEnum)]>input_min) min[nodes[i]->GetDof(0,GsetEnum)] = input_min;
		if(max[nodes[i]->GetDof(0,GsetEnum)]<input_max) max[nodes[i]->GetDof(0,GsetEnum)] = input_max;
	}
}
/*}}}*/
void       Element::GetInputValue(bool* pvalue,int inputenum){/*{{{*/

	this->inputs->GetInputValue(pvalue,inputenum,this->lid);

}/*}}}*/
void       Element::GetInputValue(int* pvalue,int inputenum){/*{{{*/

	this->inputs->GetInputValue(pvalue,inputenum,this->lid);

}/*}}}*/
void       Element::GetInputValue(IssmDouble* pvalue,int inputenum){/*{{{*/

	this->inputs->GetInputValue(pvalue,inputenum,this->lid);

}/*}}}*/
void       Element::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum){/*{{{*/

	Input* input=this->GetInput(inputenum);
	if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
	input->GetInputValue(pvalue,gauss);

}/*}}}*/
Node*      Element::GetNode(int nodeindex){/*{{{*/
	_assert_(nodeindex>=0);
	_assert_(nodeindex<this->GetNumberOfNodes(this->element_type));
	return this->nodes[nodeindex];
}/*}}}*/
int        Element::GetNodeIndex(Node* node){/*{{{*/

	_assert_(this->nodes);
	int numnodes = this->GetNumberOfNodes(this->element_type);

	for(int i=0;i<numnodes;i++){
		if(node==nodes[i]) return i;
	}
	_error_("Node provided not found among element nodes");

}
/*}}}*/
void       Element::GetNodesLidList(int* lidlist){/*{{{*/

	_assert_(lidlist);
	_assert_(nodes);
	int numnodes = this->GetNumberOfNodes();
	for(int i=0;i<numnodes;i++){
		lidlist[i]=nodes[i]->Lid();
	}
}
/*}}}*/
void       Element::GetNodesSidList(int* sidlist){/*{{{*/

	_assert_(sidlist);
	_assert_(nodes);
	int numnodes = this->GetNumberOfNodes();
	for(int i=0;i<numnodes;i++){
		sidlist[i]=nodes[i]->Sid();
	}
}
/*}}}*/
void       Element::GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity){/*{{{*/
	/*Compute deformational heating from epsilon and viscosity */

	IssmDouble epsilon_matrix[3][3];
	IssmDouble epsilon_eff;
	IssmDouble epsilon_sqr[3][3];

	/* Build epsilon matrix */
	epsilon_matrix[0][0]=epsilon[0];
	epsilon_matrix[1][0]=epsilon[3];
	epsilon_matrix[2][0]=epsilon[4];
	epsilon_matrix[0][1]=epsilon[3];
	epsilon_matrix[1][1]=epsilon[1];
	epsilon_matrix[2][1]=epsilon[5];
	epsilon_matrix[0][2]=epsilon[4];
	epsilon_matrix[1][2]=epsilon[5];
	epsilon_matrix[2][2]=epsilon[2];

	/* Effective value of epsilon_matrix */
	epsilon_sqr[0][0]=epsilon_matrix[0][0]*epsilon_matrix[0][0];
	epsilon_sqr[1][0]=epsilon_matrix[1][0]*epsilon_matrix[1][0];
	epsilon_sqr[2][0]=epsilon_matrix[2][0]*epsilon_matrix[2][0];
	epsilon_sqr[0][1]=epsilon_matrix[0][1]*epsilon_matrix[0][1];
	epsilon_sqr[1][1]=epsilon_matrix[1][1]*epsilon_matrix[1][1];
	epsilon_sqr[2][1]=epsilon_matrix[2][1]*epsilon_matrix[2][1];
	epsilon_sqr[0][2]=epsilon_matrix[0][2]*epsilon_matrix[0][2];
	epsilon_sqr[1][2]=epsilon_matrix[1][2]*epsilon_matrix[1][2];
	epsilon_sqr[2][2]=epsilon_matrix[2][2]*epsilon_matrix[2][2];
	epsilon_eff=1/sqrt(2.)*sqrt(epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]);

	/*Phi = Tr(sigma * eps)
	 *    = Tr(sigma'* eps)
	 *    = 2 * eps_eff * sigma'_eff
	 *    = 4 * mu * eps_eff ^2*/
	*phi=4.*epsilon_eff*epsilon_eff*viscosity;
}
/*}}}*/
void       Element::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/

	int        *doflist = NULL;
	IssmDouble  value;

	/*Fetch number of nodes for this finite element*/
	int numnodes = this->GetNumberOfNodes();

	/*Fetch dof list and allocate solution vector*/
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
	IssmDouble* values = xNew<IssmDouble>(numnodes);

	/*Get inputs*/
	Input* enum_input=this->GetInput(enum_type); _assert_(enum_input);

	/*Ok, we have the values, fill in the array: */
	Gauss* gauss=this->NewGauss();
	for(int i=0;i<numnodes;i++){
		gauss->GaussNode(this->element_type,i);

		enum_input->GetInputValue(&value,gauss);
		values[i]=value;
	}

	solution->SetValues(numnodes,doflist,values,INS_VAL);

	/*Free ressources:*/
	xDelete<int>(doflist);
	xDelete<IssmDouble>(values);
	delete gauss;
}
/*}}}*/
/* void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){/\*{{{*\/ */

/* 	/\*Fetch number vertices for this element and allocate arrays*\/ */
/* 	const int NUM_VERTICES = this->GetNumberOfVertices(); */
/* 	int*        vertexpidlist = xNew<int>(NUM_VERTICES); */
/* 	IssmDouble* values        = xNew<IssmDouble>(NUM_VERTICES); */

/* 	/\*Fill in values*\/ */
/* 	this->GetVerticesPidList(vertexpidlist); */
/* 	this->GetInputListOnVertices(values,input_enum); */
/* 	vector->SetValues(NUM_VERTICES,vertexpidlist,values,INS_VAL); */

/* 	/\*Clean up*\/ */
/* 	xDelete<int>(vertexpidlist); */
/* 	xDelete<IssmDouble>(values); */

/* } */
/* /\*}}}*\/ */
void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum,int type){/*{{{*/

	/*Fetch number vertices for this element and allocate arrays*/
	const int NUM_VERTICES = this->GetNumberOfVertices();

	int         numnodes    = this->GetNumberOfNodes();
	int*        doflist     = NULL;
	IssmDouble  value;
	IssmDouble* values      = NULL;
	Input*     input       = NULL;

	switch(type){
		case ElementSIdEnum:
			input=this->GetInput(input_enum); _assert_(input);
			input->GetInputAverage(&value);
			vector->SetValue(this->sid,value,INS_VAL);
			break;
		case VertexPIdEnum:
			doflist = xNew<int>(NUM_VERTICES);
			values = xNew<IssmDouble>(NUM_VERTICES);
			/*Fill in values*/
			this->GetVerticesPidList(doflist);
			this->GetInputListOnVertices(values,input_enum);
			vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
			break;
		case VertexSIdEnum:
			doflist = xNew<int>(NUM_VERTICES);
			values = xNew<IssmDouble>(NUM_VERTICES);
			/*Fill in values*/
			this->GetVerticesSidList(doflist);
			this->GetInputListOnVertices(values,input_enum);
			vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
			break;
		case NodesEnum:
			doflist = xNew<int>(numnodes);
			values = xNew<IssmDouble>(numnodes);
			/*Fill in values*/
			this->GetInputListOnNodes(values,input_enum);
			this->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
			vector->SetValues(numnodes,doflist,values,INS_VAL);
			break;
		case NodeSIdEnum:
			doflist = xNew<int>(numnodes);
			values = xNew<IssmDouble>(numnodes);
			/*Fill in values*/
			this->GetNodesSidList(doflist);
			this->GetInputListOnNodes(values,input_enum);
			vector->SetValues(numnodes,doflist,values,INS_VAL);
			break;
		default:
			_error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
	}

	/*Clean up*/
	xDelete<int>(doflist);
	xDelete<IssmDouble>(values);

}
/*}}}*/
void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum,int type, IssmDouble time){/*{{{*/

	/*Fetch number vertices for this element and allocate arrays*/
	const int NUM_VERTICES = this->GetNumberOfVertices();

	int         numnodes    = this->GetNumberOfNodes();
	int*        doflist     = NULL;
	IssmDouble* values      = NULL;

	switch(type){
		case VertexPIdEnum:
			doflist = xNew<int>(NUM_VERTICES);
			values = xNew<IssmDouble>(NUM_VERTICES);
			/*Fill in values*/
			this->GetVerticesPidList(doflist);
			this->GetInputListOnVerticesAtTime(values,input_enum,time);
			vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
			break;
		case VertexSIdEnum:
			doflist = xNew<int>(NUM_VERTICES);
			values = xNew<IssmDouble>(NUM_VERTICES);
			/*Fill in values*/
			this->GetVerticesSidList(doflist);
			this->GetInputListOnVerticesAtTime(values,input_enum,time);
			vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
			break;
		default:
			_error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
	}

	/*Clean up*/
	xDelete<int>(doflist);
	xDelete<IssmDouble>(values);

}
/*}}}*/
void       Element::GetVerticesLidList(int* lidlist){/*{{{*/

	const int NUM_VERTICES = this->GetNumberOfVertices();
	for(int i=0;i<NUM_VERTICES;i++) lidlist[i]=vertices[i]->Lid();

}
/*}}}*/
void       Element::GetVerticesPidList(int* pidlist){/*{{{*/

	const int NUM_VERTICES = this->GetNumberOfVertices();
	for(int i=0;i<NUM_VERTICES;i++) pidlist[i]=vertices[i]->Pid();

}
/*}}}*/
void       Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/

	const int NUM_VERTICES = this->GetNumberOfVertices();
	for(int i=0;i<NUM_VERTICES;i++) connectivity[i]=this->vertices[i]->Connectivity();
}
/*}}}*/
void       Element::GetVerticesCoordinates(IssmDouble** pxyz_list){/*{{{*/

	const int NUM_VERTICES = this->GetNumberOfVertices();

	IssmDouble* xyz_list    = xNew<IssmDouble>(NUM_VERTICES*3);
	::GetVerticesCoordinates(xyz_list,this->vertices,NUM_VERTICES);

	*pxyz_list = xyz_list;

}/*}}}*/
void       Element::GetVerticesSidList(int* sidlist){/*{{{*/

	const int NUM_VERTICES = this->GetNumberOfVertices();
	for(int i=0;i<NUM_VERTICES;i++) sidlist[i]=this->vertices[i]->Sid();
}
/*}}}*/
IssmDouble Element::GetXcoord(IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

	/*output*/
	IssmDouble x;

	/*Create list of x*/
	const int NUM_VERTICES = this->GetNumberOfVertices();

	IssmDouble* x_list = xNew<IssmDouble>(NUM_VERTICES);

	for(int i=0;i<NUM_VERTICES;i++) x_list[i]=xyz_list[i*3+0];
	ValueP1OnGauss(&x,x_list,gauss);

	xDelete<IssmDouble>(x_list);
	return x;
}/*}}}*/
IssmDouble Element::GetYcoord(IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

	/*output*/
	IssmDouble y;

	/*Create list of y*/
	const int NUM_VERTICES = this->GetNumberOfVertices();

	IssmDouble* y_list      = xNew<IssmDouble>(NUM_VERTICES);

	for(int i=0;i<NUM_VERTICES;i++) y_list[i]=xyz_list[i*3+1];
	ValueP1OnGauss(&y,y_list,gauss);

	xDelete<IssmDouble>(y_list);
	return y;
}/*}}}*/
IssmDouble Element::GetZcoord(IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

	/*output*/
	IssmDouble z;

	/*Create list of z*/
	const int NUM_VERTICES = this->GetNumberOfVertices();

	IssmDouble* z_list      = xNew<IssmDouble>(NUM_VERTICES);

	for(int i=0;i<NUM_VERTICES;i++) z_list[i]=xyz_list[i*3+2];
	ValueP1OnGauss(&z,z_list,gauss);

	xDelete<IssmDouble>(z_list);
	return z;
}/*}}}*/
void       Element::GradientIndexing(int* indexing,int control_index){/*{{{*/

	/*Get number of controls*/
	int num_controls;
	bool isautodiff;
	int* N=NULL;
	int* M=NULL;
	int* interp=NULL;
	parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
	parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
	parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
	parameters->FindParam(&interp,NULL,ControlInputInterpolationEnum);

	/*Get number of vertices*/
	const int NUM_VERTICES = this->GetNumberOfVertices();

	/*Get starting index of gradient for this control*/
	int start = 0;
	for(int n=0;n<control_index;n++) start+=N[n]*M[n];

	/*Create index*/
	if(interp[control_index]==P1Enum){
		for(int n=0;n<N[control_index];n++){
			for(int i=0;i<NUM_VERTICES;i++){
				indexing[i+n*NUM_VERTICES]= start + n*M[control_index] + this->vertices[i]->Sid();
			}
		}
	}
	else if(interp[control_index]==P0Enum){
		for(int n=0;n<N[control_index];n++){
				indexing[n]= start + n*M[control_index] + this->Sid();
		}
	}
	else{
		_error_("interpolation not supported");
	}

	xDelete<int>(M);
	xDelete<int>(N);
	xDelete<int>(interp);
}
/*}}}*/
IssmDouble Element::GroundedArea(IssmDouble* mask, bool scaled){/*{{{*/

	/*Retrieve values of the mask defining the element: */
	for(int i=0;i<this->GetNumberOfVertices();i++){
		if(mask[this->vertices[i]->Sid()]<=0.){
			return 0.;
		}
	}

	/*Return: */
	return this->GroundedArea(scaled);
}
/*}}}*/
bool       Element::HasNodeOnBase(){/*{{{*/
	Input* input=this->GetInput(MeshVertexonbaseEnum); _assert_(input);
	return (input->GetInputMax()>0.);
}/*}}}*/
bool       Element::HasNodeOnSurface(){/*{{{*/
	Input* input=this->GetInput(MeshVertexonsurfaceEnum); _assert_(input);
	return (input->GetInputMax()>0.);
}/*}}}*/
IssmDouble Element::IceMass(bool scaled){/*{{{*/

	IssmDouble rho_ice;

	if(!IsIceInElement())return 0.; //do not contribute to the volume of the ice!

	/*recover ice density: */
	rho_ice=FindParam(MaterialsRhoIceEnum);

	return rho_ice*this->IceVolume(scaled);
}
/*}}}*/
IssmDouble Element::IceMass(IssmDouble* mask, bool scaled){/*{{{*/

	/*Retrieve values of the mask defining the element: */
	for(int i=0;i<this->GetNumberOfVertices();i++){
		if(mask[this->vertices[i]->Sid()]<=0.){
			return 0.;
		}
	}

	/*Return: */
	return this->IceMass(scaled);
}
/*}}}*/
IssmDouble Element::IceVolume(IssmDouble* mask, bool scaled){/*{{{*/

	/*Retrieve values of the mask defining the element: */
	for(int i=0;i<this->GetNumberOfVertices();i++){
		if(mask[this->vertices[i]->Sid()]<=0.){
			return 0.;
		}
	}

	/*Return: */
	return this->IceVolume(scaled);
}
/*}}}*/
IssmDouble Element::IceVolumeAboveFloatation(IssmDouble* mask, bool scaled){/*{{{*/

	/*Retrieve values of the mask defining the element: */
	for(int i=0;i<this->GetNumberOfVertices();i++){
		if(mask[this->vertices[i]->Sid()]<=0.){
			return 0.;
		}
	}

	/*Return: */
	return this->IceVolumeAboveFloatation(scaled);
}
/*}}}*/
int        Element::Id(){/*{{{*/

	return this->id;

}
/*}}}*/
void       Element::InputCreate(IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){/*{{{*/

	/*Intermediaries*/
	int i,t;

	/*Branch on type of vector: nodal or elementary: */
	if(vector_type==1){ //nodal vector

		const int NUM_VERTICES = this->GetNumberOfVertices();

		int        *vertexids  = xNew<int>(NUM_VERTICES);
		int        *vertexlids = xNew<int>(NUM_VERTICES);
		IssmDouble *values     = xNew<IssmDouble>(NUM_VERTICES);

		/*Recover vertices ids needed to initialize inputs*/
		_assert_(iomodel->elements);
		for(i=0;i<NUM_VERTICES;i++){
			vertexids[i] =reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
			vertexlids[i]=iomodel->my_vertices_lids[vertexids[i]-1];
		}

		/*Are we in transient or static? */
		if(M==1){
			values[0]=vector[0];
			this->SetElementInput(inputs,vector_enum,vector[0]);
		}
		else if(M==iomodel->numberofvertices){
			for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
			this->SetElementInput(inputs,NUM_VERTICES,vertexlids,values,vector_enum);
		}
		else if(M==iomodel->numberofvertices+1){
			/*create transient input: */
			IssmDouble* times = xNew<IssmDouble>(N);
			for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
			inputs->SetTransientInput(vector_enum,times,N);
			TransientInput* transientinput = inputs->GetTransientInput(vector_enum);
			for(t=0;t<N;t++){
				for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
				switch(this->ObjectEnum()){
					case TriaEnum:  transientinput->AddTriaTimeInput( t,NUM_VERTICES,vertexlids,values,P1Enum); break;
					case PentaEnum: transientinput->AddPentaTimeInput(t,NUM_VERTICES,vertexlids,values,P1Enum); break;
					default: _error_("Not implemented yet");
				}
			}
			xDelete<IssmDouble>(times);
		}
		else if(M==iomodel->numberofelements){

			/*This is a Patch!*/
			xDelete<IssmDouble>(values);
			values = xNew<IssmDouble>(N);
			for(int j=0;j<N;j++) values[j]=vector[this->Sid()*N+j];

			if (N==this->GetNumberOfNodes(P1Enum)){
				this->SetElementInput(inputs,NUM_VERTICES,vertexlids,values,vector_enum);
			}
			else if(N==this->GetNumberOfNodes(P0Enum)){
				this->SetElementInput(inputs,vector_enum,values[0]);
			}
			else if(N==this->GetNumberOfNodes(P1xP2Enum)){ _assert_(this->ObjectEnum()==PentaEnum);
				inputs->SetPentaInput(vector_enum,P1xP2Enum,this->lid,N,values);
			}
			else if(N==this->GetNumberOfNodes(P1xP3Enum)){ _assert_(this->ObjectEnum()==PentaEnum);
				inputs->SetPentaInput(vector_enum,P1xP3Enum,this->lid,N,values);
			}
			else{
				_error_("Patch interpolation not supported yet");
			}

		}
		else{
			_error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
		}

		xDelete<IssmDouble>(values);
		xDelete<int>(vertexids);
		xDelete<int>(vertexlids);
	}
	else if(vector_type==2){ //element vector

		IssmDouble value;

		/*Are we in transient or static? */
		if(M==iomodel->numberofelements){
			if (code==5){ //boolean
				this->SetBoolInput(inputs,vector_enum,reCast<bool>(vector[this->Sid()]));
			}
			else if (code==6){ //integer
				this->SetIntInput(inputs,vector_enum,reCast<int>(vector[this->Sid()]));
			}
			else if (code==7){ //IssmDouble
				this->SetElementInput(inputs,vector_enum,vector[this->Sid()]);
			}
			else _error_("could not recognize nature of vector from code " << code);
		}
		else if(M==iomodel->numberofelements+1){
			/*create transient input: */
			IssmDouble* times = xNew<IssmDouble>(N);
			for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
			inputs->SetTransientInput(vector_enum,times,N);
			TransientInput* transientinput = inputs->GetTransientInput(vector_enum);
			for(t=0;t<N;t++){
				value=vector[N*this->Sid()+t];
				switch(this->ObjectEnum()){
					case TriaEnum:  transientinput->AddTriaTimeInput( t,1,&(this->lid),&value,P0Enum); break;
					case PentaEnum: transientinput->AddPentaTimeInput(t,1,&(this->lid),&value,P0Enum); break;
					default: _error_("Not implemented yet");
				}
			}
			xDelete<IssmDouble>(times);
		}
		else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
	}
	else if(vector_type==3){ //Double array matrix

		/*For right now we are static */
		if(M==iomodel->numberofelements){
			IssmDouble* layers = xNewZeroInit<IssmDouble>(N);
			for(t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
			inputs->SetArrayInput(vector_enum,this->lid,layers,N);
			xDelete<IssmDouble>(layers);
		}
		else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
	}
	else _error_("Cannot add input for vector type " << vector_type << " (not supported)");
}
/*}}}*/
void       Element::ControlInputCreate(IssmDouble* vector,IssmDouble* min_vector,IssmDouble* max_vector,Inputs* inputs,IoModel* iomodel,int M,int N,IssmDouble scale,int input_enum,int id){/*{{{*/

	/*Intermediaries*/
	const int numvertices = this->GetNumberOfVertices();
	IssmDouble  value,value_min,value_max;

	int         vertexids[MAXVERTICES];
	int         vertexlids[MAXVERTICES];
	IssmDouble  values[MAXVERTICES];
	IssmDouble  values_min[MAXVERTICES];
	IssmDouble  values_max[MAXVERTICES];

	/*Some sanity checks*/
	_assert_(vector);
	_assert_(min_vector);
	_assert_(max_vector);

	/*Recover vertices ids needed to initialize inputs*/
	_assert_(iomodel->elements);
	for(int i=0;i<numvertices;i++){
		vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
		vertexlids[i]=iomodel->my_vertices_lids[vertexids[i]-1];
	}

	/*Are we in transient or static? */
	if(M==iomodel->numberofvertices && N==1){
		for(int i=0;i<numvertices;i++){
			values[i]     = vector[vertexids[i]-1];
			values_min[i] = scale*min_vector[vertexids[i]-1];
			values_max[i] = scale*max_vector[vertexids[i]-1];
		}
		this->AddControlInput(input_enum,inputs,iomodel,&values[0],&values_min[0],&values_max[0],P1Enum,id);
	}
	else if(M==iomodel->numberofelements && N==1){
		values[0]     = vector[this->Sid()];
		values_min[0] = scale*min_vector[this->Sid()];
		values_max[0] = scale*max_vector[this->Sid()];
		this->AddControlInput(input_enum,inputs,iomodel,&values[0],&values_min[0],&values_max[0],P0Enum,id);
	}
	else if(M==iomodel->numberofelements+1 && N>1){

		/*create transient input: */
		IssmDouble* times = xNew<IssmDouble>(N);
		for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
		inputs->SetTransientControlInput(input_enum,id,times,N);

		/*Get input*/
		ControlInput* control_input = inputs->GetControlInput(input_enum); _assert_(control_input);
		TransientInput* transientinput_val = control_input->GetTransientInput("value");
		TransientInput* transientinput_min = control_input->GetTransientInput("lowerbound");
		TransientInput* transientinput_max = control_input->GetTransientInput("upperbound");
		for(int t=0;t<N;t++){
			value     = vector[N*this->Sid()+t];
			value_min = scale*min_vector[N*this->Sid()+t];
			value_max = scale*max_vector[N*this->Sid()+t];
			switch(this->ObjectEnum()){
				case TriaEnum:
					transientinput_val->AddTriaTimeInput( t,1,&(this->lid),&value,P0Enum);
					transientinput_min->AddTriaTimeInput( t,1,&(this->lid),&value_min,P0Enum);
					transientinput_max->AddTriaTimeInput( t,1,&(this->lid),&value_max,P0Enum);
					break;
				case PentaEnum:
					//transientinput->AddPentaTimeInput(t,1,&(this->lid),&value,P0Enum);
					_error_("to be implemented");
					break;
				default: _error_("Not implemented yet");
			}
		}
		xDelete<IssmDouble>(times);
	}
	else if(M==iomodel->numberofvertices+1 && N>1){

		/*create transient input: */
		IssmDouble* times = xNew<IssmDouble>(N);
		for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
		inputs->SetTransientControlInput(input_enum,id,times,N);

		/*Get input*/
		ControlInput* control_input = inputs->GetControlInput(input_enum); _assert_(control_input);
		TransientInput* transientinput_val = control_input->GetTransientInput("value");
		TransientInput* transientinput_min = control_input->GetTransientInput("lowerbound");
		TransientInput* transientinput_max = control_input->GetTransientInput("upperbound");
		for(int t=0;t<N;t++){
			for(int i=0;i<numvertices;i++){
				values[i]     = vector[N*(vertexids[i]-1)+t];
				values_min[i] = scale*min_vector[N*(vertexids[i]-1)+t];
				values_max[i] = scale*max_vector[N*(vertexids[i]-1)+t];
			}
			switch(this->ObjectEnum()){
				case TriaEnum:
					transientinput_val->AddTriaTimeInput( t,numvertices,vertexlids,&values[0],P1Enum);
					transientinput_min->AddTriaTimeInput( t,numvertices,vertexlids,&values_min[0],P1Enum);
					transientinput_max->AddTriaTimeInput( t,numvertices,vertexlids,&values_max[0],P1Enum);
					break;
				case PentaEnum:
					//transientinput->AddPentaTimeInput(t,1,&(this->lid),&value,P0Enum);
					_error_("to be implemented");
					break;
				default: _error_("Not implemented yet");
			}
		}
		xDelete<IssmDouble>(times);
	}
	else _error_("not currently supported type of M and N attempted");
}/*}}}*/
void       Element::DatasetInputAdd(int enum_type,IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int input_enum,int code,int input_id){/*{{{*/
	/*enum_type: the name of the DatasetInput (eg Outputdefinition1)
	 * vector: information being stored (eg observations)
	 * vector_type: is if by element or by vertex
	 * input_enum: is the name of the vector being stored
	 * code: what type of data is in the vector (booleans, ints, doubles)
	 */

	/*Intermediaries*/
	int i,t;

	/*Branch on type of vector: nodal or elementary: */
	if(vector_type==1){ //nodal vector

		const int NUM_VERTICES = this->GetNumberOfVertices();

		int        *vertexids  = xNew<int>(NUM_VERTICES);
		int        *vertexlids = xNew<int>(NUM_VERTICES);
		IssmDouble *values     = xNew<IssmDouble>(NUM_VERTICES);

		/*Recover vertices ids needed to initialize inputs*/
		_assert_(iomodel->elements);
		for(i=0;i<NUM_VERTICES;i++){
			vertexids[i] =reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
			vertexlids[i]=iomodel->my_vertices_lids[vertexids[i]-1];
		}

		/*Are we in transient or static? */
		if(M==1){
			values[0]=vector[0];
			//this->AddInput(vector_enum,values,P0Enum);
			_error_("not implemented yet");
		}
		else if(M==iomodel->numberofvertices){
			for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
			switch(this->ObjectEnum()){
				case TriaEnum:  inputs->SetTriaDatasetInput(enum_type,input_id,P1Enum,NUM_VERTICES,vertexlids,values); break;
				case PentaEnum: inputs->SetPentaDatasetInput(enum_type,input_id,P1Enum,NUM_VERTICES,vertexlids,values); break;
				default: _error_("Not implemented yet for "<<this->ObjectEnum());
			}
		}
		else if(M==iomodel->numberofvertices+1){
			/*create transient input: */
			IssmDouble* times = xNew<IssmDouble>(N);
			for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
			TransientInput* transientinput = inputs->SetDatasetTransientInput(enum_type,input_id,times,N);
			for(t=0;t<N;t++){
				for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
				switch(this->ObjectEnum()){
					case TriaEnum:  transientinput->AddTriaTimeInput( t,NUM_VERTICES,vertexlids,values,P1Enum); break;
					case PentaEnum: transientinput->AddPentaTimeInput(t,NUM_VERTICES,vertexlids,values,P1Enum); break;
					default: _error_("Not implemented yet");
				}
			}
			xDelete<IssmDouble>(times);
		}
		else{
			_error_("not implemented yet (M="<<M<<")");
		}

		xDelete<IssmDouble>(values);
		xDelete<int>(vertexids);
		xDelete<int>(vertexlids);
	}
	else if(vector_type==2){ //element vector
		_error_("not supported");

		IssmDouble value;

		/*Are we in transient or static? */
		if(M==iomodel->numberofelements){
			if (code==5){ //boolean
				_error_("not implemented");
				//datasetinput->AddInput(new BoolInput(input_enum,reCast<bool>(vector[this->Sid()])),input_id);
			}
			else if (code==6){ //integer
				_error_("not implemented");
				//datasetinput->AddInput(new IntInput(input_enum,reCast<int>(vector[this->Sid()])),input_id);
			}
			else if (code==7){ //IssmDouble
				_error_("not implemented");
				//datasetinput->AddInput(new DoubleInput(input_enum,vector[this->Sid()]),input_id);
			}
			else _error_("could not recognize nature of vector from code " << code);
		}
		else if(M==iomodel->numberofelements+1){
			_error_("not supported");
			///*create transient input: */
			//IssmDouble* times = xNew<IssmDouble>(N);
			//for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
			//TransientInput* transientinput=new TransientInput(input_enum,times,N);
			//TriaInput* bof=NULL;
			//for(t=0;t<N;t++){
			//	value=vector[N*this->Sid()+t];
			//	switch(this->ObjectEnum()){
			//		case TriaEnum:  transientinput->AddTimeInput(new TriaInput( input_enum,&value,P0Enum)); break;
			//		case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,&value,P0Enum)); break;
			//		case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,&value,P0Enum)); break;
			//		default: _error_("Not implemented yet");
			//	}
			//}
			//xDelete<IssmDouble>(times);
		}
		else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
	}
	else if(vector_type==3){ //element vector
		_error_("not supported");

		///*For right now we are static */
		//if(M==iomodel->numberofelements){
		//	/*create transient input: */
		//	IssmDouble* layers = xNewZeroInit<IssmDouble>(N);;
		//	for(t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
		//	DoubleArrayInput* arrayinput=new DoubleArrayInput(input_enum,layers,N);
		//	xDelete<IssmDouble>(layers);
		//}
		//else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
	}
	else{
		_error_("Cannot add input for vector type " << vector_type << " (not supported)");
	}
}
/*}}}*/
void       Element::InputUpdateFromConstant(int constant, int name){/*{{{*/

	/*Check that name is an element input*/
	if(!IsInputEnum(name)) _error_("Enum "<<EnumToStringx(name)<<" is not in IsInput");

	/*update input*/
	this->SetIntInput(this->inputs,name,constant);
}
/*}}}*/
void       Element::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/

	/*Check that name is an element input*/
	if(!IsInputEnum(name)) return;

	/*update input*/
	this->SetElementInput(name,constant);
}
/*}}}*/
void       Element::InputUpdateFromConstant(bool constant, int name){/*{{{*/

	/*Check that name is an element input*/
	if(!IsInputEnum(name)) return;

	/*update input*/
	this->SetBoolInput(this->inputs,name,constant);
}
/*}}}*/
bool       Element::IsOnSurface(){/*{{{*/
	return this->isonsurface;
}/*}}}*/
bool       Element::IsOnBase(){/*{{{*/
	return this->isonbase;
}/*}}}*/
bool       Element::IsFloating(){/*{{{*/

	bool shelf;
	int  migration_style;
	parameters->FindParam(&migration_style,GroundinglineMigrationEnum);

	Input* input = this->GetInput(MaskOceanLevelsetEnum); _assert_(input);

	if(migration_style==SubelementMigrationEnum){ //Floating if all nodes are floating
		if(input->GetInputMax() <= 0.) shelf=true;
		else shelf=false;
	}
	else if(migration_style==ContactEnum){
		if(input->GetInputMin() < 0.) shelf=true;
		else shelf=false;
	}
	else if(migration_style==NoneEnum || migration_style==AggressiveMigrationEnum || migration_style==SoftMigrationEnum || migration_style==GroundingOnlyEnum){ //Floating if all nodes are floating
		if(input->GetInputMin() > 0.) shelf=false;
		else shelf=true;
	}
	else _error_("migration_style not implemented yet");

	return shelf;
}/*}}}*/
bool       Element::IsGrounded(){/*{{{*/

	Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input);
	if(input->GetInputMax() > 0.){
		return true;
	}
	else{
		return false;
	}
}/*}}}*/
bool       Element::IsIceInElement(){/*{{{*/
	Input* input=this->GetInput(MaskIceLevelsetEnum); _assert_(input);
	return (input->GetInputMin()<0.);
}
/*}}}*/
bool       Element::IsIceOnlyInElement(){/*{{{*/
	Input* input=this->GetInput(MaskIceLevelsetEnum); _assert_(input);
	return (input->GetInputMax()<0.);
}
/*}}}*/
bool       Element::IsLandInElement(){/*{{{*/
	Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input);
	return (input->GetInputMax()>0.);
}
/*}}}*/
bool       Element::IsOceanInElement(){/*{{{*/
	Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input);
	return (input->GetInputMin()<0.);
}
/*}}}*/
void       Element::Ismip6FloatingiceMeltingRate(){/*{{{*/

	if(!this->IsIceInElement() || !this->IsFloating() || !this->IsOnBase()) return;

	int         basinid,num_basins,M,N;
	IssmDouble  tf,gamma0,base,delta_t_basin,mean_tf_basin,absval,meltanomaly;
	bool        islocal;
	IssmDouble* delta_t = NULL;
	IssmDouble* mean_tf = NULL;
	IssmDouble* depths  = NULL;

	/*Allocate some arrays*/
	const int numvertices = this->GetNumberOfVertices();
	IssmDouble* basalmeltrate = xNew<IssmDouble>(numvertices);

	/*Get variables*/
	IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum);
	IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum);
	IssmDouble lf   = this->FindParam(MaterialsLatentheatEnum);
	IssmDouble cp   = this->FindParam(MaterialsMixedLayerCapacityEnum);

	/*Hard code sea water density to be consistent with ISMIP6 documentation*/
	rhow = 1028.;

	/* Get parameters and inputs */
	this->GetInputValue(&basinid,BasalforcingsIsmip6BasinIdEnum);
	this->parameters->FindParam(&num_basins,BasalforcingsIsmip6NumBasinsEnum);
	this->parameters->FindParam(&gamma0,BasalforcingsIsmip6Gamma0Enum);
	this->parameters->FindParam(&delta_t,&M,BasalforcingsIsmip6DeltaTEnum);    _assert_(M==num_basins);
	this->parameters->FindParam(&islocal,BasalforcingsIsmip6IsLocalEnum);
	if(!islocal) {
		this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmip6AverageTfEnum); _assert_(N==num_basins);
	}
	Input* tf_input = this->GetInput(BasalforcingsIsmip6TfShelfEnum);              _assert_(tf_input);
	Input* meltanomaly_input = this->GetInput(BasalforcingsIsmip6MeltAnomalyEnum); _assert_(meltanomaly_input);
	delta_t_basin = delta_t[basinid];
	if(!islocal) mean_tf_basin = mean_tf[basinid];

	/*Compute melt rate for Local and Nonlocal parameterizations*/
	Gauss* gauss=this->NewGauss();
	for(int i=0;i<numvertices;i++){
		gauss->GaussVertex(i);
		tf_input->GetInputValue(&tf,gauss);
		meltanomaly_input->GetInputValue(&meltanomaly,gauss);
		if(!islocal) {
			absval = mean_tf_basin+delta_t_basin;
			if (absval<0) {absval = -1.*absval;}
			basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval + meltanomaly;
		}
		else{
			basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*pow(max(tf+delta_t_basin,0.),2) + meltanomaly;
		}
	}

	/*Return basal melt rate*/
	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrate,P1DGEnum);

	/*Cleanup and return*/
	delete gauss;
	xDelete<IssmDouble>(delta_t);
	xDelete<IssmDouble>(mean_tf);
	xDelete<IssmDouble>(depths);
	xDelete<IssmDouble>(basalmeltrate);

}/*}}}*/
void       Element::LinearFloatingiceMeltingRate(){/*{{{*/

	const int NUM_VERTICES = this->GetNumberOfVertices();

	IssmDouble  deepwaterel,upperwaterel,deepwatermelt,upperwatermelt;
	IssmDouble *base         = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble *values       = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble *perturbation = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble 	time;

	parameters->FindParam(&time,TimeEnum);
	parameters->FindParam(&deepwaterel,BasalforcingsDeepwaterElevationEnum,time);
	parameters->FindParam(&deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum,time);
	parameters->FindParam(&upperwaterel,BasalforcingsUpperwaterElevationEnum,time);
	parameters->FindParam(&upperwatermelt,BasalforcingsUpperwaterMeltingRateEnum,time);
	_assert_(upperwaterel>deepwaterel);

	this->GetInputListOnVertices(base,BaseEnum);
   this->GetInputListOnVertices(perturbation,BasalforcingsPerturbationMeltingRateEnum);
	for(int i=0;i<NUM_VERTICES;i++){
		if(base[i]>=upperwaterel){
			values[i]=upperwatermelt;
		}
		else if (base[i]<deepwaterel){
			values[i]=deepwatermelt;
		}
		else{
			IssmDouble alpha = (base[i]-upperwaterel)/(deepwaterel-upperwaterel);
			values[i]=deepwatermelt*alpha+(1.-alpha)*upperwatermelt;
		}

      values[i]+=perturbation[i];
	}

	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
	xDelete<IssmDouble>(base);
   xDelete<IssmDouble>(perturbation);
	xDelete<IssmDouble>(values);

}/*}}}*/
void       Element::SpatialLinearFloatingiceMeltingRate(){/*{{{*/

	const int NUM_VERTICES = this->GetNumberOfVertices();

	IssmDouble *deepwatermelt = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble *deepwaterel   = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble *upperwaterel  = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble *base          = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble *values        = xNew<IssmDouble>(NUM_VERTICES);

	this->GetInputListOnVertices(base,BaseEnum);
	this->GetInputListOnVertices(deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum);
	this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum);
	this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum);

	for(int i=0;i<NUM_VERTICES;i++){
		if(base[i]>upperwaterel[i])      values[i]=0;
		else if (base[i]<deepwaterel[i]) values[i]=deepwatermelt[i];
		else values[i]=deepwatermelt[i]*(base[i]-upperwaterel[i])/(deepwaterel[i]-upperwaterel[i]);
	}

	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
	xDelete<IssmDouble>(base);
	xDelete<IssmDouble>(deepwaterel);
	xDelete<IssmDouble>(deepwatermelt);
	xDelete<IssmDouble>(upperwaterel);
	xDelete<IssmDouble>(values);

}/*}}}*/
void       Element::MantlePlumeGeothermalFlux(){/*{{{*/

	const int NUM_VERTICES = this->GetNumberOfVertices();
	IssmDouble  mantleconductivity,nusselt,dtbg,plumeradius,topplumedepth,bottomplumedepth,plumex,plumey;
	IssmDouble  crustthickness,uppercrustthickness,uppercrustheat,lowercrustheat;
	IssmDouble  crustheat,plumeheat,dt,middleplumedepth,a,e,eprime,A0,lambda,Alambda,dAlambda;
	IssmDouble  x,y,z,c;
	IssmDouble* values   = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble *xyz_list = NULL;

	parameters->FindParam(&mantleconductivity,BasalforcingsMantleconductivityEnum);
	parameters->FindParam(&nusselt,BasalforcingsNusseltEnum);
	parameters->FindParam(&dtbg,BasalforcingsDtbgEnum);
	parameters->FindParam(&plumeradius,BasalforcingsPlumeradiusEnum);
	parameters->FindParam(&topplumedepth,BasalforcingsTopplumedepthEnum);
	parameters->FindParam(&bottomplumedepth,BasalforcingsBottomplumedepthEnum);
	parameters->FindParam(&plumex,BasalforcingsPlumexEnum);
	parameters->FindParam(&plumey,BasalforcingsPlumeyEnum);
	parameters->FindParam(&crustthickness,BasalforcingsCrustthicknessEnum);
	parameters->FindParam(&uppercrustthickness,BasalforcingsUppercrustthicknessEnum);
	parameters->FindParam(&uppercrustheat,BasalforcingsUppercrustheatEnum);
	parameters->FindParam(&lowercrustheat,BasalforcingsLowercrustheatEnum);

	this->GetVerticesCoordinates(&xyz_list);
	c=plumeradius;
	a=(bottomplumedepth-topplumedepth)/2.;
	e=pow(a*a-c*c,1./2.)/a;
	A0=(1-pow(e,2.))/pow(e,3.)*(1./2.*log((1+e)/(1-e))-e);
	for(int i=0;i<NUM_VERTICES;i++){
		y=xyz_list[i*3+0]-plumex;
		z=xyz_list[i*3+1]-plumey;
		x=-(a+topplumedepth+crustthickness);
		lambda=(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))/2;
		dAlambda=(-8*a*pow(c,2)*x*(-2*pow(a,2)+2*pow(c,2)+sqrt(2)*sqrt((a-c)*(a+c))*sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))*(pow(a,4)*(pow(y,2)+pow(z,2))+pow(c,4)*(pow(y,2)+pow(z,2))+pow(pow(x,2)+pow(y,2)+pow(z,2),2)*(pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))+pow(c,2)*(pow(x,4)-pow(x,2)*(pow(y,2)+pow(z,2))-(pow(y,2)+pow(z,2))*(2*pow(y,2)+2*pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))+pow(a,2)*(-pow(x,4)+pow(x,2)*(pow(y,2)+pow(z,2))+(pow(y,2)+pow(z,2))*(-2*pow(c,2)+2*(pow(y,2)+pow(z,2))+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))))))/(sqrt((a-c)*(a+c))*sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))*pow(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))),3.5)*pow(-(sqrt(2)*sqrt((a-c)*(a+c)))+sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2))))),2)*(sqrt(2)*sqrt((a-c)*(a+c))+sqrt(pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2)+sqrt(pow(-pow(a,2)-pow(c,2)+pow(x,2)+pow(y,2)+pow(z,2),2)+4*(pow(c,2)*pow(x,2)+pow(a,2)*(-pow(c,2)+pow(y,2)+pow(z,2)))))));
		eprime=pow((a*a-plumeradius*plumeradius)/(a*a+lambda),1./2.);
		Alambda=(1.-e*e)/(e*e*e)*(1./2.*log((1.+eprime)/(1.-eprime))-eprime);
		dt=dtbg-(nusselt-1.)/(1.+A0*(nusselt-1.))*(Alambda*dtbg+x*dtbg*dAlambda);
		plumeheat=mantleconductivity*dt;
		crustheat=uppercrustheat*uppercrustthickness+lowercrustheat*(crustthickness-uppercrustthickness);
		values[i]=crustheat+plumeheat;
	}

	this->AddInput(BasalforcingsGeothermalfluxEnum,values,P1Enum);
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(values);

}/*}}}*/
void       Element::MarshallElement2(MarshallHandle* marshallhandle,int numanalyses){/*{{{*/

	_assert_(this);
	if(marshallhandle->OperationNumber()==MARSHALLING_LOAD){
		this->nodes = NULL;
	}

	int object_enum = ElementEnum;
	marshallhandle->call(object_enum);

	marshallhandle->call(this->id);
	marshallhandle->call(this->sid);
	marshallhandle->call(this->lid);
	marshallhandle->call(this->element_type);
	marshallhandle->call(this->element_type_list,numanalyses);
}
/*}}}*/
void       Element::MigrateGroundingLine(IssmDouble* phi_ungrounding){/*{{{*/

	const int  NUM_VERTICES = this->GetNumberOfVertices();
	int        i,migration_style;
	IssmDouble bed_hydro,yts;
	IssmDouble rho_water,rho_ice,density;
	IssmDouble* melting = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* phi     = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* h       = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* s       = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* b       = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* r       = xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* sl      = xNew<IssmDouble>(NUM_VERTICES);

	/*Recover info at the vertices: */
	parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
	parameters->FindParam(&yts,ConstantsYtsEnum);
	GetInputListOnVertices(&h[0],ThicknessEnum);
	GetInputListOnVertices(&s[0],SurfaceEnum);
	GetInputListOnVertices(&b[0],BaseEnum);
	GetInputListOnVertices(&r[0],BedEnum);
	GetInputListOnVertices(&sl[0],SealevelEnum);
	GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum);
	rho_water   = FindParam(MaterialsRhoSeawaterEnum);
	rho_ice     = FindParam(MaterialsRhoIceEnum);
	density     = rho_ice/rho_water;

	/*go through vertices, and update inputs, considering them to be TriaVertex type: */
	for(i=0;i<NUM_VERTICES;i++){
		/* Contact FS*/
		if(migration_style == ContactEnum){
			phi[i]=phi_ungrounding[vertices[i]->Pid()];
			if(phi[i]>=0.) b[i]=r[i];
		}
		else if(migration_style == GroundingOnlyEnum && b[i]<r[i]) b[i]=r[i];
		/*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
		else if(phi[i]<=0.){
			if(b[i]<=r[i]){
				b[i]        = r[i];
				s[i]        = b[i]+h[i];
			}
		}
		/*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
		/*Change only if AggressiveMigration or if the ice sheet is in contact with the ocean*/
		else{ // phi>0
			bed_hydro=-density*h[i]+sl[i];
			if (bed_hydro>r[i]){
				/*Unground only if the element is connected to the ice shelf*/
				if(migration_style==AggressiveMigrationEnum || migration_style==SubelementMigrationEnum){
					s[i]        = (1-density)*h[i]+sl[i];
					b[i]        = -density*h[i]+sl[i];
				}
				else if(migration_style==SoftMigrationEnum && phi_ungrounding[vertices[i]->Pid()]<0.){
					s[i]        = (1-density)*h[i]+sl[i];
					b[i]        = -density*h[i]+sl[i];
				}
				else{
					if(migration_style!=SoftMigrationEnum && migration_style!=ContactEnum && migration_style!=GroundingOnlyEnum) _error_("Error: migration should be Aggressive, Soft, Subelement, Contact or GroundingOnly");
				}
			}
		}
	}

	/*Recalculate phi*/
	for(i=0;i<NUM_VERTICES;i++){
		if(migration_style==SoftMigrationEnum){
			bed_hydro=-density*h[i]+sl[i];
			if(phi[i]<0. || bed_hydro<=r[i] || phi_ungrounding[vertices[i]->Pid()]<0.){
				phi[i]=h[i]+(r[i]-sl[i])/density;
			}
		}
		else if(migration_style!=ContactEnum) phi[i]=h[i]+(r[i]-sl[i])/density;
		else{
			/*do nothing*/
		}
	}
	this->AddInput(MaskOceanLevelsetEnum,&phi[0],P1Enum);

	/*Update inputs*/
	this->AddInput(SurfaceEnum,&s[0],P1Enum);
	this->AddInput(BaseEnum,&b[0],P1Enum);

	/*Delete*/
	xDelete<IssmDouble>(melting);
	xDelete<IssmDouble>(phi);
	xDelete<IssmDouble>(r);
	xDelete<IssmDouble>(b);
	xDelete<IssmDouble>(s);
	xDelete<IssmDouble>(sl);
	xDelete<IssmDouble>(h);

}
/*}}}*/
void       Element::MismipFloatingiceMeltingRate(){/*{{{*/

	IssmDouble thresholdthickness,upperdepthmelt;
	IssmDouble base[MAXVERTICES];
	IssmDouble bed[MAXVERTICES];
	IssmDouble meltratefactor[MAXVERTICES];
	IssmDouble values[MAXVERTICES];

	const int NUM_VERTICES = this->GetNumberOfVertices();

	parameters->FindParam(&thresholdthickness,BasalforcingsThresholdThicknessEnum);
	parameters->FindParam(&upperdepthmelt,BasalforcingsUpperdepthMeltEnum);
	this->GetInputListOnVertices(base,BaseEnum);
	this->GetInputListOnVertices(bed,BedEnum);
	this->GetInputListOnVertices(meltratefactor,BasalforcingsMeltrateFactorEnum);
	for(int i=0;i<NUM_VERTICES;i++){
		if(base[i]>upperdepthmelt){
			values[i]=0;
		}
		else{
			values[i]=meltratefactor[i]*tanh((base[i]-bed[i])/thresholdthickness)*(upperdepthmelt-base[i]);
		}
	}

	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);

}/*}}}*/
void       Element::BeckmannGoosseFloatingiceMeltingRate(){/*{{{*/

	int numvertices      = this->GetNumberOfVertices();
	IssmDouble T_f,ocean_heat_flux;
	IssmDouble rho_water    = this->FindParam(MaterialsRhoSeawaterEnum);
	IssmDouble rho_ice      = this->FindParam(MaterialsRhoIceEnum);
	IssmDouble latentheat   = this->FindParam(MaterialsLatentheatEnum);
	IssmDouble mixed_layer_capacity = this->FindParam(MaterialsMixedLayerCapacityEnum);
	IssmDouble thermal_exchange_vel = this->FindParam(MaterialsThermalExchangeVelocityEnum);

	IssmDouble base[MAXVERTICES];
	IssmDouble values[MAXVERTICES];
	IssmDouble oceansalinity [MAXVERTICES];
	IssmDouble oceantemp[MAXVERTICES];
	IssmDouble meltratefactor[MAXVERTICES];

	this->GetInputListOnVertices(base,BaseEnum);
	this->GetInputListOnVertices(oceansalinity,BasalforcingsOceanSalinityEnum);
	this->GetInputListOnVertices(oceantemp,BasalforcingsOceanTempEnum);
	this->GetInputListOnVertices(meltratefactor,BasalforcingsMeltrateFactorEnum);

	Gauss* gauss=this->NewGauss();
	for(int i=0;i<numvertices;i++){
		T_f=(0.0939 - 0.057 * oceansalinity[i] + 7.64e-4 * base[i]); //degC

		// compute ocean_heat_flux according to beckmann_goosse2003
		// positive, if T_oc > T_ice ==> heat flux FROM ocean TO ice
		ocean_heat_flux = meltratefactor[i] * rho_water * mixed_layer_capacity * thermal_exchange_vel * (oceantemp[i] - T_f); // in W/m^2

		// shelfbmassflux is positive if ice is freezing on; here it is always negative:
		// same sign as ocean_heat_flux (positive if massflux FROM ice TO ocean)
		values[i] = ocean_heat_flux / (latentheat * rho_ice); // m s-1
	}

	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
	delete gauss;
}/*}}}*/
void       Element::MungsmtpParameterization(void){/*{{{*/
	/*Are we on the base? If not, return*/
	if(!IsOnBase()) return;

	const int NUM_VERTICES 					= this->GetNumberOfVertices();
	const int NUM_VERTICES_MONTHS_PER_YEAR 	= NUM_VERTICES * 12;

	int i;
	int*        vertexlids=xNew<int>(NUM_VERTICES);
	IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* TemperaturesLgm=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* PrecipitationsLgm=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble TdiffTime,PfacTime;

	/*Recover parameters*/
	IssmDouble time,yts,time_yr;
	this->parameters->FindParam(&time,TimeEnum);
	this->parameters->FindParam(&yts,ConstantsYtsEnum);
	this->GetVerticesLidList(vertexlids);
	time_yr=floor(time/yts)*yts;

	/*Recover present day temperature and precipitation*/
	DatasetInput* dinput1=this->GetDatasetInput(SmbTemperaturesPresentdayEnum);   _assert_(dinput1);
	DatasetInput* dinput2=this->GetDatasetInput(SmbTemperaturesLgmEnum);          _assert_(dinput2);
	DatasetInput* dinput3=this->GetDatasetInput(SmbPrecipitationsPresentdayEnum); _assert_(dinput3);
	DatasetInput* dinput4=this->GetDatasetInput(SmbPrecipitationsLgmEnum);        _assert_(dinput4);

	/*loop over vertices: */
	Gauss* gauss=this->NewGauss();
	for(int month=0;month<12;month++) {
		for(int iv=0;iv<NUM_VERTICES;iv++) {
			gauss->GaussVertex(iv);
			dinput1->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month);
			dinput2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss,month);
			dinput3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month);
			dinput4->GetInputValue(&PrecipitationsLgm[iv*12+month],gauss,month);

			PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts;
			PrecipitationsLgm[iv*12+month]=PrecipitationsLgm[iv*12+month]*yts;
		}
	}

	/*Recover interpolation parameters at time t*/
	this->parameters->FindParam(&TdiffTime,SmbTdiffEnum,time);
	this->parameters->FindParam(&PfacTime,SmbPfacEnum,time);

	/*Compute the temperature and precipitation*/
	for(int iv=0;iv<NUM_VERTICES;iv++){
		ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
					&PrecipitationsLgm[iv*12],&PrecipitationsPresentday[iv*12],
					&TemperaturesLgm[iv*12], &TemperaturesPresentday[iv*12],
					&monthlytemperatures[iv*12], &monthlyprec[iv*12]);
	}

	/*Update inputs*/
	for (int imonth=0;imonth<12;imonth++) {
		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
		switch(this->ObjectEnum()){
			case TriaEnum:  this->inputs->SetTriaDatasetInput(SmbMonthlytemperaturesEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			case PentaEnum: this->inputs->SetPentaDatasetInput(SmbMonthlytemperaturesEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			default: _error_("Not implemented yet");
		}
		for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
		switch(this->ObjectEnum()){
			case TriaEnum:  this->inputs->SetTriaDatasetInput(SmbPrecipitationEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			case PentaEnum: this->inputs->SetPentaDatasetInput(SmbPrecipitationEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
			default: _error_("Not implemented yet");
		}
	}

	switch(this->ObjectEnum()){
		case TriaEnum: break;
		case PentaEnum:
		case TetraEnum:
							this->DatasetInputExtrude(SmbMonthlytemperaturesEnum,-1);
							this->DatasetInputExtrude(SmbPrecipitationEnum,-1);
							break;
		default: _error_("Not implemented yet");
	}

	/*clean-up*/
	delete gauss;
	xDelete<IssmDouble>(monthlytemperatures);
	xDelete<IssmDouble>(monthlyprec);
	xDelete<IssmDouble>(TemperaturesPresentday);
	xDelete<IssmDouble>(TemperaturesLgm);
	xDelete<IssmDouble>(PrecipitationsPresentday);
	xDelete<IssmDouble>(PrecipitationsLgm);
	xDelete<IssmDouble>(tmp);
	xDelete<int>(vertexlids);

}
/*}}}*/
ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/
	return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
}
/*}}}*/
ElementMatrix* Element::NewElementMatrixCoupling(int number_nodes,int approximation_enum){/*{{{*/
	return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum);
}
/*}}}*/
ElementVector* Element::NewElementVector(int approximation_enum){/*{{{*/
	return new ElementVector(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
}
/*}}}*/
void       Element::PicoUpdateBoxid(int* max_boxid_basin_list){/*{{{*/

	if(!this->IsIceInElement() || !this->IsFloating()) return;

	int        basin_id;
	IssmDouble dist_gl,dist_cf;

	this->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum);
	IssmDouble boxid_max=reCast<IssmDouble>(max_boxid_basin_list[basin_id])+1.;

	Input* dist_gl_input=this->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
	Input* dist_cf_input=this->GetInput(DistanceToCalvingfrontEnum);  _assert_(dist_cf_input);

	/*Get dist_gl and dist_cf at center of element*/
	Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
	dist_gl_input->GetInputValue(&dist_gl,gauss);
	dist_cf_input->GetInputValue(&dist_cf,gauss);
	delete gauss;

	/*Ensure values are positive for floating ice*/
	dist_gl = fabs(dist_gl);
	dist_cf = fabs(dist_cf);

	/*Compute relative distance to grounding line*/
	IssmDouble rel_dist_gl=dist_gl/(dist_gl+dist_cf);

	/*Assign box numbers based on rel_dist_gl*/
	int boxid = -1;
	for(IssmDouble i=0.;i<boxid_max;i++){
		IssmDouble lowbound  = 1. -sqrt((boxid_max-i   )/boxid_max);
		IssmDouble highbound = 1. -sqrt((boxid_max-i-1.)/boxid_max);
		if(rel_dist_gl>=lowbound && rel_dist_gl<=highbound){
			boxid=reCast<int>(i);
			break;
		}
	}
	if(boxid==-1) _error_("No boxid found for element " << this->Sid() << "!");

	this->SetIntInput(this->inputs,BasalforcingsPicoBoxIdEnum, boxid);

}/*}}}*/
void       Element::PicoUpdateBox(int loopboxid){/*{{{*/

	if(!this->IsIceInElement() || !this->IsFloating()) return;

	int boxid;
	this->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
	if(loopboxid!=boxid) return;

	const int NUM_VERTICES = this->GetNumberOfVertices();

	int        basinid, maxbox, num_basins, numnodes, M;
	IssmDouble gamma_T, overturning_coeff, thickness;
	IssmDouble pressure, T_star,p_coeff, q_coeff;
	bool       isplume;

	/*Get variables*/
	IssmDouble rhoi       = this->FindParam(MaterialsRhoIceEnum);
	IssmDouble rhow       = this->FindParam(MaterialsRhoSeawaterEnum);
	IssmDouble earth_grav = this->FindParam(ConstantsGEnum);
	IssmDouble rho_star   = 1033.;             // kg/m^3
	IssmDouble nu         = rhoi/rhow;
	IssmDouble latentheat = this->FindParam(MaterialsLatentheatEnum);
	IssmDouble c_p_ocean  = this->FindParam(MaterialsMixedLayerCapacityEnum);
	IssmDouble lambda     = latentheat/c_p_ocean;
	IssmDouble a          = -0.0572;          // K/psu
	IssmDouble b          = 0.0788 + this->FindParam(MaterialsMeltingpointEnum);  //K
	IssmDouble c          = 7.77e-4;
	IssmDouble alpha      = 7.5e-5;           // 1/K
	IssmDouble Beta       = 7.7e-4;           // K

	/* Get non-box-specific parameters and inputs */
	this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum);
	this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
	this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
	this->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
	this->parameters->FindParam(&isplume, BasalforcingsPicoIsplumeEnum);
	Input *thickness_input    = this->GetInput(ThicknessEnum);                         _assert_(thickness_input);
	_assert_(basinid<=num_basins);

	IssmDouble* boxareas = NULL;
	this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
	_assert_(M==num_basins*maxbox);

	IssmDouble area_boxi        = boxareas[basinid*maxbox+boxid];
	IssmDouble g1               = area_boxi*gamma_T;

	IssmDouble* basalmeltrates_shelf 				= xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* potential_pressure_melting_point	= xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* Tocs 								= xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* Socs								= xNew<IssmDouble>(NUM_VERTICES);

	/* First box calculations */
	if(boxid==0){
		/* Get box1 parameters and inputs */
		IssmDouble time, toc_farocean, soc_farocean;
		this->parameters->FindParam(&time,TimeEnum);
		this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum);
		this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);
		IssmDouble 	s1 				= soc_farocean/(nu*lambda);
		IssmDouble* overturnings 	= xNew<IssmDouble>(NUM_VERTICES);
		Input *overturningC_input = this->GetInput(BasalforcingsPicoOverturningCoeffEnum); _assert_(overturningC_input);

		/* Start looping on the number of verticies and calculate ocean vars */
		Gauss* gauss=this->NewGauss();
		for(int i=0;i<NUM_VERTICES;i++){
			gauss->GaussVertex(i);
			thickness_input->GetInputValue(&thickness,gauss);
			overturningC_input->GetInputValue(&overturning_coeff,gauss);
			pressure = (rhoi*earth_grav*1e-4)*thickness;
			T_star   = a*soc_farocean+b-c*pressure-toc_farocean;
			p_coeff  = g1/(overturning_coeff*rho_star*(Beta*s1-alpha));
			q_coeff  = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha)));

			/* To avoid negatives under the square root */
			if((0.25*pow(p_coeff,2)-q_coeff)<0) q_coeff = 0.25*p_coeff*p_coeff;

			Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff));
			Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]);
			potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
			if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
			overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i]));
		}

		if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1DGEnum);
		this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1DGEnum);
		this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1DGEnum);
		this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1DGEnum);

		/*Cleanup and return*/
		delete gauss;
	}

	/* Subsequent box calculations */
	else {
		/* Get subsequent box parameters and inputs */
		IssmDouble* toc_weighted_avg         = NULL;
		IssmDouble* soc_weighted_avg         = NULL;
		IssmDouble* overturning_weighted_avg = NULL;
		this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum);
		this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum);
		this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum);
		IssmDouble mean_toc                  = toc_weighted_avg[basinid];
		IssmDouble mean_soc                  = soc_weighted_avg[basinid];
		IssmDouble mean_overturning          = overturning_weighted_avg[basinid];
		IssmDouble g2                        = g1/(nu*lambda);

		/* Start looping on the number of verticies and calculate ocean vars */
		Gauss* gauss=this->NewGauss();
		for(int i=0;i<NUM_VERTICES;i++){
			gauss->GaussVertex(i);
			thickness_input->GetInputValue(&thickness,gauss);
			pressure = (rhoi*earth_grav*1.e-4)*thickness;
			T_star   = a*mean_soc+b-c*pressure-mean_toc;
			Tocs[i]  = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc));
			Socs[i]  = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda));
			potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
			if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
		}

		if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1DGEnum);
		this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1DGEnum);
		this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1DGEnum);

		/*Cleanup and return*/
		xDelete<IssmDouble>(toc_weighted_avg);
		xDelete<IssmDouble>(soc_weighted_avg);
		xDelete<IssmDouble>(overturning_weighted_avg);
		delete gauss;
	}

	/*Cleanup and return*/
	xDelete<IssmDouble>(boxareas);

}/*}}}*/
void       Element::PicoComputeBasalMelt(){/*{{{*/

	if(!this->IsIceInElement() || !this->IsFloating()) return;

	const int NUM_VERTICES = this->GetNumberOfVertices();

	IssmDouble E0, Cd, CdT, YT, lam1, lam2, lam3, M0, CdTS0, y1, y2, x0;
	IssmDouble Tf_gl, YTS, CdTS, G1, G2, G3, g_alpha, M, l, X_hat, M_hat;
	IssmDouble alpha, zgl, Toc, Soc, z_base, yts, slopex, slopey;

	/*Get variables*/
	E0    = 3.6e-2;        //Entrainment coefficient
	Cd    = 2.5e-3;        //Drag coefficient
	CdT   = 1.1e-3;        //Turbulent heat exchange coefficient
	YT    = CdT/sqrt(Cd);  //Heat exchange coefficient
	lam1  = -5.73e-2;      //Freezing point-salinity coefficient (degrees C)
	lam2  = 8.32e-2;       //Freezing point offset (degrees C)
	lam3  = 7.61e-4;       //Freezing point-depth coefficient (K m-1)
	M0    = 10.;           //Melt-rate parameter (m yr-1 C-2)
	CdTS0 = 6e-4;          //Heat exchange parameter
	y1    = 0.545;         //Heat exchange parameter
	y2    = 3.5e-5;        //Heat exchange parameter
	x0    = 0.56;          //Dimentionless scaling factor

	/*Define arrays*/
	IssmDouble* basalmeltrates_shelf = xNew<IssmDouble>(NUM_VERTICES); //Basal melt-rate

	/*Polynomial coefficients*/
	IssmDouble p[12];
	p[0]  =  0.1371330075095435;
	p[1]  =  5.527656234709359E1;
	p[2]  = -8.951812433987858E2;
	p[3]  =  8.927093637594877E3;
	p[4]  = -5.563863123811898E4;
	p[5]  =  2.218596970948727E5;
	p[6]  = -5.820015295669482E5;
	p[7]  =  1.015475347943186E6;
	p[8]  = -1.166290429178556E6;
	p[9]  =  8.466870335320488E5;
	p[10] = -3.520598035764990E5;
	p[11] =  6.387953795485420E4;

	/*Get inputs*/
	Input* zgl_input         = this->GetInput(GroundinglineHeightEnum);                     _assert_(zgl_input);
	Input* toc_input         = this->GetInput(BasalforcingsPicoSubShelfOceanTempEnum);      _assert_(toc_input);
	Input* soc_input         = this->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum);  _assert_(soc_input);
	Input* base_input        = this->GetInput(BaseEnum);                                    _assert_(base_input);
	Input* baseslopex_input  = this->GetInput(BaseSlopeXEnum);                              _assert_(baseslopex_input);
	Input* baseslopey_input  = this->GetInput(BaseSlopeYEnum);                              _assert_(baseslopey_input);
	this->FindParam(&yts, ConstantsYtsEnum);

	/*Loop over the number of vertices in this element*/
	Gauss* gauss=this->NewGauss();
	for(int i=0;i<NUM_VERTICES;i++){
		gauss->GaussVertex(i);

		/*Get inputs*/
		zgl_input->GetInputValue(&zgl,gauss);
		toc_input->GetInputValue(&Toc,gauss); //(K)
		soc_input->GetInputValue(&Soc,gauss);
		base_input->GetInputValue(&z_base,gauss);
		baseslopex_input->GetInputValue(&slopex,gauss);
		baseslopey_input->GetInputValue(&slopey,gauss);

		/*Compute ice shelf base slope (radians)*/
		alpha = atan(sqrt(slopex*slopex + slopey*slopey));
		if(alpha>=M_PI) alpha = M_PI - 0.001;               //ensure sin(alpha) > 0 for meltrate calculations

		/*Make necessary conversions*/
		Toc = Toc-273.15;
		if(zgl>z_base) zgl=z_base;

		/*Low bound for Toc to ensure X_hat is between 0 and 1*/
		if(Toc<lam1*Soc+lam2) Toc=lam1*Soc+lam2;

		/*Compute parameters needed for melt-rate calculation*/
		Tf_gl = lam1*Soc+lam2+lam3*zgl;                                              //Characteristic freezing point
		YTS = YT*(y1+y2*(((Toc-Tf_gl)*E0*sin(alpha))/(lam3*(CdTS0+E0*sin(alpha))))); //Effective heat exchange coefficient
		CdTS = sqrt(Cd)*YTS;                                                         //Heat exchange coefficient
		G1 = sqrt(sin(alpha)/(Cd+E0*sin(alpha)));                                    //Geometric factor
		G2 = sqrt(CdTS/(CdTS+E0*sin(alpha)));                                        //Geometric factor
		G3 = (E0*sin(alpha))/(CdTS+E0*sin(alpha));                                   //Geometric factor
		g_alpha = G1*G2*G3;                                                          //Melt scaling factor
		M = M0*g_alpha*pow((Toc-Tf_gl),2);                                           //Melt-rate scale
		l = ((Toc-Tf_gl)*(x0*CdTS+E0*sin(alpha)))/(lam3*x0*(CdTS+E0*sin(alpha)));    //Length scale
		X_hat = (z_base-zgl)/l;                                                      //Dimentionless coordinate system

		/*Compute polynomial fit*/
		M_hat = 0.;                                                                  //Reset summation variable for each node
		for(int ii=0;ii<12;ii++) {
			M_hat += p[ii]*pow(X_hat,ii);                                               //Polynomial fit
		}

		/*Compute melt-rate*/
		basalmeltrates_shelf[i] = (M*M_hat)/yts;                                     //Basal melt-rate (m/s)
	}

	/*Save computed melt-rate*/
	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1DGEnum);

	/*Cleanup and return*/
	delete gauss;

}/*}}}*/
void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/

	const int NUM_VERTICES 		= this->GetNumberOfVertices();
	const int NUM_VERTICES_MONTHS_PER_YEAR = NUM_VERTICES * 12;

	int  		i,vertexlids[MAXVERTICES];
	IssmDouble* agd=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
	IssmDouble* melt=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
	IssmDouble* accu=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
	IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* yearlytemperatures=xNew<IssmDouble>(NUM_VERTICES); memset(yearlytemperatures, 0., NUM_VERTICES*sizeof(IssmDouble));
	IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* h=xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* s0p=xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* s0t=xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble pddsnowfac = -1.;
	IssmDouble pddicefac = -1.;
	IssmDouble rho_water,rho_ice,desfac,rlaps,rlapslgm;
	IssmDouble PfacTime,TdiffTime,sealevTime;
	IssmDouble mavg=1./12.; //factor for monthly average

	/*Get vertex Lids for later*/
	this->GetVerticesLidList(&vertexlids[0]);

	/*Get material parameters :*/
	rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
	rho_ice=this->FindParam(MaterialsRhoIceEnum);

	/*Get some pdd parameters*/
	desfac=this->FindParam(SmbDesfacEnum);
	rlaps=this->FindParam(SmbRlapsEnum);
	rlapslgm=this->FindParam(SmbRlapslgmEnum);

	IssmDouble time,yts,time_yr;
	this->parameters->FindParam(&time,TimeEnum);
	this->parameters->FindParam(&yts,ConstantsYtsEnum);
	time_yr=floor(time/yts)*yts;

	/*Get inputs*/
	DatasetInput* dinput =this->GetDatasetInput(SmbMonthlytemperaturesEnum); _assert_(dinput);
	DatasetInput* dinput2=this->GetDatasetInput(SmbPrecipitationEnum);       _assert_(dinput2);

	/*loop over vertices: */
	Gauss* gauss=this->NewGauss();
	for(int month=0;month<12;month++) {
		/*Recover monthly temperatures and precipitation and compute the yearly mean temperatures*/

		for(int iv=0;iv<NUM_VERTICES;iv++) {
			gauss->GaussVertex(iv);
			dinput->GetInputValue(&monthlytemperatures[iv*12+month],gauss,month);
			// yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv*12+month]*mavg; // Has to be in Kelvin
			monthlytemperatures[iv*12+month]=monthlytemperatures[iv*12+month]-273.15; // conversion from Kelvin to celcius for PDD module
			dinput2->GetInputValue(&monthlyprec[iv*12+month],gauss,month);
			monthlyprec[iv*12+month]=monthlyprec[iv*12+month]*yts;
		}
	}

	/*Recover Pfac, Tdiff and sealev at time t:
	 *     This parameters are used to interpolate the temperature
	 *         and precipitaton between PD and LGM when ismungsm==1 */
	if (ismungsm==1){
		this->parameters->FindParam(&TdiffTime,SmbTdiffEnum,time);
		this->parameters->FindParam(&sealevTime,SmbSealevEnum,time);
	}
	else {
		TdiffTime=0;
		sealevTime=0;
	}

	/*Recover pdd factors at time t.
	 *     This parameter is set, if the user wants to define the
	 *     pdd factors regionally, if issetpddfac==1 in the d18opdd method */
	Input* input  = NULL;
	Input* input2 = NULL;
	if(issetpddfac==1){
		input  = this->GetInput(SmbPddfacSnowEnum); _assert_(input);
		input2 = this->GetInput(SmbPddfacIceEnum); _assert_(input2);
	}

	/*Recover info at the vertices: */
	GetInputListOnVertices(&h[0],ThicknessEnum);
	GetInputListOnVertices(&s[0],SurfaceEnum);
	GetInputListOnVertices(&s0p[0],SmbS0pEnum);
	GetInputListOnVertices(&s0t[0],SmbS0tEnum);

	/*measure the surface mass balance*/
	for(int iv = 0; iv<NUM_VERTICES; iv++){
		gauss->GaussVertex(iv);
		pddsnowfac=0.;
		pddicefac=0.;
		if(issetpddfac==1){
			input->GetInputValue(&pddsnowfac,gauss);
			input2->GetInputValue(&pddicefac,gauss);
		}
		agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv*12], &monthlyprec[iv*12],
					pdds, pds, &melt[iv], &accu[iv], signorm, yts, h[iv], s[iv],
					desfac, s0t[iv], s0p[iv],rlaps,rlapslgm,TdiffTime,sealevTime,
					pddsnowfac,pddicefac,rho_water,rho_ice);
		/*Get yearlytemperatures */
		for(int month=0;month<12;month++) {
			yearlytemperatures[iv]=yearlytemperatures[iv]+(monthlytemperatures[iv*12+month]+273.15)*mavg; // Has to be in Kelvin
		}
	}

	/*Update inputs*/
	switch(this->ObjectEnum()){
		case TriaEnum:
			this->AddInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum);
			this->AddInput(SmbMassBalanceEnum,&agd[0],P1Enum);
			this->AddInput(SmbAccumulationEnum,&accu[0],P1Enum);
			this->AddInput(SmbMeltEnum,&melt[0],P1Enum);
			break;
		case PentaEnum:
			bool isthermal;
         this->parameters->FindParam(&isthermal,TransientIsthermalEnum);
         if(isthermal){
				if(IsOnSurface()){
					/*Here, we want to change the BC of the thermal model, keep
					 * the temperatures as they are for the base of the penta and
					 * yse yearlytemperatures for the top*/
					PentaInput* temp_input = xDynamicCast<PentaInput*>(this->GetInput(TemperatureEnum)); _assert_(temp_input);
					switch(temp_input->GetInputInterpolationType()){
						case P1Enum:
							temp_input->element_values[3] = yearlytemperatures[3];
							temp_input->element_values[4] = yearlytemperatures[4];
							temp_input->element_values[5] = yearlytemperatures[5];
							temp_input->SetInput(P1Enum,NUM_VERTICES,&vertexlids[0],temp_input->element_values);
							break;
						case P1xP2Enum:
						case P1xP3Enum:
						case P1xP4Enum:
							temp_input->element_values[3] = yearlytemperatures[3];
							temp_input->element_values[4] = yearlytemperatures[4];
							temp_input->element_values[5] = yearlytemperatures[5];
							temp_input->SetInput(temp_input->GetInputInterpolationType(),this->lid,this->GetNumberOfNodes(temp_input->GetInputInterpolationType()),temp_input->element_values);
							break;
						default:
							_error_("Interpolation "<<EnumToStringx(temp_input->GetInputInterpolationType())<<" not supported yet");
					}

					bool isenthalpy;
					this->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
					if(isenthalpy){
						/*Convert that to enthalpy for the enthalpy model*/
						PentaInput* enth_input = xDynamicCast<PentaInput*>(this->GetInput(EnthalpyEnum)); _assert_(enth_input);
						switch(enth_input->GetInputInterpolationType()){
							case P1Enum:
								ThermalToEnthalpy(&enth_input->element_values[3],yearlytemperatures[3],0.,0.);
								ThermalToEnthalpy(&enth_input->element_values[4],yearlytemperatures[4],0.,0.);
								ThermalToEnthalpy(&enth_input->element_values[5],yearlytemperatures[5],0.,0.);
								enth_input->SetInput(P1Enum,NUM_VERTICES,&vertexlids[0],enth_input->element_values);
								break;
							case P1xP2Enum:
							case P1xP3Enum:
							case P1xP4Enum:
								ThermalToEnthalpy(&enth_input->element_values[3],yearlytemperatures[3],0.,0.);
								ThermalToEnthalpy(&enth_input->element_values[4],yearlytemperatures[4],0.,0.);
								ThermalToEnthalpy(&enth_input->element_values[5],yearlytemperatures[5],0.,0.);
								enth_input->SetInput(enth_input->GetInputInterpolationType(),this->lid,this->GetNumberOfNodes(enth_input->GetInputInterpolationType()),enth_input->element_values);
								break;
							default:
								_error_("Interpolation "<<EnumToStringx(temp_input->GetInputInterpolationType())<<" not supported yet");
						}
					}
				}
			}
			this->AddInput(SmbMassBalanceEnum,&agd[0],P1Enum);
			this->AddInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum);
			this->AddInput(SmbAccumulationEnum,&accu[0],P1Enum);
			this->AddInput(SmbMeltEnum,&melt[0],P1Enum);
			this->InputExtrude(TemperaturePDDEnum,-1);
			this->InputExtrude(SmbMassBalanceEnum,-1);
			this->InputExtrude(SmbAccumulationEnum,-1);
			this->InputExtrude(SmbMeltEnum,-1);
			break;
		default: _error_("Not implemented yet");
	}

	/*clean-up*/
	delete gauss;
	xDelete<IssmDouble>(monthlytemperatures);
	xDelete<IssmDouble>(monthlyprec);
	xDelete<IssmDouble>(agd);
	xDelete<IssmDouble>(melt);
	xDelete<IssmDouble>(accu);
	xDelete<IssmDouble>(yearlytemperatures);
	xDelete<IssmDouble>(h);
	xDelete<IssmDouble>(s);
	xDelete<IssmDouble>(s0t);
	xDelete<IssmDouble>(s0p);
	xDelete<IssmDouble>(tmp);
}
/*}}}*/
void       Element::PositiveDegreeDaySicopolis(bool isfirnwarming){/*{{{*/

	/* General FIXMEs: get Tmelting point, pddicefactor, pddsnowfactor, sigma from parameters/user input */

	const int NUM_VERTICES 		= this->GetNumberOfVertices();
	const int NUM_VERTICES_MONTHS_PER_YEAR	= NUM_VERTICES * 12;

	int        	i,vertexlids[MAXVERTICES];;
	IssmDouble* smb=xNew<IssmDouble>(NUM_VERTICES);		// surface mass balance
	IssmDouble* melt=xNew<IssmDouble>(NUM_VERTICES);		// melting comp. of surface mass balance
	IssmDouble* accu=xNew<IssmDouble>(NUM_VERTICES);		// accuumulation comp. of surface mass balance
	IssmDouble* melt_star=xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
	IssmDouble* yearlytemperatures=xNew<IssmDouble>(NUM_VERTICES); memset(yearlytemperatures, 0., NUM_VERTICES*sizeof(IssmDouble));
	IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);			// actual surface height
	IssmDouble* s0p=xNew<IssmDouble>(NUM_VERTICES);		// reference elevation for precip.
	IssmDouble* s0t=xNew<IssmDouble>(NUM_VERTICES);		// reference elevation for temperature
	IssmDouble* smbcorr=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance correction; will be added after pdd call
	IssmDouble* p_ampl=xNew<IssmDouble>(NUM_VERTICES);	// precip anomaly
	IssmDouble* t_ampl=xNew<IssmDouble>(NUM_VERTICES);	// remperature anomaly
	IssmDouble rho_water,rho_ice,desfac,rlaps;
	IssmDouble inv_twelve=1./12.;								//factor for monthly average
	IssmDouble time,yts,time_yr;

	/*Get vertex Lids for later*/
	this->GetVerticesLidList(&vertexlids[0]);

	/*Get material parameters :*/
	rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
	rho_ice=this->FindParam(MaterialsRhoIceEnum);

	/*Get parameters for height corrections*/
	desfac=this->FindParam(SmbDesfacEnum);
	rlaps=this->FindParam(SmbRlapsEnum);

	/* Get time */
	this->parameters->FindParam(&time,TimeEnum);
	this->parameters->FindParam(&yts,ConstantsYtsEnum);
	time_yr=floor(time/yts)*yts;

	/* Set parameters for finrnwarming */
	IssmDouble MU_0         = 9.7155; //Firn-warming correction, in (d*deg C)/(mm WE)
	IssmDouble mu           = MU_0*(1000.0*86400.0)*(rho_ice/rho_water);   // (d*deg C)/(mm WE) --> (s*deg C)/(m IE)

	/*Get inputs*/
	DatasetInput* dinput =this->GetDatasetInput(SmbMonthlytemperaturesEnum); _assert_(dinput);
	DatasetInput* dinput2=this->GetDatasetInput(SmbPrecipitationEnum);       _assert_(dinput2);

	/*loop over vertices: */
	Gauss* gauss=this->NewGauss();
	for(int month=0;month<12;month++){

		for(int iv=0;iv<NUM_VERTICES;iv++){
			gauss->GaussVertex(iv);
			dinput->GetInputValue(&monthlytemperatures[iv*12+month],gauss,month);
			monthlytemperatures[iv*12+month]=monthlytemperatures[iv*12+month]-273.15; // conversion from Kelvin to celcius for PDD module
			dinput2->GetInputValue(&monthlyprec[iv*12+month],gauss,month);
			monthlyprec[iv*12+month]=monthlyprec[iv*12+month]*yts;
		}
	}

	/*Recover info at the vertices: */
	GetInputListOnVertices(&s[0],SurfaceEnum);
	GetInputListOnVertices(&s0p[0],SmbS0pEnum);
	GetInputListOnVertices(&s0t[0],SmbS0tEnum);
	GetInputListOnVertices(&smbcorr[0],SmbSmbCorrEnum);
	GetInputListOnVertices(&t_ampl[0],SmbTemperaturesAnomalyEnum);
	GetInputListOnVertices(&p_ampl[0],SmbPrecipitationsAnomalyEnum);

	/*measure the surface mass balance*/
	for (int iv = 0; iv<NUM_VERTICES; iv++){
		smb[iv]=PddSurfaceMassBalanceSicopolis(&monthlytemperatures[iv*12], &monthlyprec[iv*12],
					&melt[iv], &accu[iv], &melt_star[iv], &t_ampl[iv], &p_ampl[iv], yts, s[iv],
					desfac, s0t[iv], s0p[iv],rlaps,rho_water,rho_ice);

		/* make correction */
		smb[iv] = smb[iv]+smbcorr[iv];
		/*Get yearlytemperatures */
		for(int month=0;month<12;month++) yearlytemperatures[iv]=yearlytemperatures[iv]+((monthlytemperatures[iv*12+month]+273.15)+t_ampl[iv])*inv_twelve; // Has to be in Kelvin

		if(isfirnwarming){
			if(melt_star[iv]>=melt[iv]){
				yearlytemperatures[iv]= yearlytemperatures[iv]+mu*(melt_star[iv]-melt[iv]);
			}
			else{
				yearlytemperatures[iv]= yearlytemperatures[iv];
			}
		}
		if (yearlytemperatures[iv]>273.15) yearlytemperatures[iv]=273.15;
	}

	switch(this->ObjectEnum()){
		case TriaEnum:
			//this->AddInput(TemperatureEnum,&yearlytemperatures[0],P1Enum);
			this->AddInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum);
			this->AddInput(SmbMassBalanceEnum,&smb[0],P1Enum);
			this->AddInput(SmbAccumulationEnum,&accu[0],P1Enum);
			this->AddInput(SmbMeltEnum,&melt[0],P1Enum);
			break;
		case PentaEnum:
			bool isthermal;
			this->parameters->FindParam(&isthermal,TransientIsthermalEnum);
			if(isthermal){
				bool isenthalpy;
				this->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
				if(IsOnSurface()){
					/*Here, we want to change the BC of the thermal model, keep
					 * the temperatures as they are for the base of the penta and
					 * use yearlytemperatures for the top*/

					/*FIXME: look at other function Element::PositiveDegreeDay and propagate change! Just assert for now*/
					PentaInput* temp_input = xDynamicCast<PentaInput*>(this->GetInput(TemperatureEnum)); _assert_(temp_input);
					switch(temp_input->GetInputInterpolationType()){
						case P1Enum:
							temp_input->element_values[3] = yearlytemperatures[3];
							temp_input->element_values[4] = yearlytemperatures[4];
							temp_input->element_values[5] = yearlytemperatures[5];
							temp_input->SetInput(P1Enum,NUM_VERTICES,&vertexlids[0],temp_input->element_values);
							break;
						case P1DGEnum:
						case P1xP2Enum:
						case P1xP3Enum:
						case P1xP4Enum:
							temp_input->element_values[3] = yearlytemperatures[3];
							temp_input->element_values[4] = yearlytemperatures[4];
							temp_input->element_values[5] = yearlytemperatures[5];
							temp_input->SetInput(temp_input->GetInputInterpolationType(),this->lid,this->GetNumberOfNodes(temp_input->GetInputInterpolationType()),temp_input->element_values);
							break;
						default:
							_error_("Interpolation "<<EnumToStringx(temp_input->GetInputInterpolationType())<<" not supported yet");
					}

					if(isenthalpy){
						/*Convert that to enthalpy for the enthalpy model*/
						PentaInput* enth_input = xDynamicCast<PentaInput*>(this->GetInput(EnthalpyEnum)); _assert_(enth_input);
						switch(enth_input->GetInputInterpolationType()){
							case P1Enum:
								ThermalToEnthalpy(&enth_input->element_values[3],yearlytemperatures[3],0.,0.);
								ThermalToEnthalpy(&enth_input->element_values[4],yearlytemperatures[4],0.,0.);
								ThermalToEnthalpy(&enth_input->element_values[5],yearlytemperatures[5],0.,0.);
								enth_input->SetInput(P1Enum,NUM_VERTICES,&vertexlids[0],enth_input->element_values);
								break;
							case P1DGEnum:
							case P1xP2Enum:
							case P1xP3Enum:
							case P1xP4Enum:
								ThermalToEnthalpy(&enth_input->element_values[3],yearlytemperatures[3],0.,0.);
								ThermalToEnthalpy(&enth_input->element_values[4],yearlytemperatures[4],0.,0.);
								ThermalToEnthalpy(&enth_input->element_values[5],yearlytemperatures[5],0.,0.);
								enth_input->SetInput(enth_input->GetInputInterpolationType(),this->lid,this->GetNumberOfNodes(enth_input->GetInputInterpolationType()),enth_input->element_values);
								break;
							default:
								_error_("Interpolation "<<EnumToStringx(temp_input->GetInputInterpolationType())<<" not supported yet");
						}
					}
				}
			}
			this->AddInput(SmbMassBalanceEnum,&smb[0],P1Enum);
			this->AddInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum);
			this->AddInput(SmbAccumulationEnum,&accu[0],P1Enum);
			this->AddInput(SmbMeltEnum,&melt[0],P1Enum);
			this->InputExtrude(TemperaturePDDEnum,-1);
			this->InputExtrude(SmbMassBalanceEnum,-1);
			this->InputExtrude(SmbAccumulationEnum,-1);
			this->InputExtrude(SmbMeltEnum,-1);
			break;
		default: _error_("Not implemented yet");
	}

	/*clean-up*/
	delete gauss;
	xDelete<IssmDouble>(monthlytemperatures);
	xDelete<IssmDouble>(monthlyprec);
	xDelete<IssmDouble>(smb);
	xDelete<IssmDouble>(melt);
	xDelete<IssmDouble>(accu);
	xDelete<IssmDouble>(yearlytemperatures);
	xDelete<IssmDouble>(s);
	xDelete<IssmDouble>(s0t);
	xDelete<IssmDouble>(s0p);
	xDelete<IssmDouble>(t_ampl);
	xDelete<IssmDouble>(p_ampl);
	xDelete<IssmDouble>(smbcorr);
	xDelete<IssmDouble>(melt_star);
}
/*}}}*/
void       Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int* parray_size, int output_enum){/*{{{*/

	/*Some intputs need to be computed, even if they are already in inputs, they might not be up to date!*/
	switch(output_enum){
		case ViscousHeatingEnum: this->ViscousHeatingCreateInput(); break;
		case StressMaxPrincipalEnum: this->StressMaxPrincipalCreateInput(); break;
		case StressTensorxxEnum:
		case StressTensorxyEnum:
		case StressTensorxzEnum:
		case StressTensoryyEnum:
		case StressTensoryzEnum:
		case StressTensorzzEnum: this->ComputeStressTensor(); break;
		case StrainRatexxEnum:
		case StrainRatexyEnum:
		case StrainRatexzEnum:
		case StrainRateyyEnum:
		case StrainRateyzEnum:
		case StrainRatezzEnum:
		case StrainRateeffectiveEnum: this->ComputeStrainRate(); break;
		case DeviatoricStressxxEnum:
		case DeviatoricStressxyEnum:
		case DeviatoricStressxzEnum:
		case DeviatoricStressyyEnum:
		case DeviatoricStressyzEnum:
		case DeviatoricStresszzEnum:
		case DeviatoricStress1Enum:
		case DeviatoricStress2Enum:
		case DeviatoricStresseffectiveEnum: this->ComputeDeviatoricStressTensor(); break;
		case EsaStrainratexxEnum:
		case EsaStrainratexyEnum:
		case EsaStrainrateyyEnum:
		case EsaRotationrateEnum: this->ComputeEsaStrainAndVorticity(); break;
		case SigmaNNEnum: this->ComputeSigmaNN(); break;
		case LambdaSEnum: this->ComputeLambdaS(); break;
		case StressIntensityFactorEnum: this->StressIntensityFactor(); break;
		case CalvingratexEnum:
		case CalvingrateyEnum:
		case CalvingCalvingrateEnum:
												  this->StrainRateparallel();
												  this->StrainRateperpendicular();
												  int calvinglaw;
												  this->FindParam(&calvinglaw,CalvingLawEnum);
												  switch(calvinglaw){
													  case DefaultCalvingEnum:
														  //do nothing
														  break;
													  case CalvingLevermannEnum:
														  this->CalvingRateLevermann();
														  break;
													  case CalvingVonmisesEnum:
													  case CalvingDev2Enum:
														  this->CalvingRateVonmises();
														  break;
													  case CalvingCrevasseDepthEnum:
														  this->CalvingCrevasseDepth();
														  break;
													  default:
														  _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
												  }
												  break;
		case CalvingFluxLevelsetEnum: this->CalvingFluxLevelset(); break;
		case CalvingMeltingFluxLevelsetEnum: this->CalvingMeltingFluxLevelset(); break;
		case StrainRateparallelEnum: this->StrainRateparallel(); break;
		case StrainRateperpendicularEnum: this->StrainRateperpendicular(); break;
		case SurfaceCrevasseEnum: this->CalvingCrevasseDepth(); break;
		case SigmaVMEnum: this->CalvingRateVonmises(); break;
		case PartitioningEnum: this->inputs->SetInput(PartitioningEnum,this->lid,IssmComm::GetRank()); break;
	}

	/*If this input is not already in Inputs, maybe it needs to be computed?*/
	switch(this->inputs->GetInputObjectEnum(output_enum)){
		case TriaInputEnum:
		case PentaInputEnum:
		case ControlInputEnum:
		case TransientInputEnum:{
			Input* input2 = this->GetInput(output_enum);
			if(!input2) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
			*pinterpolation   = input2->GetResultInterpolation();
			*pnodesperelement = input2->GetResultNumberOfNodes();
			*parray_size      = input2->GetResultArraySize();
			}
			break;
		case BoolInputEnum:
			*pinterpolation   = P0Enum;
			*pnodesperelement = 1;
			*parray_size      = 1;
			break;
		case IntInputEnum:
			*pinterpolation   = P0Enum;
			*pnodesperelement = 1;
			*parray_size      = 1;
			break;
		case ArrayInputEnum:{
			int M;
			this->inputs->GetArray(output_enum,this->lid,NULL,&M);
			*pinterpolation   = P0ArrayEnum;
			*pnodesperelement = 1;
			*parray_size      = M;
			}
			break;
		default:
			_error_("Input type \""<<EnumToStringx(this->inputs->GetInputObjectEnum(output_enum))<<"\" not supported yet (While trying to return "<<EnumToStringx(output_enum)<<")");
	}


	/*Assign output pointer*/

	return;
}/*}}}*/
void       Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/

	/*Find input*/
	Input* input=this->GetInput(output_enum);
	if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");

	/*Cast to ElementInput*/
	if(input->ObjectEnum()!=TriaInputEnum && input->ObjectEnum()!=PentaInputEnum){
		_error_("Input "<<EnumToStringx(output_enum)<<" is not an ElementInput");
	}
	ElementInput* element_input = xDynamicCast<ElementInput*>(input);

	/*Get Number of nodes and make sure that it is the same as the one provided*/
	int numnodes = this->GetNumberOfNodes(element_input->GetInputInterpolationType());
	_assert_(numnodes==nodesperelement);

	/*Fill in arrays*/
	for(int i=0;i<numnodes;i++) values[this->sid*numnodes + i] = element_input->element_values[i];

} /*}}}*/
void       Element::ResultToMatrix(IssmDouble* values,int ncols,int output_enum){/*{{{*/

	IssmDouble* array = NULL;
	int         m;
	this->inputs->GetArray(output_enum,this->lid,&array,&m);
	for(int i=0;i<m;i++) values[this->Sid()*ncols + i] = array[i];
	xDelete<IssmDouble>(array);

} /*}}}*/
void       Element::ResultToVector(Vector<IssmDouble>* vector,int output_enum){/*{{{*/

	IssmDouble values[MAXVERTICES];
	int        connectivity[MAXVERTICES];
	int        sidlist[MAXVERTICES];

	switch(this->inputs->GetInputObjectEnum(output_enum)){
		case TriaInputEnum:
		case PentaInputEnum:
		case ControlInputEnum:
		case TransientInputEnum:{

			Input* input2 = this->GetInput(output_enum);
			if(!input2) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");

			switch(input2->GetResultInterpolation()){
				case P0Enum:{
					IssmDouble  value;
					bool        bvalue;
					Gauss* gauss = this->NewGauss();
					input2->GetInputValue(&value,gauss);
					delete gauss;
					vector->SetValue(this->Sid(),value,INS_VAL);
					break;
					}
				case P1Enum:{
					const int NUM_VERTICES = this->GetNumberOfVertices();



					this->GetVerticesSidList(&sidlist[0]);
					this->GetVerticesConnectivityList(&connectivity[0]);
					this->GetInputListOnVertices(&values[0],output_enum);
					for(int i=0;i<NUM_VERTICES;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);
					vector->SetValues(NUM_VERTICES,sidlist,values,ADD_VAL);
					break;
					}
				default:
					_error_("interpolation "<<EnumToStringx(input2->GetResultInterpolation())<<" not supported yet");
				}
			}
			break;
		case BoolInputEnum:
			bool bvalue;
			this->GetInputValue(&bvalue,output_enum);
			vector->SetValue(this->Sid(),reCast<IssmDouble>(bvalue),INS_VAL);
			break;
		case IntInputEnum:
			int ivalue;
			this->GetInputValue(&ivalue,output_enum);
			vector->SetValue(this->Sid(),reCast<IssmDouble>(ivalue),INS_VAL);
			break;
		default:
			_error_("Input type \""<<EnumToStringx(this->inputs->GetInputObjectEnum(output_enum))<<"\" not supported yet");
	}

} /*}}}*/
void       Element::SetBoolInput(Inputs* inputs,int enum_in,bool value){/*{{{*/

	_assert_(inputs);
	inputs->SetInput(enum_in,this->lid,value);

}
/*}}}*/
void       Element::SetIntInput(Inputs* inputs,int enum_in,int value){/*{{{*/

	_assert_(inputs);
	inputs->SetInput(enum_in,this->lid,value);

}
/*}}}*/
void       Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int* flagsindices_counter,int set1_enum,int set2_enum, int analysis_type){/*{{{*/

	/*Intermediaries*/
	const int numnodes = this->GetNumberOfNodes();
	/*Output */
	int d_nz = 0;
	int o_nz = 0;

	/*Loop over all nodes*/
	for(int i=0;i<numnodes;i++){
		int nodelid = this->nodes[i]->Lid();
		if(!flags[nodelid]){
			/*flag current node so that no other element processes it*/
			flags[nodelid]=true;

			flagsindices[flagsindices_counter[0]]=nodelid;
			flagsindices_counter[0]++;

			/*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
			switch(set2_enum){
				case FsetEnum:
					if(nodes[i]->fsize){
						if(this->nodes[i]->IsClone())
						 o_nz += 1;
						else
						 d_nz += 1;
					}
					break;
				case GsetEnum:
					if(nodes[i]->gsize){
						if(this->nodes[i]->IsClone())
						 o_nz += 1;
						else
						 d_nz += 1;
					}
					break;
				case SsetEnum:
					if(nodes[i]->ssize){
						if(this->nodes[i]->IsClone())
						 o_nz += 1;
						else
						 d_nz += 1;
					}
					break;
				default: _error_("not supported");
			}
		}
	}

	/*Special case: 2d/3d coupling, the node of this element might be connected
	 *to the basal element*/
	int approximation,numlayers;
	if(analysis_type==StressbalanceAnalysisEnum){
		this->GetInputValue(&approximation,ApproximationEnum);
		if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){
			parameters->FindParam(&numlayers,MeshNumberoflayersEnum);
			o_nz += numlayers*3;
			d_nz += numlayers*3;
		}
	}

	/*Assign output pointers: */
	*pd_nz=d_nz;
	*po_nz=o_nz;
}
/*}}}*/
#ifdef _HAVE_SEMIC_
void       Element::SmbSemic(){/*{{{*/

	/*only compute SMB at the surface: */
	if (!IsOnSurface()) return;

	const int NUM_VERTICES 					= this->GetNumberOfVertices();
	const int NUM_VERTICES_DAYS_PER_YEAR 	= NUM_VERTICES * 365;

	IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* s0gcm=xNew<IssmDouble>(NUM_VERTICES);
	IssmDouble* st=xNew<IssmDouble>(NUM_VERTICES);

	// daily forcing inputs
	IssmDouble* dailyrainfall=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
	IssmDouble* dailysnowfall=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
	IssmDouble* dailydlradiation=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
	IssmDouble* dailydsradiation=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
	IssmDouble* dailywindspeed=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
	IssmDouble* dailypressure=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
	IssmDouble* dailyairdensity=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
	IssmDouble* dailyairhumidity=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
	IssmDouble* dailytemperature=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
	// daily outputs
	IssmDouble* tsurf_out=xNew<IssmDouble>(NUM_VERTICES); memset(tsurf_out, 0., NUM_VERTICES*sizeof(IssmDouble));
	IssmDouble* smb_out=xNew<IssmDouble>(NUM_VERTICES); memset(smb_out, 0., NUM_VERTICES*sizeof(IssmDouble));
	IssmDouble* saccu_out=xNew<IssmDouble>(NUM_VERTICES); memset(saccu_out, 0., NUM_VERTICES*sizeof(IssmDouble));
	IssmDouble* smelt_out=xNew<IssmDouble>(NUM_VERTICES); memset(smelt_out, 0., NUM_VERTICES*sizeof(IssmDouble));

	IssmDouble rho_water,rho_ice,desfac,rlaps,rdl;
	IssmDouble time,yts,time_yr;

	/* Get time: */
	this->parameters->FindParam(&time,TimeEnum);
	this->parameters->FindParam(&yts,ConstantsYtsEnum);
	time_yr=floor(time/yts)*yts;

	/*Get material parameters :*/
	rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
	rho_ice=this->FindParam(MaterialsRhoIceEnum);
	desfac=this->FindParam(SmbDesfacEnum);
	rlaps=this->FindParam(SmbRlapsEnum);
	rdl=this->FindParam(SmbRdlEnum);

	/* Recover info at the vertices: */
	GetInputListOnVertices(&s[0],SurfaceEnum);
	GetInputListOnVertices(&s0gcm[0],SmbS0gcmEnum);

	/* loop over vertices and days */ //FIXME account for leap years (365 -> 366)
	Gauss* gauss=this->NewGauss();
	for (int iday = 0; iday < 365; iday++){
		/* Retrieve inputs: */
		Input* dailysnowfall_input    = this->GetInput(SmbDailysnowfallEnum,time_yr+(iday+1)/365.*yts); _assert_(dailysnowfall_input);
		Input* dailyrainfall_input    = this->GetInput(SmbDailyrainfallEnum,time_yr+(iday+1)/365.*yts); _assert_(dailyrainfall_input);
		Input* dailydlradiation_input = this->GetInput(SmbDailydlradiationEnum,time_yr+(iday+1)/365.*yts); _assert_(dailydlradiation_input);
		Input* dailydsradiation_input = this->GetInput(SmbDailydsradiationEnum,time_yr+(iday+1)/365.*yts); _assert_(dailydsradiation_input);
		Input* dailywindspeed_input   = this->GetInput(SmbDailywindspeedEnum,time_yr+(iday+1)/365.*yts); _assert_(dailywindspeed_input);
		Input* dailypressure_input    = this->GetInput(SmbDailypressureEnum,time_yr+(iday+1)/365.*yts); _assert_(dailypressure_input);
		Input* dailyairdensity_input  = this->GetInput(SmbDailyairdensityEnum,time_yr+(iday+1)/365.*yts); _assert_(dailyairdensity_input);
		Input* dailyairhumidity_input = this->GetInput(SmbDailyairhumidityEnum,time_yr+(iday+1)/365.*yts); _assert_(dailyairhumidity_input);
		Input* dailytemperature_input = this->GetInput(SmbDailytemperatureEnum,time_yr+(iday+1)/365.*yts); _assert_(dailytemperature_input);

		for(int iv=0;iv<NUM_VERTICES;iv++){
			gauss->GaussVertex(iv);
			/* get forcing */
			dailyrainfall_input->GetInputValue(&dailyrainfall[iv*365+iday],gauss);
			dailysnowfall_input->GetInputValue(&dailysnowfall[iv*365+iday],gauss);
			dailydlradiation_input->GetInputValue(&dailydlradiation[iv*365+iday],gauss);
			dailydsradiation_input->GetInputValue(&dailydsradiation[iv*365+iday],gauss);
			dailywindspeed_input->GetInputValue(&dailywindspeed[iv*365+iday],gauss);
			dailypressure_input->GetInputValue(&dailypressure[iv*365+iday],gauss);
			dailyairdensity_input->GetInputValue(&dailyairdensity[iv*365+iday],gauss);
			dailyairhumidity_input->GetInputValue(&dailyairhumidity[iv*365+iday],gauss);
			dailytemperature_input->GetInputValue(&dailytemperature[iv*365+iday],gauss);

			/* Surface temperature correction */
			st[iv]=(s[iv]-s0gcm[iv])/1000.;
			dailytemperature[iv*365+iday]=dailytemperature[iv*365+iday]-rlaps *st[iv];

			/* Precipitation correction (Vizcaino et al. 2010) */
			if (s0gcm[iv] < 2000.0) {
				dailysnowfall[iv*365+iday] = dailysnowfall[iv*365+iday]*exp(desfac*(max(s[iv],2000.0)-2000.0));
				dailyrainfall[iv*365+iday] = dailyrainfall[iv*365+iday]*exp(desfac*(max(s[iv],2000.0)-2000.0));
			}else{
				dailysnowfall[iv*365+iday] = dailysnowfall[iv*365+iday]*exp(desfac*(max(s[iv],2000.0)-s0gcm[iv]));
				dailyrainfall[iv*365+iday] = dailyrainfall[iv*365+iday]*exp(desfac*(max(s[iv],2000.0)-s0gcm[iv]));
			}

			/* downward longwave radiation correction (Marty et al. 2002) */
			st[iv]=(s[iv]-s0gcm[iv])/1000.;
			dailydlradiation[iv*365+iday]=dailydlradiation[iv*365+iday]+rdl*st[iv];
		}
	}

	for (int iv = 0; iv<NUM_VERTICES; iv++){
		/* call semic */
		run_semic_(&dailysnowfall[iv*365], &dailyrainfall[iv*365], &dailydsradiation[iv*365], &dailydlradiation[iv*365],
					&dailywindspeed[iv*365], &dailypressure[iv*365], &dailyairdensity[iv*365], &dailyairhumidity[iv*365], &dailytemperature[iv*365],
					&tsurf_out[iv], &smb_out[iv], &saccu_out[iv], &smelt_out[iv]);
	}

	switch(this->ObjectEnum()){
		case TriaEnum:
			this->AddInput(TemperatureSEMICEnum,&tsurf_out[0],P1Enum); // TODO add TemperatureSEMICEnum to EnumDefinitions
			this->AddInput(SmbMassBalanceEnum,&smb_out[0],P1Enum);
			this->AddInput(SmbAccumulationEnum,&saccu_out[0],P1Enum);
			this->AddInput(SmbMeltEnum,&smelt_out[0],P1Enum);
			break;
		case PentaEnum:
			// TODO
			break;
		case TetraEnum:
			// TODO
			break;
		default: _error_("Not implemented yet");
	}

	/*clean-up*/
	delete gauss;
	xDelete<IssmDouble>(dailysnowfall);
	xDelete<IssmDouble>(dailyrainfall);
	xDelete<IssmDouble>(dailydlradiation);
	xDelete<IssmDouble>(dailydsradiation);
	xDelete<IssmDouble>(dailywindspeed);
	xDelete<IssmDouble>(dailypressure);
	xDelete<IssmDouble>(dailyairdensity);
	xDelete<IssmDouble>(dailyairhumidity);
	xDelete<IssmDouble>(dailypressure);
	xDelete<IssmDouble>(dailytemperature);
	xDelete<IssmDouble>(smb_out);
	xDelete<IssmDouble>(smelt_out);
	xDelete<IssmDouble>(saccu_out);
	xDelete<IssmDouble>(tsurf_out);
	xDelete<IssmDouble>(s);
	xDelete<IssmDouble>(st);
	xDelete<IssmDouble>(s0gcm);
}
/*}}}*/
#endif // _HAVE_SEMIC_
int        Element::Sid(){/*{{{*/

	return this->sid;

}
/*}}}*/
void       Element::SmbGemb(IssmDouble timeinputs, int count){/*{{{*/

	/*only compute SMB at the surface: */
	if (!IsOnSurface()) return;

	/*Intermediary variables: {{{*/
	bool       isinitialized;
	IssmDouble zTop=0.0;
	IssmDouble dzTop=0.0;
	IssmDouble zMax=0.0;
	IssmDouble zMin=0.0;
	IssmDouble zY=0.0;
	IssmDouble dzMin=0.0;
	IssmDouble Tmean=0.0;
	IssmDouble Vmean=0.0;
	IssmDouble C=0.0;
	IssmDouble Tz,Vz=0.0;
	IssmDouble yts;
	IssmDouble Ta=0.0;
	IssmDouble V=0.0;
	IssmDouble dlw=0.0;
	IssmDouble dsw=0.0;
	IssmDouble P=0.0;
	IssmDouble eAir=0.0;
	IssmDouble pAir=0.0;
	IssmDouble teValue=1.0;
	IssmDouble aValue=0.0;
	IssmDouble dt,time,smb_dt;
	int        aIdx=0;
	int        denIdx=0;
	int        dsnowIdx=0;
	int        swIdx=0;
	IssmDouble cldFrac,t0wet, t0dry, K;
	IssmDouble lhf=0.0;
	IssmDouble shf=0.0;
	IssmDouble dayEC=0.0;
	IssmDouble initMass=0.0;
   IssmDouble sumR=0.0;
   IssmDouble sumM=0.0;
	IssmDouble sumMsurf=0.0;
   IssmDouble sumEC=0.0;
   IssmDouble sumP=0.0;
   IssmDouble sumW=0.0;
   IssmDouble sumMassAdd=0.0;
	IssmDouble fac=0.0;
	IssmDouble sumMass=0.0;
	IssmDouble dMass=0.0;
	IssmDouble accsumR=0.0;
	IssmDouble accsumM=0.0;
	IssmDouble accsumSMB=0.0;
	IssmDouble accsumEC=0.0;
	IssmDouble accsumP=0.0;
	bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux;
	bool isclimatology=false;
	bool isconstrainsurfaceT=false;
	IssmDouble init_scaling=0.0;
	IssmDouble thermo_scaling=1.0;
	IssmDouble adThresh=1023.0;
	/*}}}*/
	/*Output variables:{{{ */
	IssmDouble* dz=NULL;
	IssmDouble* d = NULL;
	IssmDouble* re = NULL;
	IssmDouble* gdn = NULL;
	IssmDouble* gsp = NULL;
	IssmDouble  EC = 0.0;
	IssmDouble  ulw = 0.0;
	IssmDouble  netSW=0.0;
	IssmDouble  netLW=0.0;
	IssmDouble  meanULW=0.0;
	IssmDouble  meanLHF=0.0;
	IssmDouble  meanSHF=0.0;
	IssmDouble* W = NULL;
	IssmDouble* a = NULL;
	IssmDouble* swf=NULL;
	IssmDouble* T = NULL;
	IssmDouble  T_bottom = 0.0;
	IssmDouble  M = 0.0;
	IssmDouble  Msurf = 0.0;
	IssmDouble  R = 0.0;
	IssmDouble  mAdd = 0.0;
	IssmDouble  dz_add = 0.0;
	IssmDouble* dzini=NULL;
	IssmDouble* dini = NULL;
	IssmDouble* reini = NULL;
	IssmDouble* gdnini = NULL;
	IssmDouble* gspini = NULL;
	IssmDouble* Wini = NULL;
	IssmDouble* aini = NULL;
	IssmDouble* Tini = NULL;
	int         m=0;
	/*}}}*/

	/*Retrieve material properties and parameters:{{{ */
	IssmDouble rho_ice   = FindParam(MaterialsRhoIceEnum);
	IssmDouble rho_water = FindParam(MaterialsRhoFreshwaterEnum);
	IssmDouble aSnow     = parameters->FindParam(SmbASnowEnum);
	IssmDouble aIce      = parameters->FindParam(SmbAIceEnum);
	parameters->FindParam(&time,TimeEnum);                        /*transient core time at which we run the smb core*/
	parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
	parameters->FindParam(&yts,ConstantsYtsEnum);
	parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/
	parameters->FindParam(&aIdx,SmbAIdxEnum);
	parameters->FindParam(&denIdx,SmbDenIdxEnum);
	parameters->FindParam(&swIdx,SmbSwIdxEnum);
	parameters->FindParam(&dsnowIdx,SmbDsnowIdxEnum);
	parameters->FindParam(&cldFrac,SmbCldFracEnum);
	parameters->FindParam(&t0wet,SmbT0wetEnum);
	parameters->FindParam(&t0dry,SmbT0dryEnum);
	parameters->FindParam(&K,SmbKEnum);
	parameters->FindParam(&isgraingrowth,SmbIsgraingrowthEnum);
	parameters->FindParam(&isalbedo,SmbIsalbedoEnum);
	parameters->FindParam(&isshortwave,SmbIsshortwaveEnum);
	parameters->FindParam(&isthermal,SmbIsthermalEnum);
	parameters->FindParam(&isaccumulation,SmbIsaccumulationEnum);
	parameters->FindParam(&ismelt,SmbIsmeltEnum);
	parameters->FindParam(&isdensification,SmbIsdensificationEnum);
	parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum);
	parameters->FindParam(&isconstrainsurfaceT,SmbIsconstrainsurfaceTEnum);
	parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum);
	parameters->FindParam(&thermo_scaling,SmbThermoDeltaTScalingEnum);
	parameters->FindParam(&adThresh,SmbAdThreshEnum);
	/*}}}*/
	/*Retrieve inputs: {{{*/
	Input *zTop_input          = this->GetInput(SmbZTopEnum);         _assert_(zTop_input);
	Input *dzTop_input         = this->GetInput(SmbDzTopEnum);        _assert_(dzTop_input);
	Input *dzMin_input         = this->GetInput(SmbDzMinEnum);        _assert_(dzMin_input);
	Input *zMax_input          = this->GetInput(SmbZMaxEnum);         _assert_(zMax_input);
	Input *zMin_input          = this->GetInput(SmbZMinEnum);         _assert_(zMin_input);
	Input *zY_input            = this->GetInput(SmbZYEnum);           _assert_(zY_input);
	Input *Tmean_input         = this->GetInput(SmbTmeanEnum);        _assert_(Tmean_input);
	Input *Vmean_input         = this->GetInput(SmbVmeanEnum);        _assert_(Vmean_input);
	Input *C_input             = this->GetInput(SmbCEnum);            _assert_(C_input);
	Input *Tz_input            = this->GetInput(SmbTzEnum);           _assert_(Tz_input);
	Input *Vz_input            = this->GetInput(SmbVzEnum);           _assert_(Vz_input);
	Input *EC_input            = NULL;

	/*Retrieve input values:*/
	Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);

	this->GetInputValue(&isinitialized,SmbIsInitializedEnum);
	zTop_input->GetInputValue(&zTop,gauss);
	dzTop_input->GetInputValue(&dzTop,gauss);
	dzMin_input->GetInputValue(&dzMin,gauss);
	zMax_input->GetInputValue(&zMax,gauss);
	zMin_input->GetInputValue(&zMin,gauss);
	zY_input->GetInputValue(&zY,gauss);
	Tmean_input->GetInputValue(&Tmean,gauss);
	Vmean_input->GetInputValue(&Vmean,gauss);
	C_input->GetInputValue(&C,gauss);
	Tz_input->GetInputValue(&Tz,gauss);
	Vz_input->GetInputValue(&Vz,gauss);
	/*}}}*/

	/*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/
	if(!isinitialized){
		if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n");
		//if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" <<
		//dz[i] << "\n");

		this->inputs->GetArray(SmbDziniEnum,this->lid,&dzini,&m);
		this->inputs->GetArray(SmbDiniEnum,this->lid,&dini,&m);
		this->inputs->GetArray(SmbReiniEnum,this->lid,&reini,&m);
		this->inputs->GetArray(SmbGdniniEnum,this->lid,&gdnini,&m);
		this->inputs->GetArray(SmbGspiniEnum,this->lid,&gspini,&m);
		this->inputs->GetArray(SmbWiniEnum,this->lid,&Wini,&m);
		this->inputs->GetArray(SmbAiniEnum,this->lid,&aini,&m);
		this->inputs->GetArray(SmbTiniEnum,this->lid,&Tini,&m);
		EC_input = this->GetInput(SmbECiniEnum);  _assert_(EC_input);
		EC_input->GetInputAverage(&EC);

		/*Retrieve the correct value of m (without the zeroes at the end)*/
		this->GetInputValue(&m,SmbSizeiniEnum);

		if(m==2){ //Snow properties are initialized with default values. Vertical grid has to be initialized too
			//            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w DEFAULT values\n");

			/*initialize profile variables:*/
			GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);

			d = xNew<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=dini[0]; //ice density [kg m-3]
			re = xNew<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=reini[0];         //set grain size to old snow [mm]
			gdn = xNew<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=gdnini[0];         //set grain dentricity to old snow
			gsp = xNew<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=gspini[0];         //set grain sphericity to old snow
			W = xNew<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0];             //set water content to zero [kg m-2]
			a = xNew<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0];         //set albedo equal to fresh snow [fraction]
			T = xNew<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
			/*/!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties)
			 *    because don't know Tmean yet when set default values.
			 *    Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/

			//fixed lower temperature bounday condition - T is fixed
			T_bottom=T[m-1];
		}
		else{ //Retrieve snow properties from previous run. Need to provide values for all layers
			//            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n");

			dz = xNew<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzini[i];
			d = xNew<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=dini[i];
			re = xNew<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=reini[i];
			gdn = xNew<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnini[i];
			gsp = xNew<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gspini[i];
			W = xNew<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wini[i];
			a = xNew<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=aini[i];
			T = xNew<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Tini[i];

			//fixed lower temperature bounday condition - T is fixed
			_assert_(m>0);
			T_bottom=T[m-1];
		}

		this->SetElementInput(SmbAccumulatedECEnum,0.0);
		this->SetElementInput(SmbAccumulatedMassBalanceEnum,0.0);
		this->SetElementInput(SmbAccumulatedMeltEnum,0.0);
		this->SetElementInput(SmbAccumulatedRunoffEnum,0.0);
		this->SetElementInput(SmbAccumulatedPrecipitationEnum,0.0);

		/*Flag the initialization:*/
		this->SetBoolInput(this->inputs,SmbIsInitializedEnum,true);
	}
	else{
		/*Recover inputs: */
		this->inputs->GetArray(SmbDzEnum,this->lid,&dz,&m);
		this->inputs->GetArray(SmbDEnum,this->lid,&d,&m);
		this->inputs->GetArray(SmbReEnum,this->lid,&re,&m);
		this->inputs->GetArray(SmbGdnEnum,this->lid,&gdn,&m);
		this->inputs->GetArray(SmbGspEnum,this->lid,&gsp,&m);
		this->inputs->GetArray(SmbWEnum,this->lid,&W,&m);
		this->inputs->GetArray(SmbAEnum,this->lid,&a,&m);
		this->inputs->GetArray(SmbTEnum,this->lid,&T,&m);
		EC_input = this->GetInput(SmbECDtEnum);  _assert_(EC_input);
		EC_input->GetInputAverage(&EC);

		//fixed lower temperature bounday condition - T is fixed
		_assert_(m>0);
		T_bottom=T[m-1];
	} /*}}}*/

	// determine initial mass [kg]
	initMass=0; for(int i=0;i<m;i++) initMass += dz[i]*d[i] + W[i];

	// initialize cumulative variables
	sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0; sumMsurf = 0;

	//before starting loop, realize that the transient core runs this smb_core at time = time +deltaT.
   //go back to time - deltaT:
   time-=dt;

	if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << timeinputs/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n");

	/*Get daily accumulated inputs {{{*/
	if (count>1){
		Input *sumEC_input         = this->GetInput(SmbECEnum);  _assert_(sumEC_input);
		Input *sumM_input          = this->GetInput(SmbMeltEnum);  _assert_(sumM_input);
		Input *sumR_input          = this->GetInput(SmbRunoffEnum);  _assert_(sumR_input);
		Input *sumP_input          = this->GetInput(SmbPrecipitationEnum);  _assert_(sumP_input);
		Input *ULW_input           = this->GetInput(SmbMeanULWEnum);  _assert_(ULW_input);
		Input *LW_input            = this->GetInput(SmbNetLWEnum);  _assert_(LW_input);
		Input *SW_input            = this->GetInput(SmbNetSWEnum);  _assert_(SW_input);
		Input *LHF_input           = this->GetInput(SmbMeanLHFEnum);  _assert_(LHF_input);
		Input *SHF_input           = this->GetInput(SmbMeanSHFEnum);  _assert_(SHF_input);
		Input *DzAdd_input         = this->GetInput(SmbDzAddEnum);  _assert_(DzAdd_input);
		Input *MassAdd_input       = this->GetInput(SmbMAddEnum);  _assert_(MassAdd_input);
		Input *InitMass_input      = this->GetInput(SmbMInitnum);  _assert_(InitMass_input);
		Input *sumMsurf_input         = this->GetInput(SmbMSurfEnum);  _assert_(sumMsurf_input);

		ULW_input->GetInputAverage(&meanULW);
		LW_input->GetInputAverage(&netLW);
		SW_input->GetInputAverage(&netSW);
		LHF_input->GetInputAverage(&meanLHF);
		SHF_input->GetInputAverage(&meanSHF);
		DzAdd_input->GetInputAverage(&dz_add);
		MassAdd_input->GetInputAverage(&sumMassAdd);
		sumMassAdd=sumMassAdd*dt;
		InitMass_input->GetInputAverage(&initMass);
		sumEC_input->GetInputAverage(&sumEC);
		sumEC=sumEC*dt*rho_ice;
		sumM_input->GetInputAverage(&sumM);
		sumM=sumM*dt*rho_ice;
		sumMsurf_input->GetInputAverage(&sumMsurf);
      sumMsurf=sumMsurf*dt*rho_ice;
		sumR_input->GetInputAverage(&sumR);
		sumR=sumR*dt*rho_ice;
		sumP_input->GetInputAverage(&sumP);
		sumP=sumP*dt*rho_ice;
	}
	/*}}}*/

	// Get time forcing inputs
	Input *Ta_input  = this->GetInput(SmbTaEnum,timeinputs);    _assert_(Ta_input);
	Input *V_input   = this->GetInput(SmbVEnum,timeinputs);     _assert_(V_input);
	Input *Dlwr_input= this->GetInput(SmbDlwrfEnum,timeinputs); _assert_(Dlwr_input);
	Input *Dswr_input= this->GetInput(SmbDswrfEnum,timeinputs); _assert_(Dswr_input);
	Input *P_input   = this->GetInput(SmbPEnum,timeinputs);     _assert_(P_input);
	Input *eAir_input= this->GetInput(SmbEAirEnum,timeinputs);  _assert_(eAir_input);
	Input *pAir_input= this->GetInput(SmbPAirEnum,timeinputs);  _assert_(pAir_input);
	Input *teValue_input= this->GetInput(SmbTeValueEnum,timeinputs); _assert_(teValue_input);
	Input *aValue_input= this->GetInput(SmbAValueEnum,timeinputs); _assert_(aValue_input);

	/*extract daily data:{{{*/
	Ta_input->GetInputValue(&Ta,gauss);//screen level air temperature [K]
	V_input->GetInputValue(&V,gauss);  //wind speed [m s-1]
	Dlwr_input->GetInputValue(&dlw,gauss);   //downward longwave radiation flux [W m-2]
	Dswr_input->GetInputValue(&dsw,gauss);   //downward shortwave radiation flux [W m-2]
	P_input->GetInputValue(&P,gauss);        //precipitation [kg m-2]
	eAir_input->GetInputValue(&eAir,gauss);  //screen level vapor pressure [Pa]
	pAir_input->GetInputValue(&pAir,gauss);  // screen level air pressure [Pa]
	teValue_input->GetInputValue(&teValue,gauss);  // Emissivity [0-1]
	aValue_input->GetInputValue(&aValue,gauss);  // Albedo [0 1]
	//_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
	/*}}}*/

	/*Snow grain metamorphism:*/
	if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid());

	/*Snow, firn and ice albedo:*/
	if(isalbedo)albedo(&a,aIdx,re,dz,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,Msurf,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());

	/*Distribution of absorbed short wave radation with depth:*/
	if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());

	/*Calculate net shortwave [W m-2]*/
	netSW = netSW + cellsum(swf,m)*smb_dt/dt;

	if(isconstrainsurfaceT){
		if (m>0) T[0]=Ta;
		if (m>1) T[1]=Ta;
	}
	/*Thermal profile computation:*/
	if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid(),isconstrainsurfaceT);

	/*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.
	 * need to fix this in case all or more of cell evaporates */
	dz[0] = dz[0] + EC / d[0];

	/*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
	if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, dsnowIdx, Tmean, Ta, P, dzMin, aSnow, C, V, Vmean, rho_ice,this->Sid());

	/*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
	 * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
	if(ismelt)melt(&M, &Msurf, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop, zY, rho_ice,this->Sid());

	/*Allow non-melt densification and determine compaction [m]*/
	if(isdensification)densification(&d,&dz, T, re, denIdx, aIdx, swIdx, adThresh, C, smb_dt, Tmean,rho_ice,m,this->Sid());

	/*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
	 * sub-time step in thermo equations*/
	//ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here

	/*Calculate net longwave [W m-2]*/
	meanULW = meanULW + ulw*smb_dt/dt;
	netLW = netLW + (dlw - ulw)*smb_dt/dt;

	/*Calculate turbulent heat fluxes [W m-2]*/
	if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());

	/*Verbose some results in debug mode: {{{*/
	if(VerboseSmb() && 0){
		_printf_("smb log: count[" << count << "] m[" << m << "] "
					<< setprecision(16)   << "T[" << cellsum(T,m)  << "] "
					<< "d[" << cellsum(d,m)  << "] "
					<< "dz[" << cellsum(dz,m)  << "] "
					<< "a[" << cellsum(a,m)  << "] "
					<< "W[" << cellsum(W,m)  << "] "
					<< "re[" << cellsum(re,m)  << "] "
					<< "gdn[" << cellsum(gdn,m)  << "] "
					<< "gsp[" << cellsum(gsp,m)  << "] "
					<< "swf[" << netSW << "] "
					<< "lwf[" << netLW << "] "
					<< "a[" << a << "] "
					<< "te[" << teValue << "] "
					<< "\n");
	} /*}}}*/

	meanLHF = meanLHF + lhf*smb_dt/dt;
	meanSHF = meanSHF + shf*smb_dt/dt;

	/*Sum component mass changes [kg m-2]*/
	sumMassAdd = mAdd + sumMassAdd;
	sumM = M + sumM;
	sumMsurf = Msurf + sumMsurf;
	sumR = R + sumR;
	sumW = cellsum(W,m);
	sumP = P +  sumP;
	sumEC = sumEC + EC;  // evap (-)/cond(+)

	/*Calculate total system mass:*/
	sumMass=0;
	fac=0;
	for(int i=0;i<m;i++){
		sumMass += dz[i]*d[i];
		if (d[i] > 0) fac += dz[i]*(rho_ice - fmin(d[i],rho_ice));
	}

	#if defined(_HAVE_AD_)
	/*we want to avoid the round operation at all cost. Not differentiable.*/
	_error_("not implemented yet");
	#else
	dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
	dMass = round(dMass * 100.0)/100.0;

	/*Check mass conservation:*/
	if (dMass != 0.0){
		_printf_("total system mass not conserved in MB function \n");
	}
	#endif

	/*Check bottom grid cell T is unchanged:*/
	if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0){
		if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
	}

	/*Save accumulated inputs {{{*/
   Input *accsumEC_input         = this->GetInput(SmbAccumulatedECEnum);  _assert_(accsumEC_input);
   Input *accsumM_input          = this->GetInput(SmbAccumulatedMeltEnum);  _assert_(accsumM_input);
   Input *accsumR_input          = this->GetInput(SmbAccumulatedRunoffEnum);  _assert_(accsumR_input);
   Input *accsumP_input          = this->GetInput(SmbAccumulatedPrecipitationEnum);  _assert_(accsumP_input);
   Input *accsumSMB_input        = this->GetInput(SmbAccumulatedMassBalanceEnum);  _assert_(accsumP_input);

   accsumEC_input->GetInputAverage(&accsumEC);
   accsumM_input->GetInputAverage(&accsumM);
   accsumR_input->GetInputAverage(&accsumR);
   accsumP_input->GetInputAverage(&accsumP);
   accsumSMB_input->GetInputAverage(&accsumSMB);

	this->SetElementInput(SmbAccumulatedECEnum,accsumEC+EC/rho_ice);
	this->SetElementInput(SmbAccumulatedMassBalanceEnum,accsumSMB+(P + EC -R)/rho_ice);
	this->SetElementInput(SmbAccumulatedMeltEnum,accsumM+M/rho_ice);
	this->SetElementInput(SmbAccumulatedRunoffEnum,accsumR+R/rho_ice);
	this->SetElementInput(SmbAccumulatedPrecipitationEnum,accsumP+P/rho_ice);
   /*}}}*/

	/*Save generated inputs: */
	this->inputs->SetArrayInput(SmbDzEnum,this->lid,dz,m);
	this->inputs->SetArrayInput(SmbDEnum,this->lid,d,m);
	this->inputs->SetArrayInput(SmbReEnum,this->lid,re,m);
	this->inputs->SetArrayInput(SmbGdnEnum,this->lid,gdn,m);
	this->inputs->SetArrayInput(SmbGspEnum,this->lid,gsp,m);
	this->inputs->SetArrayInput(SmbTEnum,this->lid,T,m);
	this->inputs->SetArrayInput(SmbWEnum,this->lid,W,m);
	this->inputs->SetArrayInput(SmbAEnum,this->lid,a,m);
	this->SetElementInput(SmbECEnum,sumEC/dt/rho_ice);
	this->SetElementInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/dt/rho_ice);
	this->SetElementInput(SmbMeltEnum,sumM/dt/rho_ice);
	this->SetElementInput(SmbRunoffEnum,sumR/dt/rho_ice);
	this->SetElementInput(SmbPrecipitationEnum,sumP/dt/rho_ice);
	this->SetElementInput(SmbMeanULWEnum,meanULW);
	this->SetElementInput(SmbNetLWEnum,netLW);
	this->SetElementInput(SmbNetSWEnum,netSW);
	this->SetElementInput(SmbMeanLHFEnum,meanLHF);
	this->SetElementInput(SmbMeanSHFEnum,meanSHF);
	this->SetElementInput(SmbDzAddEnum,dz_add);
	this->SetElementInput(SmbMInitnum,initMass);
	this->SetElementInput(SmbMAddEnum,sumMassAdd/dt);
	this->SetElementInput(SmbMSurfEnum,sumMsurf/dt/rho_ice);
	this->SetElementInput(SmbWAddEnum,sumW/dt);
	this->SetElementInput(SmbFACEnum,fac/1000.); // output in meters
	this->SetElementInput(SmbECDtEnum,EC);

	/*Free allocations:{{{*/
	if(dz) xDelete<IssmDouble>(dz);
	if(d) xDelete<IssmDouble>(d);
	if(re) xDelete<IssmDouble>(re);
	if(gdn) xDelete<IssmDouble>(gdn);
	if(gsp) xDelete<IssmDouble>(gsp);
	if(W) xDelete<IssmDouble>(W);
	if(a) xDelete<IssmDouble>(a);
	if(T) xDelete<IssmDouble>(T);
	if(dzini) xDelete<IssmDouble>(dzini);
	if(dini) xDelete<IssmDouble>(dini);
	if(reini) xDelete<IssmDouble>(reini);
	if(gdnini) xDelete<IssmDouble>(gdnini);
	if(gspini) xDelete<IssmDouble>(gspini);
	if(Wini) xDelete<IssmDouble>(Wini);
	if(aini) xDelete<IssmDouble>(aini);
	if(Tini) xDelete<IssmDouble>(Tini);
	if(swf) xDelete<IssmDouble>(swf);

	delete gauss;
	/*}}}*/
}
/*}}}*/
void       Element::StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble dvx[3];
	IssmDouble dvy[3];

	/*Check that both inputs have been found*/
	if(!vx_input || !vy_input){
		_error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n");
	}

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
	epsilon[0] = dvx[0];	// normal strain rate x-direction
	epsilon[1] = dvy[1]; // normal strain rate y-direction
	epsilon[2] = 0.5*(dvx[1] + dvy[0]); // shear strain rate
	epsilon[3] = 0.5*(dvx[1] - dvy[0]); // rotation rate

}/*}}}*/
void       Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
	/*Compute the 3d Strain Rate (6 components):
	 *
	 * epsilon=[exx eyy ezz exy exz eyz]
	 */

	/*Intermediaries*/
	IssmDouble dvx[3];
	IssmDouble dvy[3];
	IssmDouble dvz[3];

	/*Check that both inputs have been found*/
	if (!vx_input || !vy_input || !vz_input){
		_error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << ", vz: " << vz_input << "\n");
	}

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
	vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);

	epsilon[0] = dvx[0];
	epsilon[1] = dvy[1];
	epsilon[2] = dvz[2];
	epsilon[3] = 0.5*(dvx[1] + dvy[0]);
	epsilon[4] = 0.5*(dvx[2] + dvz[0]);
	epsilon[5] = 0.5*(dvy[2] + dvz[1]);

}/*}}}*/
void       Element::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
	/*Compute the 3d Blatter/HOStrain Rate (5 components):
	 *
	 * epsilon=[exx eyy exy exz eyz]
	 *
	 * with exz=1/2 du/dz
	 *      eyz=1/2 dv/dz
	 *
	 * the contribution of vz is neglected
	 */

	/*Intermediaries*/
	IssmDouble dvx[3];
	IssmDouble dvy[3];

	/*Check that both inputs have been found*/
	if (!vx_input || !vy_input){
		_error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n");
	}

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
	epsilon[0] = dvx[0];
	epsilon[1] = dvy[1];
	epsilon[2] = 0.5*(dvx[1] + dvy[0]);
	epsilon[3] = 0.5*dvx[2];
	epsilon[4] = 0.5*dvy[2];

}/*}}}*/
void       Element::StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
	/*Compute the 2d Blatter/HOStrain Rate (2 components):
	 *
	 * epsilon=[exx exz]
	 *
	 * with exz=1/2 du/dz
	 *
	 * the contribution of vz is neglected
	 */

	/*Intermediaries*/
	IssmDouble dvx[3];

	/*Check that both inputs have been found*/
	if (!vx_input){
		_error_("Input missing. Here are the input pointers we have for vx: " << vx_input <<"\n");
	}

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	epsilon[0] = dvx[0];
	epsilon[1] = 0.5*dvx[1];

}/*}}}*/
void       Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble dvx[3];
	IssmDouble dvy[3];

	/*Check that both inputs have been found*/
	if(!vx_input || !vy_input){
		_error_("Input missing. Here are the input pointers we have for vx: " << vx_input << ", vy: " << vy_input << "\n");
	}

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
	epsilon[0] = dvx[0];
	epsilon[1] = dvy[1];
	epsilon[2] = 0.5*(dvx[1] + dvy[0]);

}/*}}}*/
void       Element::StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble dvx[3];

	/*Check that both inputs have been found*/
	if (!vx_input){
		_error_("Input missing. Here are the input pointers we have for vx: " << vx_input << "\n");
	}

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	*epsilon = dvx[0];

}/*}}}*/
void       Element::StressMaxPrincipalCreateInput(void){/*{{{*/

	/*Intermediaries*/
	IssmDouble *xyz_list = NULL;
	IssmDouble  sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz;
	IssmDouble  a,b,c,d,x[3],max;
	int         dim,numroots;

	/*First: get stress tensor*/
	this->ComputeStressTensor();

	/*Get domain dimension*/
	this->FindParam(&dim,DomainDimensionEnum);

	/*Fetch number vertices and allocate memory*/
	const int NUM_VERTICES  = this->GetNumberOfVertices();

	IssmDouble* maxprincipal = xNew<IssmDouble>(NUM_VERTICES);

	/*Retrieve all inputs and parameters*/
	this->GetVerticesCoordinatesBase(&xyz_list);
	Input* sigma_xx_input  = this->GetInput(StressTensorxxEnum); _assert_(sigma_xx_input);
	Input* sigma_yy_input  = this->GetInput(StressTensoryyEnum); _assert_(sigma_yy_input);
	Input* sigma_xy_input  = this->GetInput(StressTensorxyEnum); _assert_(sigma_xy_input);
	Input* sigma_xz_input  = NULL;
	Input* sigma_yz_input  = NULL;
	Input* sigma_zz_input  = NULL;
	if(dim==3){
		sigma_xz_input  = this->GetInput(StressTensorxzEnum); _assert_(sigma_xz_input);
		sigma_yz_input  = this->GetInput(StressTensoryzEnum); _assert_(sigma_yz_input);
		sigma_zz_input  = this->GetInput(StressTensorzzEnum); _assert_(sigma_zz_input);
	}

	/*loop over vertices: */
	Gauss* gauss=this->NewGauss();
	for (int iv=0;iv<NUM_VERTICES;iv++){
		gauss->GaussVertex(iv);

		sigma_xx_input->GetInputValue(&sigma_xx,gauss);
		sigma_yy_input->GetInputValue(&sigma_yy,gauss);
		sigma_xy_input->GetInputValue(&sigma_xy,gauss);
		if(dim==3){
			sigma_xz_input->GetInputValue(&sigma_xz,gauss);
			sigma_yz_input->GetInputValue(&sigma_yz,gauss);
			sigma_zz_input->GetInputValue(&sigma_zz,gauss);
		}

		if(dim==2){
			a = 0.;
			b = 1.;
			c = -sigma_yy -sigma_xx;
			d = sigma_xx*sigma_yy - sigma_xy*sigma_xy;
		}
		else{
			a = -1.;
			b = sigma_xx+sigma_yy+sigma_zz;
			c = -sigma_xx*sigma_yy -sigma_xx*sigma_zz -sigma_yy*sigma_zz + sigma_xy*sigma_xy +sigma_xz*sigma_xz +sigma_yz*sigma_yz;
			d = sigma_xx*sigma_yy*sigma_zz - sigma_xx*sigma_yz*sigma_yz -sigma_yy*sigma_xz*sigma_xz - sigma_zz*sigma_xy*sigma_xy + 2.*sigma_xy*sigma_xz*sigma_yz;
		}

		/*Get roots of polynomials*/
		cubic(a,b,c,d,x,&numroots);

		/*Initialize maximum eigne value*/
		if(numroots>0){
			max = fabs(x[0]);
		}
		else{
			_error_("No eigen value found");
		}

		/*Get max*/
		for(int i=1;i<numroots;i++){
			if(fabs(x[i])>max) max = fabs(x[i]);
		}

		maxprincipal[iv]=max;
	}

	/*Create input*/
	this->AddInput(StressMaxPrincipalEnum,maxprincipal,P1Enum);

	/*Clean up and return*/
	xDelete<IssmDouble>(maxprincipal);
	xDelete<IssmDouble>(xyz_list);
	delete gauss;
}
/*}}}*/
IssmDouble Element::TotalFloatingBmb(IssmDouble* mask, bool scaled){/*{{{*/

	/*Retrieve values of the mask defining the element: */
	for(int i=0;i<this->GetNumberOfVertices();i++){
		if(mask[this->vertices[i]->Sid()]<=0.){
			return 0.;
		}
	}

	/*Return: */
	return this->TotalFloatingBmb(scaled);
}
/*}}}*/
IssmDouble Element::TotalGroundedBmb(IssmDouble* mask, bool scaled){/*{{{*/

	/*Retrieve values of the mask defining the element: */
	for(int i=0;i<this->GetNumberOfVertices();i++){
		if(mask[this->vertices[i]->Sid()]<=0.){
			return 0.;
		}
	}

	/*Return: */
	return this->TotalGroundedBmb(scaled);
}
/*}}}*/
IssmDouble Element::TotalSmb(IssmDouble* mask, bool scaled){/*{{{*/

	/*Retrieve values of the mask defining the element: */
	for(int i=0;i<this->GetNumberOfVertices();i++){
		if(mask[this->vertices[i]->Sid()]<=0.){
			return 0.;
		}
	}

	/*Return: */
	return this->TotalSmb(scaled);
}
/*}}}*/
void       Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/

	/*All nodes have the same Coordinate System*/
	int numnodes  = this->GetNumberOfNodes();
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	TransformInvStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void       Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/

	int         i,j;
	int         numdofs   = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current stiffness matrix*/
	values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
	for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);

	/*Transform matrix: R*Ke*R^T */
	TripleMultiply(transform,numdofs,numdofs,0,
				values,Ke->nrows,Ke->ncols,0,
				transform,numdofs,numdofs,1,
				&Ke->values[0],0);

	/*Free Matrix*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}/*}}}*/
void       Element::TransformLoadVectorCoord(ElementVector* pe,int transformenum){/*{{{*/

	/*All nodes have the same Coordinate System*/
	int  numnodes = this->GetNumberOfNodes();
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	this->TransformLoadVectorCoord(pe,this->nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void       Element::TransformLoadVectorCoord(ElementVector* pe,int* cs_array){/*{{{*/

	this->TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),cs_array);

}/*}}}*/
void       Element::TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/

	int         i;
	int         numdofs   = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current load vector*/
	values=xNew<IssmDouble>(pe->nrows);
	for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);

	/*Transform matrix: R^T*pe */
	MatrixMultiply(transform,numdofs,numdofs,1,
				values,pe->nrows,1,0,
				&pe->values[0],0);

	/*Free Matrices*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}/*}}}*/
void       Element::TransformSolutionCoord(IssmDouble* values,int transformenum){/*{{{*/

	/*All nodes have the same Coordinate System*/
	int  numnodes = this->GetNumberOfNodes();
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	this->TransformSolutionCoord(values,this->nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void       Element::TransformSolutionCoord(IssmDouble* values,int* transformenum_list){/*{{{*/
	this->TransformSolutionCoord(values,this->nodes,this->GetNumberOfNodes(),transformenum_list);
}/*}}}*/
void       Element::TransformSolutionCoord(IssmDouble* values,int numnodes,int transformenum){/*{{{*/

	/*All nodes have the same Coordinate System*/
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	this->TransformSolutionCoord(values,this->nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void       Element::TransformSolutionCoord(IssmDouble* solution,int numnodes,int* cs_array){/*{{{*/
	this->TransformSolutionCoord(solution,this->nodes,numnodes,cs_array);
}/*}}}*/
void       Element::TransformSolutionCoord(IssmDouble* values,Node** nodes_list,int numnodes,int transformenum){/*{{{*/
	/*NOT NEEDED*/
	/*All nodes have the same Coordinate System*/
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	this->TransformSolutionCoord(values,nodes_list,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void       Element::TransformSolutionCoord(IssmDouble* solution,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/

	int         i;
	int         numdofs   = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current solution vector*/
	values=xNew<IssmDouble>(numdofs);
	for(i=0;i<numdofs;i++) values[i]=solution[i];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);

	/*Transform matrix: R*U */
	MatrixMultiply(transform,numdofs,numdofs,0,
				values,numdofs,1,0,
				&solution[0],0);

	/*Free Matrices*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}/*}}}*/
void       Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/

	/*All nodes have the same Coordinate System*/
	int  numnodes = this->GetNumberOfNodes();
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	this->TransformStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void       Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){/*{{{*/
	this->TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum_list);
}/*}}}*/
void       Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/

	int         numdofs = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(int i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current stiffness matrix*/
	values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
	for(int i=0;i<Ke->nrows*Ke->ncols;i++) values[i]=Ke->values[i];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);

	/*Transform matrix: R^T*Ke*R */
	TripleMultiply(transform,numdofs,numdofs,1,
				values,Ke->nrows,Ke->ncols,0,
				transform,numdofs,numdofs,0,
				&Ke->values[0],0);

	/*Free Matrix*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}/*}}}*/
void       Element::ViscousHeatingCreateInput(void){/*{{{*/

	/*Intermediaries*/
	IssmDouble phi;
	IssmDouble thickness;
	IssmDouble *xyz_list = NULL;

	/*Fetch number vertices and allocate memory*/
	const int NUM_VERTICES = this->GetNumberOfVertices();

	IssmDouble* viscousheating = xNew<IssmDouble>(NUM_VERTICES);

	/*Retrieve all inputs and parameters*/
	this->GetVerticesCoordinates(&xyz_list);
	Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input);
	Input* vz_input = this->GetInput(VzEnum); _assert_(vz_input);

	/*loop over vertices: */
	Gauss* gauss=this->NewGauss();
	for (int iv=0;iv<NUM_VERTICES;iv++){
		gauss->GaussVertex(iv);

		this->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);

		viscousheating[iv]=phi;
	}

	/*Create PentaVertex input, which will hold the basal friction:*/
	this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum);

	/*Clean up and return*/
	xDelete<IssmDouble>(viscousheating);
	xDelete<IssmDouble>(xyz_list);
	delete gauss;
}
/*}}}*/

/*Enthalpy*/
void       Element::ThermalToEnthalpy(IssmDouble * penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/

	/*Ouput*/
	IssmDouble enthalpy;

	/*Get necessary parameters*/
	IssmDouble latentheat,referencetemperature,heatcapacity;
	parameters->FindParam(&latentheat,MaterialsLatentheatEnum);
	parameters->FindParam(&referencetemperature,ConstantsReferencetemperatureEnum);
	parameters->FindParam(&heatcapacity,MaterialsHeatcapacityEnum);

	if(temperature<TMeltingPoint(pressure)){
		enthalpy=heatcapacity*(temperature-referencetemperature);
	}
	else{
		enthalpy=PureIceEnthalpy(pressure)+latentheat*waterfraction;
	}

	/*Assign output pointers:*/
	*penthalpy=enthalpy;
}
/*}}}*/
IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/

	/*Get necessary parameters*/
	IssmDouble beta,meltingpoint;
	parameters->FindParam(&beta,MaterialsBetaEnum);
	parameters->FindParam(&meltingpoint,MaterialsMeltingpointEnum);

	return meltingpoint-beta*pressure;
}
/*}}}*/
void       Element::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/

	/*Ouput*/
	IssmDouble temperature,waterfraction;

	/*Get necessary parameters*/
	IssmDouble latentheat,referencetemperature,heatcapacity;
	parameters->FindParam(&latentheat,MaterialsLatentheatEnum);
	parameters->FindParam(&referencetemperature,ConstantsReferencetemperatureEnum);
	parameters->FindParam(&heatcapacity,MaterialsHeatcapacityEnum);

	if(enthalpy<PureIceEnthalpy(pressure)){
		temperature=referencetemperature+enthalpy/heatcapacity;
		waterfraction=0.;
	}
	else{
		temperature=TMeltingPoint(pressure);
		waterfraction=(enthalpy-PureIceEnthalpy(pressure))/latentheat;
	}

	/*Assign output pointers:*/
	*pwaterfraction=waterfraction;
	*ptemperature=temperature;
}
/*}}}*/
IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/

	/*Get necessary parameters*/
	IssmDouble heatcapacity,thermalconductivity,temperateiceconductivity;
	parameters->FindParam(&heatcapacity,MaterialsHeatcapacityEnum);
	parameters->FindParam(&thermalconductivity,MaterialsThermalconductivityEnum);
	parameters->FindParam(&temperateiceconductivity,MaterialsTemperateiceconductivityEnum);

	if(enthalpy<PureIceEnthalpy(pressure)){
		return thermalconductivity/heatcapacity;
	}
	else{
		return temperateiceconductivity/heatcapacity;
	}
}
/*}}}*/
IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/

	IssmDouble  lambda;                 // fraction of cold ice
	IssmDouble  kappa,kappa_c,kappa_t;  //enthalpy conductivities
	IssmDouble  Hc,Ht;
	IssmDouble* PIE   = xNew<IssmDouble>(numvertices);
	IssmDouble* dHpmp = xNew<IssmDouble>(numvertices);

	for(int iv=0; iv<numvertices; iv++){
		PIE[iv]=PureIceEnthalpy(pressure[iv]);
		dHpmp[iv]=enthalpy[iv]-PIE[iv];
	}

	bool allequalsign=true;
	if(dHpmp[0]<0){
		for(int iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]<0));
	}
	else{
		for(int iv=1; iv<numvertices;iv++) allequalsign=(allequalsign && (dHpmp[iv]>=0));
	}

	if(allequalsign){
		kappa=EnthalpyDiffusionParameter(enthalpy[0], pressure[0]);
	}
	else {
		/* return harmonic mean of thermal conductivities, weighted by fraction of cold/temperate ice,
			cf Patankar 1980, pp44 */
		kappa_c=EnthalpyDiffusionParameter(PureIceEnthalpy(0.)-1.,0.);
		kappa_t=EnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.);
		Hc=0.; Ht=0.;
		for(int iv=0; iv<numvertices;iv++){
			if(enthalpy[iv]<PIE[iv])
			 Hc+=(PIE[iv]-enthalpy[iv]);
			else
			 Ht+=(enthalpy[iv]-PIE[iv]);
		}
		_assert_((Hc+Ht)>0.);
		lambda = Hc/(Hc+Ht);
		kappa  = 1./(lambda/kappa_c + (1.-lambda)/kappa_t);
	}

	/*Clean up and return*/
	xDelete<IssmDouble>(PIE);
	xDelete<IssmDouble>(dHpmp);
	return kappa;
}
/*}}}*/
IssmDouble Element::PureIceEnthalpy(IssmDouble pressure){/*{{{*/

	/*Get necessary parameters*/
	IssmDouble referencetemperature,heatcapacity;
	parameters->FindParam(&referencetemperature,ConstantsReferencetemperatureEnum);
	parameters->FindParam(&heatcapacity,MaterialsHeatcapacityEnum);

	return heatcapacity*(TMeltingPoint(pressure)-referencetemperature);
}
/*}}}*/
