Element.cpp

Go to the documentation of this file.

 
/*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 "../Inputs2/BoolInput2.h"
#include "../Inputs2/TransientInput2.h"
#include "../Inputs2/ElementInput2.h"
#include "../Inputs2/PentaInput2.h"
#include "../Inputs2/DatasetInput2.h"
#include "../Inputs2/ArrayInput2.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->inputs2    = 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;
   int         dim;
   IssmDouble *xyz_list = NULL;
 
   /*Retrieve all inputs we will be needing: */
   this->GetVerticesCoordinates(&xyz_list);
   parameters->FindParam(&dim,DomainDimensionEnum);
   Input2* vx_input=this->GetInput2(VxEnum); _assert_(vx_input);
   Input2* vy_input=this->GetInput2(VyEnum); _assert_(vy_input);
   Input2* vz_input=NULL;
   if(dim==3){vz_input=this->GetInput2(VzEnum); _assert_(vz_input);}
 
   /*Allocate arrays*/
   const int NUM_VERTICES = this->GetNumberOfVertices();
 
   IssmDouble* lambdas = xNew<IssmDouble>(NUM_VERTICES);
 
   /* 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->AddInput2(LambdaSEnum,lambdas,P1Enum);
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(lambdas);
 
}
/*}}}*/
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 */
   Input2* eps_xx_input=this->GetInput2(StrainRatexxEnum); _assert_(eps_xx_input);
   Input2* eps_yy_input=this->GetInput2(StrainRateyyEnum); _assert_(eps_yy_input);
   Input2* eps_xy_input=this->GetInput2(StrainRatexyEnum); _assert_(eps_xy_input);
   Input2* eps_xz_input=NULL;
   Input2* eps_yz_input=NULL;
   Input2* eps_zz_input=NULL;
   if(dim==3){
      eps_xz_input=this->GetInput2(StrainRatexzEnum); _assert_(eps_xz_input);
      eps_yz_input=this->GetInput2(StrainRateyzEnum); _assert_(eps_yz_input);
      eps_zz_input=this->GetInput2(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 */
   Input2* n_input=this->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
   Input2* damage_input = NULL;
   Input2* B_input = NULL;
   int domaintype;
   parameters->FindParam(&domaintype,DomainTypeEnum);
   if(domaintype==Domain2DhorizontalEnum){
      damage_input = this->GetInput2(DamageDbarOldEnum);  _assert_(damage_input);
      B_input=this->GetInput2(MaterialsRheologyBbarEnum); _assert_(B_input);
   }
   else{
      damage_input = this->GetInput2(DamageDOldEnum);   _assert_(damage_input);
      B_input=this->GetInput2(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->AddInput2(NewDamageEnum,newD,P1DGEnum);
   if(domaintype==Domain2DhorizontalEnum){
      this->AddInput2(DamageDbarEnum,newD,this->GetElementType());
   }
   else{
      this->AddInput2(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);
   Input2* vx_input=this->GetInput2(VxEnum); _assert_(vx_input);
   Input2* vy_input=this->GetInput2(VyEnum); _assert_(vy_input);
   Input2* vz_input=NULL;
   if(dim==3){vz_input=this->GetInput2(VzEnum); _assert_(vz_input);}
 
   /*Allocate arrays*/
   const int NUM_VERTICES = this->GetNumberOfVertices();
 
   IssmDouble* eps_xx = xNew<IssmDouble>(NUM_VERTICES);
   IssmDouble* eps_yy = xNew<IssmDouble>(NUM_VERTICES);
   IssmDouble* eps_zz = xNew<IssmDouble>(NUM_VERTICES);
   IssmDouble* eps_xy = xNew<IssmDouble>(NUM_VERTICES);
   IssmDouble* eps_xz = xNew<IssmDouble>(NUM_VERTICES);
   IssmDouble* eps_yz = xNew<IssmDouble>(NUM_VERTICES);
   IssmDouble* eps_ef = xNew<IssmDouble>(NUM_VERTICES);
 
   /* 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->AddInput2(StrainRatexxEnum,eps_xx,P1Enum);
   this->AddInput2(StrainRatexyEnum,eps_xy,P1Enum);
   this->AddInput2(StrainRatexzEnum,eps_xz,P1Enum);
   this->AddInput2(StrainRateyyEnum,eps_yy,P1Enum);
   this->AddInput2(StrainRateyzEnum,eps_yz,P1Enum);
   this->AddInput2(StrainRatezzEnum,eps_zz,P1Enum);
   this->AddInput2(StrainRateeffectiveEnum,eps_ef,P1Enum);
 
   /*Clean up and return*/
   delete gauss;
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(eps_xx);
   xDelete<IssmDouble>(eps_yy);
   xDelete<IssmDouble>(eps_zz);
   xDelete<IssmDouble>(eps_xy);
   xDelete<IssmDouble>(eps_xz);
   xDelete<IssmDouble>(eps_yz);
   xDelete<IssmDouble>(eps_ef);
 
}
/*}}}*/
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(inputs2) inputs2->DeepEcho();
   else _printf_("inputs2=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*/
   DatasetInput2* dinput1=this->GetDatasetInput2(SmbTemperaturesPresentdayEnum);   _assert_(dinput1);
   DatasetInput2* dinput2=this->GetDatasetInput2(SmbTemperaturesLgmEnum);          _assert_(dinput2);
   DatasetInput2* dinput3=this->GetDatasetInput2(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->inputs2->SetTriaDatasetInput(SmbMonthlytemperaturesEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
         case PentaEnum: this->inputs2->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->inputs2->SetTriaDatasetInput(SmbPrecipitationEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
         case PentaEnum: this->inputs2->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*/
   DatasetInput2 *dinput3 = NULL;
   DatasetInput2 *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->GetDatasetInput2(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->GetDatasetInput2(SmbPrecipitationsReconstructedEnum); _assert_(dinput4);
   }
 
   /*Get present day temp and precip (monthly)*/
   DatasetInput2 *dinput1 = this->GetDatasetInput2(SmbTemperaturesPresentdayEnum);   _assert_(dinput1);
   DatasetInput2 *dinput2 = this->GetDatasetInput2(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->inputs2->SetTriaDatasetInput(SmbMonthlytemperaturesEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
         case PentaEnum: this->inputs2->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->inputs2->SetTriaDatasetInput(SmbPrecipitationEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
         case PentaEnum: this->inputs2->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(&accugrad,SmbAccugradEnum);
   parameters->FindParam(&runoffalti,SmbRunoffaltiEnum);
   parameters->FindParam(&runoffgrad,SmbRunoffgradEnum);
 
   /*Recover reference values at current 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->AddInput2(SmbMassBalanceSubstepEnum,&smb[0],P1Enum);
      this->AddInput2(SmbRunoffSubstepEnum,&runoff[0],P1Enum);
      break;
   case PentaEnum:
      this->AddInput2(SmbMassBalanceSubstepEnum,&smb[0],P1Enum);
      this->AddInput2(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);
   Input2* vx_input = this->GetInput2(VxEnum); _assert_(vx_input);
   Input2* vy_input = this->GetInput2(VyEnum); _assert_(vy_input);
   Input2* vz_input = NULL;
   if(dim==3){
      vz_input = this->GetInput2(VzEnum); _assert_(vz_input);
   }
   this->GetVerticesCoordinates(&xyz_list);
 
   Gauss* gauss=this->NewGauss(5);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
      gauss->GaussPoint(ig);
      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,Input2* vx_input,Input2* vy_input,Input2* 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,Input2* vx_input,Input2* 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,Input2* vx_input,Input2* 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,Input2* vx_input,Input2* 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 (inputs2) inputs2->Echo();
   else _printf_("inputs2=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){/*{{{*/
   Input2 *input    = this->GetInput2(enumtype);
   this->GetInputListOnNodes(pvalue,input,defaultvalue);
}
/*}}}*/
void       Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype){/*{{{*/
 
   Input2 *input    = this->GetInput2(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();
   Input2 *input    = this->GetInput2(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*/
   Input2* input2=this->GetInput2(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*/
   Input2* input=this->GetInput2(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){/*{{{*/
   Input2* input=this->GetInput2(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->inputs2->GetInputValue(pvalue,inputenum,this->lid);
 
}/*}}}*/
void       Element::GetInputValue(int* pvalue,int inputenum){/*{{{*/
   this->inputs2->GetInputValue(pvalue,inputenum,this->lid);
}/*}}}*/
void       Element::GetInput2Value(bool* pvalue,int inputenum){/*{{{*/
 
   this->inputs2->GetInputValue(pvalue,inputenum,this->lid);
 
}/*}}}*/
void       Element::GetInput2Value(int* pvalue,int inputenum){/*{{{*/
 
   this->inputs2->GetInputValue(pvalue,inputenum,this->lid);
 
}/*}}}*/
void       Element::GetInput2Value(IssmDouble* pvalue,int inputenum){/*{{{*/
 
   this->inputs2->GetInputValue(pvalue,inputenum,this->lid);
 
}/*}}}*/
void       Element::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum){/*{{{*/
 
   Input2* input=this->GetInput2(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*/
   Input2* enum_input=this->GetInput2(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;
   Input2*     input       = NULL;
 
   switch(type){
      case ElementSIdEnum:
         input=this->GetInput2(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;
   parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
   parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
 
   /*Get number of vertices*/
   const int NUM_VERTICES = this->GetNumberOfVertices();
   if(isautodiff){
      int* N=NULL;
      int* M=NULL;
      int start = 0;
      parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
      parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
      if(control_index>0) {
         for(int n=0;n<control_index;n++){
            start+=N[n]*M[n];
         }
      }
 
      for(int n=0;n<N[control_index];n++){
         for(int i=0;i<NUM_VERTICES;i++){
            indexing[i+n*NUM_VERTICES]=this->vertices[i]->Sid() + start + n*M[control_index];
         }
      }
   }
   else{
      int M;
      parameters->FindParam(&M,ControlInputSizeMEnum);
      /*get gradient indices*/
      for(int i=0;i<NUM_VERTICES;i++){
         indexing[i]=this->vertices[i]->Sid() + (control_index)*M;
      }
   }
}
/*}}}*/
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(){/*{{{*/
   Input2* input=this->GetInput2(MeshVertexonbaseEnum); _assert_(input);
   return (input->GetInputMax()>0.);
}/*}}}*/
bool       Element::HasNodeOnSurface(){/*{{{*/
   Input2* input=this->GetInput2(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,Inputs2* inputs2,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(inputs2,vector_enum,vector[0]);
      }
      else if(M==iomodel->numberofvertices){
         for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
         this->SetElementInput(inputs2,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];
         inputs2->SetTransientInput(vector_enum,times,N);
         TransientInput2* transientinput = inputs2->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(inputs2,NUM_VERTICES,vertexlids,values,vector_enum);
         }
         else if(N==this->GetNumberOfNodes(P0Enum)){
            this->SetElementInput(inputs2,vector_enum,values[0]);
         }
         else if(N==this->GetNumberOfNodes(P1xP2Enum)){ _assert_(this->ObjectEnum()==PentaEnum);
            inputs2->SetPentaInput(vector_enum,P1xP2Enum,this->lid,N,values);
         }
         else if(N==this->GetNumberOfNodes(P1xP3Enum)){ _assert_(this->ObjectEnum()==PentaEnum);
            inputs2->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(inputs2,vector_enum,reCast<bool>(vector[this->Sid()]));
         }
         else if (code==6){ //integer
            this->SetIntInput(inputs2,vector_enum,reCast<int>(vector[this->Sid()]));
         }
         else if (code==7){ //IssmDouble
            this->SetElementInput(inputs2,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];
         inputs2->SetTransientInput(vector_enum,times,N);
         TransientInput2* transientinput = inputs2->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];
         inputs2->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,Inputs2* inputs2,IoModel* iomodel,int M,int N,IssmDouble scale,int input_enum,int id){/*{{{*/
 
   /*Intermediaries*/
   const int numvertices = this->GetNumberOfVertices();
 
   int        *vertexids   = xNew<int>(numvertices);
   IssmDouble *values      = xNew<IssmDouble>(numvertices);
   IssmDouble *values_min  = xNew<IssmDouble>(numvertices);
   IssmDouble *values_max  = xNew<IssmDouble>(numvertices);
 
   /*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
   }
 
   /*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,inputs2,iomodel,values,values_min,values_max,P1Enum,id);
   }
 
   else if(M==iomodel->numberofvertices+1 && N>1){
      _error_("not supported tet");
      //IssmDouble* times = xNew<IssmDouble>(N);
      //for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
      //TransientInput* values_input=new TransientInput(input_enum,times,N);
      //TransientInput* mins_input = new TransientInput(ControlInputMinsEnum,times,N);
      //TransientInput* maxs_input = new TransientInput(ControlInputMaxsEnum,times,N);
      //TransientInput* grad_input = new TransientInput(ControlInputGradEnum);
      //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] = min_vector[N*(vertexids[i]-1)+t];
      //    values_max[i] = max_vector[N*(vertexids[i]-1)+t];
      // }
      // switch(this->ObjectEnum()){
      //    case TriaEnum:
      //       values_input->AddTimeInput(new TriaInput(input_enum,values,P1Enum));
      //       mins_input->AddTimeInput(new TriaInput(ControlInputMinsEnum,values_min,P1Enum));
      //       maxs_input->AddTimeInput(new TriaInput(ControlInputMaxsEnum,values_max,P1Enum));
      //       break;
      //    case PentaEnum:
      //       values_input->AddTimeInput(new PentaInput(input_enum,values,P1Enum));
      //       mins_input->AddTimeInput(new PentaInput(ControlInputMinsEnum,values_min,P1Enum));
      //       maxs_input->AddTimeInput(new PentaInput(ControlInputMaxsEnum,values_max,P1Enum));
      //       break;
      //    case TetraEnum:
      //       values_input->AddTimeInput(new TetraInput(input_enum,values,P1Enum));
      //       mins_input->AddTimeInput(new TetraInput(ControlInputMinsEnum,values_min,P1Enum));
      //       maxs_input->AddTimeInput(new TetraInput(ControlInputMaxsEnum,values_max,P1Enum));
      //       break;
      //    default: _error_("Not implemented yet");
      // }
      //}
      //this->inputs->AddInput(new ControlInput(input_enum,TransientInputEnum,values_input,mins_input,maxs_input,grad_input,P1Enum,id));
      //xDelete<IssmDouble>(times);
   }
   else _error_("not currently supported type of M and N attempted");
 
   /*clean up*/
   xDelete<IssmDouble>(values);
   xDelete<IssmDouble>(values_min);
   xDelete<IssmDouble>(values_max);
   xDelete<int>(vertexids);
}
/*}}}*/
void       Element::DatasetInputAdd(int enum_type,IssmDouble* vector,Inputs2* inputs2,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->AddInput2(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:  inputs2->SetTriaDatasetInput(enum_type,input_id,P1Enum,NUM_VERTICES,vertexlids,values); break;
            case PentaEnum: inputs2->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];
         TransientInput2* transientinput = inputs2->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");
         //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");
 
      //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->inputs2,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->inputs2,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);
 
   Input2* input = this->GetInput2(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(){/*{{{*/
 
   Input2* input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(input);
   if(input->GetInputMax() > 0.){
      return true;
   }
   else{
      return false;
   }
}/*}}}*/
bool       Element::IsIceInElement(){/*{{{*/
   Input2* input=this->GetInput2(MaskIceLevelsetEnum); _assert_(input);
   return (input->GetInputMin()<0.);
}
/*}}}*/
bool       Element::IsIceOnlyInElement(){/*{{{*/
   Input2* input=this->GetInput2(MaskIceLevelsetEnum); _assert_(input);
   return (input->GetInputMax()<0.);
}
/*}}}*/
bool       Element::IsLandInElement(){/*{{{*/
   Input2* input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(input);
   return (input->GetInputMax()>0.);
}
/*}}}*/
bool       Element::IsOceanInElement(){/*{{{*/
   Input2* input=this->GetInput2(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);
   }
   Input2* tf_input = this->GetInput2(BasalforcingsIsmip6TfShelfEnum);              _assert_(tf_input);
   Input2* meltanomaly_input = this->GetInput2(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->AddInput2(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->AddInput2(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->AddInput2(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->AddInput2(BasalforcingsGeothermalfluxEnum,values,P1Enum);
   xDelete<IssmDouble>(xyz_list);
   xDelete<IssmDouble>(values);
 
}/*}}}*/
void       Element::MarshallElement(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction,int numanalyses){/*{{{*/
 
   _assert_(this);
   if(marshall_direction==MARSHALLING_BACKWARD){
      nodes = NULL;
   }
 
   MARSHALLING_ENUM(ElementEnum);
 
   MARSHALLING(id);
   MARSHALLING(sid);
   MARSHALLING(lid);
   MARSHALLING(element_type);
   MARSHALLING_DYNAMIC(element_type_list,int,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->AddInput2(MaskOceanLevelsetEnum,&phi[0],P1Enum);
 
   /*Update inputs*/
   this->AddInput2(SurfaceEnum,&s[0],P1Enum);
   this->AddInput2(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(){/*{{{*/
   const int NUM_VERTICES = this->GetNumberOfVertices();
 
   IssmDouble  meltratefactor,thresholdthickness,upperdepthmelt;
   IssmDouble* base     = xNew<IssmDouble>(NUM_VERTICES);
   IssmDouble* bed      = xNew<IssmDouble>(NUM_VERTICES);
   IssmDouble* values   = xNew<IssmDouble>(NUM_VERTICES);
 
   parameters->FindParam(&meltratefactor,BasalforcingsMeltrateFactorEnum);
   parameters->FindParam(&thresholdthickness,BasalforcingsThresholdThicknessEnum);
   parameters->FindParam(&upperdepthmelt,BasalforcingsUpperdepthMeltEnum);
 
   this->GetInputListOnVertices(base,BaseEnum);
   this->GetInputListOnVertices(bed,BedEnum);
   for(int i=0;i<NUM_VERTICES;i++){
      if(base[i]>upperdepthmelt){
         values[i]=0;
      }
      else{
         values[i]=meltratefactor*tanh((base[i]-bed[i])/thresholdthickness)*(upperdepthmelt-base[i]);
      }
   }
 
   this->AddInput2(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
   xDelete<IssmDouble>(base);
   xDelete<IssmDouble>(bed);
   xDelete<IssmDouble>(values);
 
}/*}}}*/
void       Element::BeckmannGoosseFloatingiceMeltingRate(){/*{{{*/
 
   int numvertices      = this->GetNumberOfVertices();
   IssmDouble meltratefactor,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    = xNew<IssmDouble>(numvertices);
   IssmDouble* values  = xNew<IssmDouble>(numvertices);
   IssmDouble* oceansalinity   = xNew<IssmDouble>(numvertices);
   IssmDouble* oceantemp       = xNew<IssmDouble>(numvertices);
 
   this->GetInputListOnVertices(base,BaseEnum);
   this->GetInputListOnVertices(oceansalinity,BasalforcingsOceanSalinityEnum);
   this->GetInputListOnVertices(oceantemp,BasalforcingsOceanTempEnum);
   parameters->FindParam(&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 * 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->AddInput2(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
   xDelete<IssmDouble>(base);
   xDelete<IssmDouble>(values);
   xDelete<IssmDouble>(oceantemp);
   xDelete<IssmDouble>(oceansalinity);
   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*/
   DatasetInput2* dinput1=this->GetDatasetInput2(SmbTemperaturesPresentdayEnum);   _assert_(dinput1);
   DatasetInput2* dinput2=this->GetDatasetInput2(SmbTemperaturesLgmEnum);          _assert_(dinput2);
   DatasetInput2* dinput3=this->GetDatasetInput2(SmbPrecipitationsPresentdayEnum); _assert_(dinput3);
   DatasetInput2* dinput4=this->GetDatasetInput2(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->inputs2->SetTriaDatasetInput(SmbMonthlytemperaturesEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
         case PentaEnum: this->inputs2->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->inputs2->SetTriaDatasetInput(SmbPrecipitationEnum,imonth,P1Enum,NUM_VERTICES,vertexlids,tmp); break;
         case PentaEnum: this->inputs2->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.;
 
   Input2* dist_gl_input=this->GetInput2(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
   Input2* dist_cf_input=this->GetInput2(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->inputs2,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);
   Input2 *thickness_input    = this->GetInput2(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);
      Input2 *overturningC_input = this->GetInput2(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->AddInput2(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1DGEnum);
      this->AddInput2(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1DGEnum);
      this->AddInput2(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1DGEnum);
      this->AddInput2(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->AddInput2(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1DGEnum);
      this->AddInput2(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1DGEnum);
      this->AddInput2(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*/
   Input2* zgl_input         = this->GetInput2(GroundinglineHeightEnum);                     _assert_(zgl_input);
   Input2* toc_input         = this->GetInput2(BasalforcingsPicoSubShelfOceanTempEnum);      _assert_(toc_input);
   Input2* soc_input         = this->GetInput2(BasalforcingsPicoSubShelfOceanSalinityEnum);  _assert_(soc_input);
   Input2* base_input        = this->GetInput2(BaseEnum);                                    _assert_(base_input);
   Input2* baseslopex_input  = this->GetInput2(BaseSlopeXEnum);                              _assert_(baseslopex_input);
   Input2* baseslopey_input  = this->GetInput2(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->AddInput2(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*/
   DatasetInput2* dinput =this->GetDatasetInput2(SmbMonthlytemperaturesEnum); _assert_(dinput);
   DatasetInput2* dinput2=this->GetDatasetInput2(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 */
   Input2* input  = NULL;
   Input2* input2 = NULL;
   if(issetpddfac==1){
      input  = this->GetInput2(SmbPddfacSnowEnum); _assert_(input);
      input2 = this->GetInput2(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->AddInput2(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum);
         this->AddInput2(SmbMassBalanceEnum,&agd[0],P1Enum);
         this->AddInput2(SmbAccumulationEnum,&accu[0],P1Enum);
         this->AddInput2(SmbMeltEnum,&melt[0],P1Enum);
         break;
      case PentaEnum:
         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*/
            PentaInput2* temp_input = xDynamicCast<PentaInput2*>(this->GetInput2(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*/
               PentaInput2* enth_input = xDynamicCast<PentaInput2*>(this->GetInput2(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->AddInput2(SmbMassBalanceEnum,&agd[0],P1Enum);
         this->AddInput2(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum);
         this->InputExtrude(TemperaturePDDEnum,-1);
         this->InputExtrude(SmbMassBalanceEnum,-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*/
   DatasetInput2* dinput =this->GetDatasetInput2(SmbMonthlytemperaturesEnum); _assert_(dinput);
   DatasetInput2* dinput2=this->GetDatasetInput2(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->AddInput2(TemperatureEnum,&yearlytemperatures[0],P1Enum);
         this->AddInput2(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum);
         this->AddInput2(SmbMassBalanceEnum,&smb[0],P1Enum);
         this->AddInput2(SmbAccumulationEnum,&accu[0],P1Enum);
         this->AddInput2(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*/
               PentaInput2* temp_input = xDynamicCast<PentaInput2*>(this->GetInput2(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*/
                  PentaInput2* enth_input = xDynamicCast<PentaInput2*>(this->GetInput2(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->AddInput2(SmbMassBalanceEnum,&smb[0],P1Enum);
         this->AddInput2(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum);
         this->AddInput2(SmbAccumulationEnum,&accu[0],P1Enum);
         this->AddInput2(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->inputs2->SetInput(PartitioningEnum,this->lid,IssmComm::GetRank()); break;
   }
 
   /*If this input is not already in Inputs, maybe it needs to be computed?*/
   switch(this->inputs2->GetInputObjectEnum(output_enum)){
      case TriaInput2Enum:
      case PentaInput2Enum:
      case TransientInput2Enum:{
         Input2* input2 = this->GetInput2(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 BoolInput2Enum:
         *pinterpolation   = P0Enum;
         *pnodesperelement = 1;
         *parray_size      = 1;
         break;
      case IntInput2Enum:
         *pinterpolation   = P0Enum;
         *pnodesperelement = 1;
         *parray_size      = 1;
         break;
      case ArrayInput2Enum:{
         int M;
         this->inputs2->GetArray(output_enum,this->lid,NULL,&M);
         *pinterpolation   = P0ArrayEnum;
         *pnodesperelement = 1;
         *parray_size      = M;
         }
         break;
      default:
         _error_("Input type \""<<EnumToStringx(this->inputs2->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*/
   Input2* input=this->GetInput2(output_enum);
   if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
 
   /*Cast to ElementInput*/
   if(input->ObjectEnum()!=TriaInput2Enum && input->ObjectEnum()!=PentaInput2Enum){
      _error_("Input "<<EnumToStringx(output_enum)<<" is not an ElementInput2");
   }
   ElementInput2* element_input = xDynamicCast<ElementInput2*>(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->inputs2->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->inputs2->GetInputObjectEnum(output_enum)){
      case TriaInput2Enum:
      case PentaInput2Enum:
      case TransientInput2Enum:{
 
         Input2* input2 = this->GetInput2(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 BoolInput2Enum:
         bool bvalue;
         this->GetInput2Value(&bvalue,output_enum);
         vector->SetValue(this->Sid(),reCast<IssmDouble>(bvalue),INS_VAL);
         break;
      case IntInput2Enum:
         int ivalue;
         this->GetInput2Value(&ivalue,output_enum);
         vector->SetValue(this->Sid(),reCast<IssmDouble>(ivalue),INS_VAL);
         break;
      default:
         _error_("Input type \""<<EnumToStringx(this->inputs2->GetInputObjectEnum(output_enum))<<"\" not supported yet");
   }
 
} /*}}}*/
void       Element::SetBoolInput(Inputs2* inputs2,int enum_in,bool value){/*{{{*/
 
   _assert_(inputs2);
   inputs2->SetInput(enum_in,this->lid,value);
 
}
/*}}}*/
void       Element::SetIntInput(Inputs2* inputs2,int enum_in,int value){/*{{{*/
 
   _assert_(inputs2);
   inputs2->SetInput(enum_in,this->lid,value);
 
}
/*}}}*/
void       Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
 
   /*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++){
 
      if(!flags[this->nodes[i]->Lid()]){
 
         /*flag current node so that no other element processes it*/
         flags[this->nodes[i]->Lid()]=true;
 
         int counter=0;
         while(flagsindices[counter]>=0) counter++;
         flagsindices[counter]=this->nodes[i]->Lid();
 
         /*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 analysis_type,approximation,numlayers;
   parameters->FindParam(&analysis_type,AnalysisTypeEnum);
   if(analysis_type==StressbalanceAnalysisEnum){
      this->GetInput2Value(&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: */
      Input2* dailysnowfall_input    = this->GetInput2(SmbDailysnowfallEnum,time_yr+(iday+1)/365.*yts); _assert_(dailysnowfall_input);
      Input2* dailyrainfall_input    = this->GetInput2(SmbDailyrainfallEnum,time_yr+(iday+1)/365.*yts); _assert_(dailyrainfall_input);
      Input2* dailydlradiation_input = this->GetInput2(SmbDailydlradiationEnum,time_yr+(iday+1)/365.*yts); _assert_(dailydlradiation_input);
      Input2* dailydsradiation_input = this->GetInput2(SmbDailydsradiationEnum,time_yr+(iday+1)/365.*yts); _assert_(dailydsradiation_input);
      Input2* dailywindspeed_input   = this->GetInput2(SmbDailywindspeedEnum,time_yr+(iday+1)/365.*yts); _assert_(dailywindspeed_input);
      Input2* dailypressure_input    = this->GetInput2(SmbDailypressureEnum,time_yr+(iday+1)/365.*yts); _assert_(dailypressure_input);
      Input2* dailyairdensity_input  = this->GetInput2(SmbDailyairdensityEnum,time_yr+(iday+1)/365.*yts); _assert_(dailyairdensity_input);
      Input2* dailyairhumidity_input = this->GetInput2(SmbDailyairhumidityEnum,time_yr+(iday+1)/365.*yts); _assert_(dailyairhumidity_input);
      Input2* dailytemperature_input = this->GetInput2(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->AddInput2(TemperatureSEMICEnum,&tsurf_out[0],P1Enum); // TODO add TemperatureSEMICEnum to EnumDefinitions
         this->AddInput2(SmbMassBalanceEnum,&smb_out[0],P1Enum);
         this->AddInput2(SmbAccumulationEnum,&saccu_out[0],P1Enum);
         this->AddInput2(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;
   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: {{{*/
   Input2 *zTop_input          = this->GetInput2(SmbZTopEnum);         _assert_(zTop_input);
   Input2 *dzTop_input         = this->GetInput2(SmbDzTopEnum);        _assert_(dzTop_input);
   Input2 *dzMin_input         = this->GetInput2(SmbDzMinEnum);        _assert_(dzMin_input);
   Input2 *zMax_input          = this->GetInput2(SmbZMaxEnum);         _assert_(zMax_input);
   Input2 *zMin_input          = this->GetInput2(SmbZMinEnum);         _assert_(zMin_input);
   Input2 *zY_input            = this->GetInput2(SmbZYEnum);           _assert_(zY_input);
   Input2 *Tmean_input         = this->GetInput2(SmbTmeanEnum);        _assert_(Tmean_input);
   Input2 *Vmean_input         = this->GetInput2(SmbVmeanEnum);        _assert_(Vmean_input);
   Input2 *C_input             = this->GetInput2(SmbCEnum);            _assert_(C_input);
   Input2 *Tz_input            = this->GetInput2(SmbTzEnum);           _assert_(Tz_input);
   Input2 *Vz_input            = this->GetInput2(SmbVzEnum);           _assert_(Vz_input);
   Input2 *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->inputs2->GetArray(SmbDziniEnum,this->lid,&dzini,&m);
      this->inputs2->GetArray(SmbDiniEnum,this->lid,&dini,&m);
      this->inputs2->GetArray(SmbReiniEnum,this->lid,&reini,&m);
      this->inputs2->GetArray(SmbGdniniEnum,this->lid,&gdnini,&m);
      this->inputs2->GetArray(SmbGspiniEnum,this->lid,&gspini,&m);
      this->inputs2->GetArray(SmbWiniEnum,this->lid,&Wini,&m);
      this->inputs2->GetArray(SmbAiniEnum,this->lid,&aini,&m);
      this->inputs2->GetArray(SmbTiniEnum,this->lid,&Tini,&m);
      EC_input = this->GetInput2(SmbECiniEnum);  _assert_(EC_input);
      EC_input->GetInputAverage(&EC);
 
      /*Retrieve the correct value of m (without the zeroes at the end)*/
      this->GetInput2Value(&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];
      }
 
      /*Flag the initialization:*/
      this->SetBoolInput(this->inputs2,SmbIsInitializedEnum,true);
   }
   else{
      /*Recover inputs: */
      this->inputs2->GetArray(SmbDzEnum,this->lid,&dz,&m);
      this->inputs2->GetArray(SmbDEnum,this->lid,&d,&m);
      this->inputs2->GetArray(SmbReEnum,this->lid,&re,&m);
      this->inputs2->GetArray(SmbGdnEnum,this->lid,&gdn,&m);
      this->inputs2->GetArray(SmbGspEnum,this->lid,&gsp,&m);
      this->inputs2->GetArray(SmbWEnum,this->lid,&W,&m);
      this->inputs2->GetArray(SmbAEnum,this->lid,&a,&m);
      this->inputs2->GetArray(SmbTEnum,this->lid,&T,&m);
      EC_input = this->GetInput2(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){
      Input2 *sumEC_input         = this->GetInput2(SmbECEnum);  _assert_(sumEC_input);
      Input2 *sumM_input          = this->GetInput2(SmbMeltEnum);  _assert_(sumM_input);
      Input2 *sumR_input          = this->GetInput2(SmbRunoffEnum);  _assert_(sumR_input);
      Input2 *sumP_input          = this->GetInput2(SmbPrecipitationEnum);  _assert_(sumP_input);
      Input2 *ULW_input           = this->GetInput2(SmbMeanULWEnum);  _assert_(ULW_input);
      Input2 *LW_input            = this->GetInput2(SmbNetLWEnum);  _assert_(LW_input);
      Input2 *SW_input            = this->GetInput2(SmbNetSWEnum);  _assert_(SW_input);
      Input2 *LHF_input           = this->GetInput2(SmbMeanLHFEnum);  _assert_(LHF_input);
      Input2 *SHF_input           = this->GetInput2(SmbMeanSHFEnum);  _assert_(SHF_input);
      Input2 *DzAdd_input         = this->GetInput2(SmbDzAddEnum);  _assert_(DzAdd_input);
      Input2 *MassAdd_input       = this->GetInput2(SmbMAddEnum);  _assert_(MassAdd_input);
      Input2 *InitMass_input      = this->GetInput2(SmbMInitnum);  _assert_(InitMass_input);
      Input2 *sumMsurf_input         = this->GetInput2(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
   Input2 *Ta_input  = this->GetInput2(SmbTaEnum,timeinputs);    _assert_(Ta_input);
   Input2 *V_input   = this->GetInput2(SmbVEnum,timeinputs);     _assert_(V_input);
   Input2 *Dlwr_input= this->GetInput2(SmbDlwrfEnum,timeinputs); _assert_(Dlwr_input);
   Input2 *Dswr_input= this->GetInput2(SmbDswrfEnum,timeinputs); _assert_(Dswr_input);
   Input2 *P_input   = this->GetInput2(SmbPEnum,timeinputs);     _assert_(P_input);
   Input2 *eAir_input= this->GetInput2(SmbEAirEnum,timeinputs);  _assert_(eAir_input);
   Input2 *pAir_input= this->GetInput2(SmbPAirEnum,timeinputs);  _assert_(pAir_input);
   Input2 *teValue_input= this->GetInput2(SmbTeValueEnum,timeinputs); _assert_(teValue_input);
   Input2 *aValue_input= this->GetInput2(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, 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 generated inputs: */
   this->inputs2->SetArrayInput(SmbDzEnum,this->lid,dz,m);
   this->inputs2->SetArrayInput(SmbDEnum,this->lid,d,m);
   this->inputs2->SetArrayInput(SmbReEnum,this->lid,re,m);
   this->inputs2->SetArrayInput(SmbGdnEnum,this->lid,gdn,m);
   this->inputs2->SetArrayInput(SmbGspEnum,this->lid,gsp,m);
   this->inputs2->SetArrayInput(SmbTEnum,this->lid,T,m);
   this->inputs2->SetArrayInput(SmbWEnum,this->lid,W,m);
   this->inputs2->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,Input2* vx_input,Input2* 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,Input2* vx_input,Input2* vy_input,Input2* 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,Input2* vx_input,Input2* 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,Input2* vx_input,Input2* 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,Input2* vx_input,Input2* 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,Input2* 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);
   Input2* sigma_xx_input  = this->GetInput2(StressTensorxxEnum); _assert_(sigma_xx_input);
   Input2* sigma_yy_input  = this->GetInput2(StressTensoryyEnum); _assert_(sigma_yy_input);
   Input2* sigma_xy_input  = this->GetInput2(StressTensorxyEnum); _assert_(sigma_xy_input);
   Input2* sigma_xz_input  = NULL;
   Input2* sigma_yz_input  = NULL;
   Input2* sigma_zz_input  = NULL;
   if(dim==3){
      sigma_xz_input  = this->GetInput2(StressTensorxzEnum); _assert_(sigma_xz_input);
      sigma_yz_input  = this->GetInput2(StressTensoryzEnum); _assert_(sigma_yz_input);
      sigma_zz_input  = this->GetInput2(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->AddInput2(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);
   Input2* vx_input = this->GetInput2(VxEnum); _assert_(vx_input);
   Input2* vy_input = this->GetInput2(VyEnum); _assert_(vy_input);
   Input2* vz_input = this->GetInput2(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->AddInput2(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);
}
/*}}}*/