Changeset 22955

Timestamp:

07/16/18 15:39:26 (7 years ago)

Author:

Eric.Larour

Message:

CHG: new reorganization of the slr capabilities.

Location:

issm/trunk-jpl/src

Files:

• c/Makefile.am (modified) (1 diff)
• c/analyses/MasstransportAnalysis.cpp (modified) (9 diffs)
• c/analyses/SealevelriseAnalysis.cpp (modified) (5 diffs)
• c/classes/Elements/Element.cpp (modified) (5 diffs)
• c/classes/Elements/Element.h (modified) (4 diffs)
• c/classes/Elements/Penta.h (modified) (1 diff)
• c/classes/Elements/Seg.h (modified) (1 diff)
• c/classes/Elements/Tetra.h (modified) (1 diff)
• c/classes/Elements/Tria.cpp (modified) (15 diffs)
• c/classes/Elements/Tria.h (modified) (1 diff)
• c/classes/FemModel.cpp (modified) (141 diffs)
• c/classes/FemModel.h (modified) (1 diff)
• c/classes/Materials/Matpar.h (modified) (1 diff)
• c/cores/cores.h (modified) (2 diffs)
• c/cores/sealevelrise_core.cpp (modified) (5 diffs)
• c/cores/sealevelrise_core_eustatic.cpp (deleted)
• c/cores/sealevelrise_core_noneustatic.cpp (deleted)
• c/cores/transient_core.cpp (modified) (1 diff)
• m/classes/maskpsl.m (modified) (2 diffs)
• m/classes/mesh3dsurface.m (modified) (1 diff)
• m/classes/model.js (modified) (2 diffs)
• m/classes/model.m (modified) (13 diffs)
• m/classes/model.py (modified) (2 diffs)
• m/classes/sealevelmodel.m (modified) (5 diffs)
• m/classes/slr.js (modified) (1 diff)
• m/classes/slr.m (modified) (15 diffs)
• m/classes/slr.py (modified) (4 diffs)

issm/trunk-jpl/src/c/Makefile.am

@@ -509 +509 @@
 if SEALEVELRISE
 issm_sources +=  ./cores/sealevelrise_core.cpp\
-                                 ./cores/sealevelrise_core_eustatic.cpp\
-                                 ./cores/sealevelrise_core_noneustatic.cpp\
                                  ./analyses/SealevelriseAnalysis.cpp
 endif

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

@@ -124 +124 @@
         bool   isoceancoupling;
         bool   issmb;
+        int    basalforcingsmodel;
 
         /*Fetch data needed: */
@@ -132 +133 @@
         iomodel->FindConstant(&isoceancoupling,"md.transient.isoceancoupling");
         iomodel->FindConstant(&issmb,"md.transient.issmb");
+        iomodel->FindConstant(&basalforcingsmodel,"md.basalforcings.model");
 
         /*Finite element type*/
@@ -156 +158 @@
         iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
         iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
+        if (basalforcingsmodel==SpatialLinearFloatingMeltRateEnum){
+                iomodel->FetchDataToInput(elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
+                iomodel->FetchDataToInput(elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
+                iomodel->FetchDataToInput(elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
+        }
         iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
+        
+        /*Initialize cumdeltalthickness input: unfortunately, we don't have femmodel, so we need to iterate on 
+         *elements, otherwise we would have used InputUpdateFromConstantx: */
+        counter=0;
+        for(int i=0;i<iomodel->numberofelements;i++){
+                if(iomodel->my_elements[i]){
+                        Element* element=(Element*)elements->GetObjectByOffset(counter);
+                        element->InputUpdateFromConstant(0.0,SealevelriseCumDeltathicknessEnum);
+                        counter++;
+                }
+        }
 
         /*Get what we need for ocean-induced basal melting*/
@@ -706 +724 @@
         basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
         IssmDouble* newthickness   = xNew<IssmDouble>(numnodes);
-        IssmDouble* deltathickness = xNew<IssmDouble>(numnodes);
+        IssmDouble* cumdeltathickness = xNew<IssmDouble>(numnodes);
+        IssmDouble* deltathickness    = xNew<IssmDouble>(numnodes);
         IssmDouble* newbase        = xNew<IssmDouble>(numnodes);
+        IssmDouble* bed            = xNew<IssmDouble>(numnodes);
         IssmDouble* newsurface     = xNew<IssmDouble>(numnodes);
         IssmDouble* oldthickness   = xNew<IssmDouble>(numnodes);
@@ -725 +745 @@
         }
 
-        /*Get previous base, thickness, surfac and current sealevel:*/
+        /*Get previous base, thickness, surfac and current sealevel and bed:*/
         basalelement->GetInputListOnNodes(&oldbase[0],BaseEnum);
         basalelement->GetInputListOnNodes(&oldsurface[0],SurfaceEnum);
@@ -731 +751 @@
         basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum);
         basalelement->GetInputListOnNodes(&sealevel[0],SealevelEnum);
-
-        /*What is the delta thickness forcing the sea-level rise core: */
-        for(i=0;i<numnodes;i++) deltathickness[i]=newthickness[i]-oldthickness[i];
+        basalelement->GetInputListOnNodes(&bed[0],BedEnum);
+        basalelement->GetInputListOnNodes(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum);
+
+        /*What is the delta thickness forcing the sea-level rise core: cumulated over time, hence the +=:*/
+        for(i=0;i<numnodes;i++){
+                cumdeltathickness[i]+=newthickness[i]-oldthickness[i];
+                deltathickness[i]=newthickness[i]-oldthickness[i];
+        }
 
         /*Find MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
@@ -742 +767 @@
         for(i=0;i<numnodes;i++) {
                 if (phi[i]>0.){ //this is an ice sheet: just add thickness to base.
-                        newsurface[i] = oldbase[i]+newthickness[i]; //surface = oldbase + newthickness
-                        newbase[i]     = oldbase[i];                 //same base: do nothing
+                        /*Update! actually, the bed has moved too!:*/
+                        newsurface[i] = bed[i]+newthickness[i]; //surface = bed + newthickness
+                        newbase[i]     = bed[i];                 //new base at new bed
                 }
                 else{ //this is an ice shelf: hydrostatic equilibrium*/
@@ -760 +786 @@
         /*Add input to the element: */
         element->AddBasalInput(ThicknessEnum,newthickness,P1Enum);
-        element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum);
         element->AddBasalInput(SurfaceEnum,newsurface,P1Enum);
         element->AddBasalInput(BaseEnum,newbase,P1Enum);
+        element->AddBasalInput(SealevelriseCumDeltathicknessEnum,cumdeltathickness,P1Enum);
+        element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum);
 
         /*Free ressources:*/
         xDelete<IssmDouble>(newthickness);
+        xDelete<IssmDouble>(cumdeltathickness);
         xDelete<IssmDouble>(newbase);
         xDelete<IssmDouble>(newsurface);
@@ -774 +802 @@
         xDelete<IssmDouble>(phi);
         xDelete<IssmDouble>(sealevel);
+        xDelete<IssmDouble>(bed);
+        
         xDelete<int>(doflist);
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp

@@ -19 +19 @@
 }/*}}}*/
 void SealevelriseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
+
+        int geodetic=0;
 
         /*Update elements: */
@@ -31 +33 @@
 
         /*Create inputs: */
+        iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
         iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
-        iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
-        iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum);
+        //those only if we have requested geodetic computations:
+        iomodel->FetchData(&geodetic,"md.slr.geodetic");
+        if (geodetic){
+                char* masktype=NULL;
+                iomodel->FetchData(&masktype,"md.mask.type");
+                if (strcmp(masktype,"maskpsl")==0){
+                        iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
+                        iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum);
+                }
+                xDelete<char>(masktype);
+        }
         iomodel->FetchDataToInput(elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum);
+        iomodel->FetchDataToInput(elements,"md.slr.spcthickness",SealevelriseSpcthicknessEnum);
         iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
+        iomodel->FetchDataToInput(elements,"md.geometry.bed",BedEnum);
+        iomodel->FetchDataToInput(elements,"md.slr.Ngia",SealevelNGiaRateEnum);
+        iomodel->FetchDataToInput(elements,"md.slr.Ugia",SealevelUGiaRateEnum);
+
+        /*Initialize cumdeltalthickness and sealevel rise rate input: unfortunately, we don't have femmodel, so we 
+         * need to iterate on elements, otherwise we would have used InputUpdateFromConstantx: */
+        counter=0;
+        for(int i=0;i<iomodel->numberofelements;i++){
+                if(iomodel->my_elements[i]){
+                        Element* element=(Element*)elements->GetObjectByOffset(counter);
+                        element->InputUpdateFromConstant(0.0,SealevelriseCumDeltathicknessEnum);
+                        element->InputUpdateFromConstant(0.0,SealevelNEsaRateEnum);
+                        element->InputUpdateFromConstant(0.0,SealevelUEsaRateEnum);
+                        element->InputUpdateFromConstant(0.0,SealevelRSLRateEnum);
+                        element->InputUpdateFromConstant(0.0,SealevelEustaticMaskEnum);
+                        element->InputUpdateFromConstant(0.0,SealevelEustaticOceanMaskEnum);
+                        counter++;
+                }
+        }
+        
         iomodel->FetchDataToInput(elements,"md.slr.steric_rate",SealevelriseStericRateEnum);
 
@@ -69 +102 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.abstol",SealevelriseAbstolEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.maxiter",SealevelriseMaxiterEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.slr.loop_increment",SealevelriseLoopIncrementEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.slr.horiz",SealevelriseHorizEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum));
@@ -79 +114 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum));
 
         iomodel->FetchData(&elastic,"md.slr.elastic");
@@ -249 +285 @@
 }/*}}}*/
 void           SealevelriseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-        
-        IssmDouble *deltaS  = NULL;
-        IssmDouble *S  = NULL;
-        int*        sidlist = NULL;
-        int         numvertices;
-        
-        numvertices= element->GetNumberOfVertices();
-        sidlist=xNew<int>(numvertices);
-        
-        element->GetVerticesSidList(sidlist);
-
-        deltaS = xNew<IssmDouble>(numvertices);
-        for(int i=0;i<numvertices;i++){
-                deltaS[i]=solution[sidlist[i]];
-        }
-
-        S = xNew<IssmDouble>(numvertices);
-        element->GetInputListOnVertices(S,SealevelEnum,0);
-
-        /*Add deltaS to S:*/
-        for (int i=0;i<numvertices;i++)S[i]+=deltaS[i];
-
-        /*Add S back into inputs: */
-        element->AddInput(SealevelEnum,S,P1Enum);
-
-        /*Free ressources:*/
-        xDelete<int>(sidlist);
-        xDelete<IssmDouble>(deltaS);
-        xDelete<IssmDouble>(S);
+        _error_("not implemeneted yet!");
 
 }/*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -100 +100 @@
                         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]);
@@ -1050 +1049 @@
                 gauss->GaussVertex(iv);
                 input->GetInputValue(&pvalue[iv],gauss);
+        }
+
+        /*clean-up*/
+        delete gauss;
+}
+/*}}}*/
+void       Element::GetInputListOnVerticesAtTime(IssmDouble* pvalue, int enumtype, IssmDouble time){/*{{{*/
+
+        /*Recover input*/
+        Input* input=this->GetInput(enumtype);
+        if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
+
+        /*Fetch number vertices for this element*/
+        int numvertices = this->GetNumberOfVertices();
+
+        /*Checks in debugging mode*/
+        _assert_(pvalue);
+
+        /* Start looping on the number of vertices: */
+        Gauss*gauss=this->NewGauss();
+        for(int iv=0;iv<numvertices;iv++){
+                gauss->GaussVertex(iv);
+                input->GetInputValue(&pvalue[iv],gauss,time);
         }
 
@@ -1333 +1355 @@
         }
 
+        /*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*/
+        int         numvertices = this->GetNumberOfVertices();
+        int         numnodes    = this->GetNumberOfNodes();
+        int*        doflist     = NULL;
+        IssmDouble* values      = NULL;
+
+        switch(type){
+        case VertexPIdEnum:
+                doflist = xNew<int>(numvertices);
+                values = xNew<IssmDouble>(numvertices);
+                /*Fill in values*/
+                this->GetVertexPidList(doflist);
+                this->GetInputListOnVerticesAtTime(values,input_enum,time);
+                vector->SetValues(numvertices,doflist,values,INS_VAL);
+                break;
+        case VertexSIdEnum:
+                doflist = xNew<int>(numvertices);
+                values = xNew<IssmDouble>(numvertices);
+                /*Fill in values*/
+                this->GetVerticesSidList(doflist);
+                this->GetInputListOnVerticesAtTime(values,input_enum,time);
+                vector->SetValues(numvertices,doflist,values,INS_VAL);
+                break;
+        default:
+                _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
+        }
+        
         /*Clean up*/
         xDelete<int>(doflist);
@@ -2039 +2096 @@
                                 name==MaterialsRheologyEsEnum ||
                                 name==MaterialsRheologyEsbarEnum ||
-                                name==SealevelEnum ||
-                                name==SealevelUmotionEnum ||
-                                name==SealevelNmotionEnum ||
-                                name==SealevelEmotionEnum ||
-                                name==SealevelAbsoluteEnum ||
-                                name==SealevelEustaticEnum ||
-                                name==SealevelriseDeltathicknessEnum ||
-                                name==EsaUmotionEnum ||
-                                name==EsaNmotionEnum ||
-                                name==EsaEmotionEnum ||
-                                name==EsaXmotionEnum ||
-                                name==EsaYmotionEnum ||
+                                name==SealevelEnum || 
+                                name==SealevelUEsaEnum || 
+                                name==SealevelUEsaRateEnum || 
+                                name==SealevelNEsaEnum || 
+                                name==SealevelNEsaRateEnum || 
+                                name==SealevelUNorthEsaEnum ||
+                                name==SealevelUEastEsaEnum ||
+                                name==SealevelRSLEustaticEnum || 
+                                name==SealevelRSLEustaticRateEnum || 
+                                name==SealevelRSLEnum || 
+                                name==SealevelRSLRateEnum || 
+                                name==SealevelEustaticMaskEnum || 
+                                name==SealevelEustaticOceanMaskEnum || 
+                                name==SealevelUGiaEnum || 
+                                name==SealevelUGiaRateEnum || 
+                                name==SealevelNGiaEnum || 
+                                name==SealevelNGiaRateEnum || 
+                                name==SealevelriseDeltathicknessEnum || 
+                                name==SealevelriseCumDeltathicknessEnum || 
+                                name==EsaUmotionEnum || 
+                                name==EsaNmotionEnum || 
+                                name==EsaEmotionEnum || 
                                 name==EsaStrainratexxEnum ||
                                 name==EsaStrainratexyEnum ||
@@ -2125 +2192 @@
         this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
         xDelete<IssmDouble>(base);
+        xDelete<IssmDouble>(values);
+
+}/*}}}*/
+void       Element::SpatialLinearFloatingiceMeltingRate(){/*{{{*/
+
+        int numvertices      = this->GetNumberOfVertices();
+        IssmDouble* deepwatermelt     = xNew<IssmDouble>(numvertices);
+        IssmDouble* deepwaterel     = xNew<IssmDouble>(numvertices);
+        IssmDouble* upperwaterel     = xNew<IssmDouble>(numvertices);
+        IssmDouble* base     = xNew<IssmDouble>(numvertices);
+        IssmDouble* values   = xNew<IssmDouble>(numvertices);
+
+        this->GetInputListOnVertices(base,BaseEnum);
+        this->GetInputListOnVertices(deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum);
+        this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum);
+        this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum);
+        
+        for(int i=0;i<numvertices;i++){
+                if(base[i]>upperwaterel[i])      values[i]=0;
+                else if (base[i]<deepwaterel[i]) values[i]=deepwatermelt[i];
+                else values[i]=deepwatermelt[i]*(base[i]-upperwaterel[i])/(deepwaterel[i]-upperwaterel[i]);
+        }
+
+        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
+        xDelete<IssmDouble>(base);
+        xDelete<IssmDouble>(deepwaterel);
+        xDelete<IssmDouble>(deepwatermelt);
+        xDelete<IssmDouble>(upperwaterel);
         xDelete<IssmDouble>(values);
 

issm/trunk-jpl/src/c/classes/Elements/Element.h

@@ -90 +90 @@
                 void               GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype);
                 void               GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
+                void               GetInputListOnVerticesAtTime(IssmDouble* pvalue,int enumtype,IssmDouble time);
                 void               GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
                 void               GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug);
@@ -104 +105 @@
                 void               GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int solutionenum);
                 void               GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum, int type);
+                void               GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum, int type,IssmDouble time);
                 void                 GetVertexPidList(int* pidlist);
                 void               GetVerticesConnectivityList(int* connectivitylist);
@@ -134 +136 @@
                 bool               IsWaterInElement();
                 void               LinearFloatingiceMeltingRate(); 
+                void               SpatialLinearFloatingiceMeltingRate(); 
                 void               MantlePlumeGeothermalFlux(); 
                 void               MarshallElement(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction,int numanalyses);
@@ -336 +339 @@
                 virtual void          SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
                 virtual void          SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0;
-                virtual void          SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0;
+                virtual void          SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0;
                 #endif
 

issm/trunk-jpl/src/c/classes/Elements/Penta.h

@@ -201 +201 @@
                 void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
                 void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
-                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
+                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
                 #endif
 

issm/trunk-jpl/src/c/classes/Elements/Seg.h

@@ -182 +182 @@
                 void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
                 void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
-                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
+                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
                 IssmDouble    OceanArea(void){_error_("not implemented yet!");};
                 IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};

issm/trunk-jpl/src/c/classes/Elements/Tetra.h

@@ -189 +189 @@
                 void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
                 void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
-                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
+                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
                 IssmDouble    OceanArea(void){_error_("not implemented yet!");};
                 IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -1354 +1354 @@
 
         }
-        else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum){
+        else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
                 /*Check that not all nodes are grounded or floating*/
                 if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
@@ -1905 +1905 @@
 
         if(!IsIceInElement())return 0.;
+        if(!(this->inputs->Max(MaskIceLevelsetEnum)<0)) return 0;
 
         int domaintype;
@@ -4676 +4677 @@
 /*}}}*/
 void    Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){/*{{{*/
-
         /*early return if we are not on an ice cap OR ocean:*/
         if(!(this->inputs->Max(MaskIceLevelsetEnum)<0) && !IsWaterInElement()){
@@ -4778 +4778 @@
         IssmDouble late,longe,re;
         IssmDouble lati,longi,ri;
+        bool notfullygrounded=false;
 
         /*elastic green function:*/
@@ -4785 +4786 @@
         /*ice properties: */
         IssmDouble rho_ice,rho_water,rho_earth;
+
+        /*constants:*/
+        IssmDouble constant=0;
 
         /*Initialize eustatic component: do not skip this step :):*/
@@ -4795 +4799 @@
         
         /*early return if we are not on an ice cap:*/
-        if(!(this->inputs->Max(MaskIceLevelsetEnum)<0)){
+        if(!(this->inputs->Max(MaskIceLevelsetEnum)<=0)){
+                constant=0; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum));
                 *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
                 return;
         }
+
+        /*early return if we are fully floating: */
+        if (this->inputs->Max(MaskGroundediceLevelsetEnum)<=0){
+                constant=0; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum));
+                *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
+                return;
+        }
+        
+        /*If we are here, we are on ice that is fully grounded or half-way to floating: */
+        if ((this->inputs->Min(MaskGroundediceLevelsetEnum))<0){
+                notfullygrounded=true; //used later on.
+        }
+                
+        /*Inform mask: */
+        constant=1; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum));
 
         /*recover material parameters: */
@@ -4859 +4879 @@
         /*}}}*/
 
-        /*Compute area of element:*/
+        /*Compute area of element. Scale it by grounded fraction if not fully grounded: */
         area=GetAreaSpherical();
-
+        if(notfullygrounded){
+                IssmDouble phi=0;
+                IssmDouble xyz_list[NUMVERTICES][3];
+                ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+                phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
+                area*=phi;
+        }
+        
         /*Compute ice thickness change: */
         Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); 
         if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
-        deltathickness_input->GetInputAverage(&I);
+
+        /*If we are fully grounded, take the average over the element: */
+        if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
+        else{ 
+                IssmDouble total_weight=0;
+                bool mainlyfloating = true;
+                int         point1;
+                IssmDouble  fraction1,fraction2;
+                
+                /*Recover portion of element that is grounded*/
+                this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+                Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
+
+                /* Start  looping on the number of gaussian points and average over these gaussian points: */
+                total_weight=0;
+                I=0;
+                for(int ig=gauss->begin();ig<gauss->end();ig++){
+                        IssmDouble Ig=0;
+                        gauss->GaussPoint(ig);
+                        deltathickness_input->GetInputValue(&Ig,gauss);
+                        I+=Ig*gauss->weight;
+                        total_weight+=gauss->weight;
+                }
+                I=I/total_weight;
+                delete gauss;
+        }
 
         /*Compute eustatic compoent:*/
@@ -4927 +4980 @@
         IssmDouble minlong=400;
         IssmDouble maxlong=-20;
+        IssmDouble constant=0;
 
         /*precomputed elastic green functions:*/
@@ -4947 +5001 @@
 
         /*early return if we are not on the ocean:*/
-        if (!IsWaterInElement())return;
+        if (!IsWaterInElement()){
+                constant=0; this->inputs->AddInput(new TriaInput(SealevelEustaticOceanMaskEnum,&constant,P0Enum));
+                return;
+        }
+        constant=1; this->inputs->AddInput(new TriaInput(SealevelEustaticOceanMaskEnum,&constant,P0Enum));
 
         /*recover computational flags: */
@@ -5058 +5116 @@
 }
 /*}}}*/
-void    Tria::SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/
+void    Tria::SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/
 
         /*diverse:*/
@@ -5082 +5140 @@
         IssmDouble* N_elastic= NULL;
         IssmDouble* E_elastic= NULL;
+        DoubleVecParam* U_parameter = NULL;
+        DoubleVecParam* H_parameter = NULL;
+        IssmDouble* U_values=NULL;
+        IssmDouble* N_values=NULL;
+        IssmDouble* E_values=NULL;
 
         /*optimization:*/
@@ -5092 +5155 @@
         /*early return if we are not on the ocean or on an ice cap:*/
         if(!(this->inputs->Max(MaskIceLevelsetEnum)<0) && !IsWaterInElement()) return; 
+
+        /*early return if we are fully floating: */
+        if (this->inputs->Max(MaskGroundediceLevelsetEnum)<=0)return;
 
         /*recover computational flags: */
@@ -5156 +5222 @@
 
         /*recover elastic Green's functions for displacement:*/
-        DoubleVecParam* U_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter);
-        DoubleVecParam* H_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter);
+        U_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter);
         U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M);
-        H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);
+        if(horiz){
+                H_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter);
+                H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);
+        }
 
         /*From Sg, recover water sea level rise:*/
@@ -5171 +5239 @@
         /*initialize: */
         U_elastic=xNewZeroInit<IssmDouble>(gsize);
-        N_elastic=xNewZeroInit<IssmDouble>(gsize);
-        E_elastic=xNewZeroInit<IssmDouble>(gsize);
+        if(horiz){
+                N_elastic=xNewZeroInit<IssmDouble>(gsize);
+                E_elastic=xNewZeroInit<IssmDouble>(gsize);
+        }
 
         int* indices=xNew<int>(gsize);
-        IssmDouble* U_values=xNewZeroInit<IssmDouble>(gsize);
-        IssmDouble* N_values=xNewZeroInit<IssmDouble>(gsize);
-        IssmDouble* E_values=xNewZeroInit<IssmDouble>(gsize);
+        U_values=xNewZeroInit<IssmDouble>(gsize);
+        if(horiz){
+                N_values=xNewZeroInit<IssmDouble>(gsize); 
+                E_values=xNewZeroInit<IssmDouble>(gsize);
+        }
         IssmDouble alpha;
         IssmDouble delPhi,delLambda;
@@ -5202 +5274 @@
                 }
                 dx = x_element-x; dy = y_element-y; dz = z_element-z; 
-                N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-                E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
+                if(horiz){
+                        N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
+                        E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
+                }
                 
                 /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
                 int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
                 U_elastic[i] += U_elastic_precomputed[index];
-                N_elastic[i] += H_elastic_precomputed[index]*N_azim;
-                E_elastic[i] += H_elastic_precomputed[index]*E_azim;
+                if(horiz){
+                        N_elastic[i] += H_elastic_precomputed[index]*N_azim;
+                        E_elastic[i] += H_elastic_precomputed[index]*E_azim;
+                }
 
                 /*Add all components to the pUp solution vectors:*/
                 if(this->inputs->Max(MaskIceLevelsetEnum)<0){
                         U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i];
-                        N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i];
-                        E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i];
+                        if(horiz){
+                                N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i];
+                                E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i];
+                        }
                 }
                 else if(IsWaterInElement()) {
                         U_values[i]+=3*rho_water/rho_earth*area/eartharea*S*U_elastic[i];
-                        N_values[i]+=3*rho_water/rho_earth*area/eartharea*S*N_elastic[i];
-                        E_values[i]+=3*rho_water/rho_earth*area/eartharea*S*E_elastic[i];
+                        if(horiz){
+                                N_values[i]+=3*rho_water/rho_earth*area/eartharea*S*N_elastic[i];
+                                E_values[i]+=3*rho_water/rho_earth*area/eartharea*S*E_elastic[i];
+                        }
                 }
         }
         pUp->SetValues(gsize,indices,U_values,ADD_VAL);
-        pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
-        pEast->SetValues(gsize,indices,E_values,ADD_VAL);
+        if(horiz){
+                pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
+                pEast->SetValues(gsize,indices,E_values,ADD_VAL);
+        }
 
         /*free ressources:*/
         xDelete<int>(indices); 
-        xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
-        xDelete<IssmDouble>(U_elastic); xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
+        xDelete<IssmDouble>(U_values); 
+        xDelete<IssmDouble>(U_elastic); 
+        if(horiz){
+                xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
+                xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
+        }
 
         return;

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -162 +162 @@
                 void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
                 void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea);
-                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea);
+                void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz);
                 #endif
                 /*}}}*/

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -1454 +1454 @@
         /*Figure out maximum across the cluster: */
         ISSM_MPI_Reduce(&maxvy,&node_maxvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&node_maxvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
         maxvy=node_maxvy;
 
@@ -1478 +1478 @@
         /*Figure out maximum across the cluster: */
         ISSM_MPI_Reduce(&maxvz,&node_maxvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&node_maxvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
         maxvz=node_maxvz;
 
@@ -1502 +1502 @@
         /*Figure out minimum across the cluster: */
         ISSM_MPI_Reduce(&minvel,&node_minvel,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&node_minvel,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
         minvel=node_minvel;
 
@@ -1526 +1526 @@
         /*Figure out minimum across the cluster: */
         ISSM_MPI_Reduce(&minvx,&node_minvx,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&node_minvx,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
         minvx=node_minvx;
 
@@ -1550 +1550 @@
         /*Figure out minimum across the cluster: */
         ISSM_MPI_Reduce(&minvy,&node_minvy,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&node_minvy,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
         minvy=node_minvy;
 
@@ -1574 +1574 @@
         /*Figure out minimum across the cluster: */
         ISSM_MPI_Reduce(&minvz,&node_minvz,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&node_minvz,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());   
         minvz=node_minvz;
 
@@ -1619 +1619 @@
                         omega_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
 
-                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
+                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
                         //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
                         J+=weight*1/2*pow(dp[0]*dp[0]+dp[1]*dp[1],2)*Jdet*gauss->weight;
@@ -1955 +1955 @@
                                                                 IssmDouble* values    = xNewZeroInit<IssmDouble>(nlines*ncols);
                                                                 IssmDouble* allvalues = xNew<IssmDouble>(nlines*ncols);
-
+                                                                
                                                                 /*Fill-in matrix*/
                                                                 for(int j=0;j<elements->Size();j++){
@@ -1964 +1964 @@
                                                                 ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
                                                                 xDelete<IssmDouble>(values);
-
+                                                                
                                                                 if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time));
                                                                 xDelete<IssmDouble>(allvalues);
@@ -2012 +2012 @@
         int domaintype;
         this->parameters->FindParam(&domaintype,DomainTypeEnum);
-
+        
         /*1: go throug all elements of this partition and figure out how many
          * segments we have (corresopnding to levelset = 0)*/
@@ -2132 +2132 @@
                 case VelEnum:                            this->ElementResponsex(responses,VelEnum); break;
                 case FrictionCoefficientEnum:            NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break;
-                default:
+                default: 
                         if(response_descriptor_enum>=Outputdefinition1Enum && response_descriptor_enum <=Outputdefinition100Enum){
                                 IssmDouble double_result = OutputDefinitionsResponsex(this,response_descriptor_enum);
                                 *responses=double_result;
                         }
-                        else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!");
-                        break;
+                        else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); 
+                        break; 
         }
 
@@ -2269 +2269 @@
                         thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
 
-                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
+                        /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
                         J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight;
                 }
@@ -2460 +2460 @@
         analysis->UpdateConstraints(this);
         delete analysis;
-
+        
         /*Second, constraints might be time dependent: */
-        SpcNodesx(nodes,constraints,parameters,analysis_type);
+        SpcNodesx(nodes,constraints,parameters,analysis_type); 
 
         /*Now, update degrees of freedoms: */
@@ -2473 +2473 @@
         IssmDouble         *surface = NULL;
         IssmDouble         *bed     = NULL;
-
+                        
         if(VerboseSolution()) _printf0_("   updating vertices positions\n");
 
@@ -2512 +2512 @@
 
 /*AMR*/
-#if !defined(_HAVE_ADOLC_)
+#if !defined(_HAVE_ADOLC_) 
 void FemModel::ReMesh(void){/*{{{*/
 
@@ -2522 +2522 @@
         int newnumberofvertices = -1;
         int newnumberofelements = -1;
-        bool* my_elements                       = NULL;
+        bool* my_elements                       = NULL; 
         int* my_vertices                        = NULL;
         int elementswidth       = this->GetElementsWidth();//just tria elements in this version
@@ -2528 +2528 @@
         bool isgroundingline;
 
-        /*Branch to specific amr depending on requested method*/
+        /*Branch to specific amr depending on requested method*/        
         this->parameters->FindParam(&amrtype,AmrTypeEnum);
         switch(amrtype){
@@ -2541 +2541 @@
                 default: _error_("not implemented yet");
         }
-
+        
         /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/
         this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices);
@@ -2572 +2572 @@
 
                 /*As the domain is 2D, it is not necessary to create nodes for this analysis*/
-                if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;
+                if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue;    
 
                 this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes);
@@ -2602 +2602 @@
                 //SetCurrentConfiguration(analysis_type);
 
-                this->analysis_counter=i;
+                this->analysis_counter=i;       
                 /*Now, plug analysis_counter and analysis_type inside the parameters: */
                 this->parameters->SetParam(this->analysis_counter,AnalysisCounterEnum);
@@ -2619 +2619 @@
 
                 ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters);
-                if(i==0){
+                if(i==0){ 
                         VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices
                 }
@@ -2663 +2663 @@
         /*Insert bedrock from mismip+ setup*/
         /*This was used to Misomip project/simulations*/
-
+        
         if(VerboseSolution())_printf0_("        call Mismip bedrock adjust module\n");
 
@@ -2680 +2680 @@
                         by              = 500./(1.+exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+exp((2./4000.)*(y-80000./2.+24000.)));
                         r[i]    = max(bx+by,-720.);
-                }
+                }       
                 /*insert new bedrock*/
                 element->AddInput(BedEnum,&r[0],P1Enum);
@@ -2693 +2693 @@
 
         if(VerboseSolution())_printf0_("        call adjust base and thickness module\n");
-
+        
         int     numvertices = this->GetElementsWidth();
    IssmDouble rho_water,rho_ice,density,base_float;
@@ -2705 +2705 @@
         for(int el=0;el<this->elements->Size();el++){
                 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(el));
-
+        
                 element->GetInputListOnVertices(&s[0],SurfaceEnum);
                 element->GetInputListOnVertices(&r[0],BedEnum);
@@ -2717 +2717 @@
                         base_float = rho_ice*s[i]/(rho_ice-rho_water);
                         if(r[i]>base_float){
-                                b[i] = r[i];
-                        }
+                                b[i] = r[i];                    
+                        } 
                         else {
-                                b[i] = base_float;
-                        }
+                                b[i] = base_float;      
+                        } 
 
                         if(fabs(sl[i])>0) _error_("Sea level value "<<sl[i]<<" not supported!");
                         /*update thickness and mask grounded ice level set*/
                         h[i]      = s[i]-b[i];
-                        phi[i]  = h[i]+r[i]/density;
+                        phi[i]  = h[i]+r[i]/density;    
                 }
 
@@ -2734 +2734 @@
                 element->AddInput(BaseEnum,&b[0],P1Enum);
         }
-
+        
    /*Delete*/
    xDelete<IssmDouble>(phi);
@@ -2762 +2762 @@
         Vector<IssmDouble>* input_interpolations        = NULL;
         IssmDouble* input_interpolations_serial = NULL;
-   int* pos                                                                                             = NULL;
+   int* pos                                                                                             = NULL; 
         IssmDouble value                                                                        = 0;
 
@@ -2782 +2782 @@
                         P0input_interp = xNew<int>(numP0inputs);
                         P1input_enums  = xNew<int>(numP1inputs);
-                        P1input_interp = xNew<int>(numP1inputs);
+                        P1input_interp = xNew<int>(numP1inputs);        
                 }
                 numP0inputs = 0;
@@ -2814 +2814 @@
                 }
         }
-
-        /*Get P0 and P1 inputs over the elements*/
+        
+        /*Get P0 and P1 inputs over the elements*/      
         pos             = xNew<int>(elementswidth);
         vP0inputs= new Vector<IssmDouble>(numberofelements*numP0inputs);
@@ -2821 +2821 @@
         for(int i=0;i<this->elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
-
+                
                 /*Get P0 inputs*/
                 for(int j=0;j<numP0inputs;j++){
-                        TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j]));
+                        TriaInput* input=xDynamicCast<TriaInput*>(element->GetInput(P0input_enums[j]));         
                         input->GetInputAverage(&value);
                         pos[0]=element->Sid()*numP0inputs+j;
-                        /*Insert input in the vector*/
+                        /*Insert input in the vector*/  
                         vP0inputs->SetValues(1,pos,&value,INS_VAL);
                 }
-
+                
                 /*Get P1 inputs*/
                 for(int j=0;j<numP1inputs;j++){
@@ -2837 +2837 @@
                         pos[1]=element->vertices[1]->Sid()*numP1inputs+j;
                         pos[2]=element->vertices[2]->Sid()*numP1inputs+j;
-                        /*Insert input in the vector*/
-                        vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL);
+                        /*Insert input in the vector*/  
+                        vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL);  
                 }
         }
@@ -2847 +2847 @@
         P0inputs=vP0inputs->ToMPISerial();
         P1inputs=vP1inputs->ToMPISerial();
-
+        
         /*Assign pointers*/
-        *pnumP0inputs           = numP0inputs;
-        *pP0inputs                      = P0inputs;
+        *pnumP0inputs           = numP0inputs; 
+        *pP0inputs                      = P0inputs; 
         *pP0input_enums = P0input_enums;
         *pP0input_interp        = P0input_interp;
-        *pnumP1inputs           = numP1inputs;
-        *pP1inputs                      = P1inputs;
+        *pnumP1inputs           = numP1inputs; 
+        *pP1inputs                      = P1inputs; 
         *pP1input_enums = P1input_enums;
-        *pP1input_interp        = P1input_interp;
+        *pP1input_interp        = P1input_interp;       
 
         /*Cleanup*/
@@ -2867 +2867 @@
 /*}}}*/
 void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/
-
+        
         int numberofelements                    = this->elements->NumberOfElements();   //global, entire old mesh
         int newnumberofelements         = newfemmodel_elements->Size();                 //just on the new partition
@@ -2883 +2883 @@
         int* P1input_enums                      = NULL;
         int* P1input_interp                     = NULL;
-        IssmDouble* values                      = NULL;
+        IssmDouble* values                      = NULL; 
    IssmDouble* vector           = NULL;
         IssmDouble* x                                   = NULL;//global, entire old mesh
@@ -2895 +2895 @@
         IssmDouble* newyc                               = NULL;//just on the new partition
         int* newelementslist                    = NULL;//just on the new partition
-        int* sidtoindex                         = NULL;//global vertices sid to partition index
+        int* sidtoindex                         = NULL;//global vertices sid to partition index 
 
         /*Get old P0 and P1  inputs (entire mesh)*/
@@ -2928 +2928 @@
 
         /*Insert P0 and P1 inputs into the new elements (just on the new partition)*/
-        values=xNew<IssmDouble>(elementswidth);
+        values=xNew<IssmDouble>(elementswidth); 
         for(int i=0;i<newfemmodel_elements->Size();i++){//just on the new partition
                 Element* element=xDynamicCast<Element*>(newfemmodel_elements->GetObjectByOffset(i));
                 /*newP0inputs is just on the new partition*/
                 for(int j=0;j<numP0inputs;j++){
-                        switch(P0input_interp[j]){
+                        switch(P0input_interp[j]){      
                                 case P0Enum:
                                 case DoubleInputEnum:
                                         element->AddInput(new DoubleInput(P0input_enums[j],newP0inputs[i*numP0inputs+j]));
                                         break;
-                                case IntInputEnum:
+                                case IntInputEnum: 
                                         element->AddInput(new IntInput(P0input_enums[j],reCast<int>(newP0inputs[i*numP0inputs+j])));
                                         break;
@@ -2956 +2956 @@
                 }
         }
-
+        
         /*Cleanup*/
         xDelete<IssmDouble>(P0inputs);
@@ -2995 +2995 @@
 
         if(!this->elements || !this->vertices || !this->results || !this->parameters) return;
-
+         
         parameters->FindParam(&step,StepEnum);
         parameters->FindParam(&time,TimeEnum);
@@ -3013 +3013 @@
         this->results->AddResult(new GenericExternalResult<IssmDouble*>(this->results->Size()+1,MeshYEnum,
                                         y,numberofvertices,1,step,time));
-
+        
         /*Cleanup*/
         xDelete<IssmDouble>(x);
@@ -3074 +3074 @@
    /*Assemble*/
         vmasklevelset->Assemble();
-
+        
         /*Serialize and set output*/
         (*pmasklevelset)=vmasklevelset->ToMPISerial();
@@ -3088 +3088 @@
 
         /*newelementslist is in Matlab indexing*/
-
+        
         /*Creating connectivity table*/
         int* connectivity=NULL;
@@ -3099 +3099 @@
                         connectivity[vertexid-1]+=1;//Matlab to C indexing
                 }
-        }
+        }       
 
         /*Create vertex and insert in vertices*/
         for(int i=0;i<newnumberofvertices;i++){
                 if(my_vertices[i]){
-                        Vertex *newvertex=new Vertex();
+                        Vertex *newvertex=new Vertex(); 
                         newvertex->id=i+1;
                         newvertex->sid=i;
@@ -3115 +3115 @@
                         newvertex->connectivity=connectivity[i];
                         newvertex->clone=false;//itapopo check this
-                        vertices->AddObject(newvertex);
-                }
+                        vertices->AddObject(newvertex); 
+                } 
         }
 
@@ -3143 +3143 @@
                         }
                         else newtria->element_type_list=NULL;
-
+                        
                         /*Element hook*/
                         int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials)
@@ -3149 +3149 @@
                         /*retrieve vertices ids*/
                         int* vertex_ids=xNew<int>(elementswidth);
-                        for(int j=0;j<elementswidth;j++)        vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing
+                        for(int j=0;j<elementswidth;j++)        vertex_ids[j]=reCast<int>(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing 
                         /*Setting the hooks*/
                         newtria->numanalyses =this->nummodels;
@@ -3161 +3161 @@
                         /*Clean up*/
                         xDelete<int>(vertex_ids);
-                        elements->AddObject(newtria);
-                }
+                        elements->AddObject(newtria);   
+                } 
         }
 
@@ -3172 +3172 @@
         for(int i=0;i<newnumberofelements;i++){
                 if(my_elements[i]){
-                        materials->AddObject(new Matice(i+1,i,MaticeEnum));
-                }
-        }
-
+                        materials->AddObject(new Matice(i+1,i,MaticeEnum));     
+                } 
+        }
+        
         /*Add new constant material property to materials, at the end: */
         Matpar *newmatpar=static_cast<Matpar*>(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy());
         newmatpar->SetMid(newnumberofelements+1);
-        materials->AddObject(newmatpar);//put it at the end of the materials
+        materials->AddObject(newmatpar);//put it at the end of the materials        
 }
 /*}}}*/
@@ -3186 +3186 @@
         int lid=0;
         for(int j=0;j<newnumberofvertices;j++){
-                if(my_vertices[j]){
-
-                        Node* newnode=new Node();
-
+                if(my_vertices[j]){                             
+                        
+                        Node* newnode=new Node();       
+                        
                         /*id: */
                         newnode->id=nodecounter+j+1;
@@ -3195 +3195 @@
                         newnode->lid=lid++;
                         newnode->analysis_enum=analysis_enum;
-
+                        
                         /*Initialize coord_system: Identity matrix by default*/
                         for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0;
                         for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0;
-
+                        
                         /*indexing:*/
                         newnode->indexingupdate=true;
-
+                        
                         Analysis* analysis=EnumToAnalysis(analysis_enum);
                         int *doftypes=NULL;
@@ -3216 +3216 @@
                         /*Stressbalance Horiz*/
                         if(analysis_enum==StressbalanceAnalysisEnum){
-                                // itapopo this code is rarely used.
+                                // itapopo this code is rarely used. 
                                 /*Coordinate system provided, convert to coord_system matrix*/
                                 //XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]);
@@ -3231 +3231 @@
         if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
         if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
-
+        
         int numberofvertices = femmodel_vertices->NumberOfVertices();
         int numberofelements = femmodel_elements->NumberOfElements();
@@ -3237 +3237 @@
         IssmDouble* x                   = NULL;
         IssmDouble* y                   = NULL;
-        IssmDouble* z                   = NULL;
+        IssmDouble* z                   = NULL; 
         int* elementslist       = NULL;
         int* elem_vertices      = NULL;
@@ -3246 +3246 @@
         /*Get vertices coordinates*/
         VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ;
-
+        
         /*Get element vertices*/
         elem_vertices                           = xNew<int>(elementswidth);
@@ -3261 +3261 @@
         vid3->SetValue(element->sid,elem_vertices[2],INS_VAL);
    }
-
+                
         /*Assemble*/
    vid1->Assemble();
@@ -3271 +3271 @@
    id2 = vid2->ToMPISerial();
         id3 = vid3->ToMPISerial();
-
+        
         /*Construct elements list*/
         elementslist=xNew<int>(numberofelements*elementswidth);
@@ -3280 +3280 @@
                 elementslist[elementswidth*i+2] = reCast<int>(id3[i])+1; //InterpMesh wants Matlab indexing
         }
-
+        
         /*Assign pointers*/
         *px                             = x;
@@ -3301 +3301 @@
         if(!femmodel_vertices) _error_("GetMesh: vertices are NULL.");
         if(!femmodel_elements) _error_("GetMesh: elements are NULL.");
-
+        
         int numberofvertices                    = femmodel_vertices->Size();    //number of vertices of this partition
         int numbertotalofvertices       = femmodel_vertices->NumberOfVertices();        //number total of vertices (entire mesh)
@@ -3308 +3308 @@
         IssmDouble* x                                   = NULL;
         IssmDouble* y                                   = NULL;
-        IssmDouble* z                                   = NULL;
+        IssmDouble* z                                   = NULL; 
         int* elementslist                               = NULL;
         int* sidtoindex                         = NULL;
         int* elem_vertices                      = NULL;
-
+        
         /*Get vertices coordinates of this partition*/
         sidtoindex      = xNewZeroInit<int>(numbertotalofvertices);//entire mesh, all vertices
@@ -3318 +3318 @@
         y                               = xNew<IssmDouble>(numberofvertices);//just this partitio;
         z                               = xNew<IssmDouble>(numberofvertices);//just this partitio;
-
+        
         /*Go through in this partition (vertices)*/
         for(int i=0;i<numberofvertices;i++){//just this partition
-                Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i);
+                Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i);        
                 /*Attention: no spherical coordinates*/
                 x[i]=vertex->GetX();
@@ -3334 +3334 @@
         elementslist = xNew<int>(numberofelements*elementswidth);
         if(numberofelements*elementswidth<0) _error_("numberofelements negative.");
-
+        
         for(int i=0;i<numberofelements;i++){//just this partition
         Element* element=xDynamicCast<Element*>(femmodel_elements->GetObjectByOffset(i));
@@ -3341 +3341 @@
                 elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing
                 elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing
-        }
-
+        }       
+                
         /*Assign pointers*/
         *px                             = x;
@@ -3348 +3348 @@
         *pz                             = z;
         *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type.
-        *psidtoindex    = sidtoindex;  //it is ncessary to insert inputs
+        *psidtoindex    = sidtoindex;  //it is ncessary to insert inputs 
 
         /*Cleanup*/
@@ -3359 +3359 @@
         /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/
         if(analysis_enum!=StressbalanceAnalysisEnum) return;
-
+        
         int numberofnodes_analysistype= this->nodes->NumberOfNodes(analysis_enum);
-        int dofpernode                                          = 2;                                                                                                            //vx and vy
+        int dofpernode                                          = 2;                                                                                                            //vx and vy 
         int numberofcols                                        = dofpernode*2;                                                                         //to keep dofs and flags in the vspc vector
         int numberofvertices                            = this->vertices->NumberOfVertices();                   //global, entire old mesh
@@ -3385 +3385 @@
         newy=xNew<IssmDouble>(newnumberofvertices);//just the new partition
         for(int i=0;i<newnumberofvertices;i++){//just the new partition
-                Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);
+                Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i);     
                 /*Attention: no spherical coordinates*/
                 newx[i]=vertex->GetX();
@@ -3393 +3393 @@
         /*Get spcvx and spcvy of old mesh*/
         for(int i=0;i<this->constraints->Size();i++){
-
+                
                 Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i);
                 if(!constraint->InAnalysis(analysis_enum)) _error_("AMR create constraints for "<<EnumToStringx(analysis_enum)<<" not supported yet!\n");
@@ -3400 +3400 @@
                 int dof                                 = spcstatic->GetDof();
                 int node                                        = spcstatic->GetNodeId();
-                IssmDouble spcvalue     = spcstatic->GetValue();
+                IssmDouble spcvalue     = spcstatic->GetValue(); 
                 int nodeindex                   = node-1;
-
+                
                 /*vx and vx flag insertion*/
                 if(dof==0) {//vx
                         vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL);    //vx
                         vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag
-                }
+                } 
                 /*vy and vy flag insertion*/
                 if(dof==1){//vy
@@ -3423 +3423 @@
                                                                 spc,numberofvertices,numberofcols,
                                                                 newx,newy,newnumberofvertices,NULL);
-
+        
         /*Now, insert the interpolated constraints in the data set (constraints)*/
         count=0;
@@ -3440 +3440 @@
                 /*spcvy*/
                 if(!xIsNan<IssmDouble>(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){
-                        newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum));
+                        newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); 
                         //add count'th spc, on node i+1, setting dof 1 to vx.
                         count++;
@@ -3497 +3497 @@
         bool *my_elements = NULL;
         int *my_vertices  = NULL;
-
-        _assert_(newnumberofvertices>0);
-        _assert_(newnumberofelements>0);
+        
+        _assert_(newnumberofvertices>0); 
+        _assert_(newnumberofelements>0); 
         epart=xNew<int>(newnumberofelements);
         npart=xNew<int>(newnumberofvertices);
    index=xNew<int>(elementswidth*newnumberofelements);
-
+   
         for (int i=0;i<newnumberofelements;i++){
         for (int j=0;j<elementswidth;j++){
@@ -3523 +3523 @@
                 for (int i=0;i<newnumberofvertices;i++) npart[i]=0;
         }
-        else _error_("At least one processor is required");
+        else _error_("At least one processor is required");         
 
         my_vertices=xNew<int>(newnumberofvertices);
@@ -3533 +3533 @@
         for(int i=0;i<newnumberofelements;i++){
                 /*!All elements have been partitioned above, only deal with elements for this cpu: */
-                if(my_rank==epart[i]){
+                if(my_rank==epart[i]){ 
                         my_elements[i]=true;
-                        /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use
-                         *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
-                         into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices
+                        /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use 
+                         *the  element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) 
+                         into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices 
                          will hold which vertices belong to this partition*/
                         for(int j=0;j<elementswidth;j++){
-                                _assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing
+                                _assert_(newelementslist[elementswidth*i+j]-1<newnumberofvertices);//newelementslist is in Matlab indexing 
                                 my_vertices[newelementslist[elementswidth*i+j]-1]=1;//newelementslist is in Matlab indexing
                         }
@@ -3552 +3552 @@
         /*Free ressources:*/
         xDelete<int>(epart);
-        xDelete<int>(npart);
+        xDelete<int>(npart);        
         xDelete<int>(index);
 }
 /*}}}*/
 void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/
-
+        
         int elementswidth                                                       = this->GetElementsWidth();//just 2D mesh, tria elements
    int numberofvertices                                         = this->vertices->NumberOfVertices();
    IssmDouble weight                                            = 0.;
-        IssmDouble*     tauxx                                                   = NULL;
-        IssmDouble*     tauyy                                                   = NULL;
-        IssmDouble*     tauxy                                                   = NULL;
+        IssmDouble*     tauxx                                                   = NULL; 
+        IssmDouble*     tauyy                                                   = NULL; 
+        IssmDouble*     tauxy                                                   = NULL; 
    IssmDouble* totalweight                              = NULL;
         IssmDouble* deviatoricstressxx          = xNew<IssmDouble>(elementswidth);
@@ -3573 +3573 @@
    Vector<IssmDouble>* vectauxy                 = new Vector<IssmDouble>(numberofvertices);
    Vector<IssmDouble>* vectotalweight   = new Vector<IssmDouble>(numberofvertices);
-
+        
         /*Update the Deviatoric Stress tensor over the elements*/
         this->DeviatoricStressx();
-
+        
    /*Calculate the Smoothed Deviatoric Stress tensor*/
         for(int i=0;i<this->elements->Size();i++){
@@ -3621 +3621 @@
         /*Divide for the total weight*/
         for(int i=0;i<numberofvertices;i++){
-                _assert_(totalweight[i]>0);
+                _assert_(totalweight[i]>0);     
                 tauxx[i] = tauxx[i]/totalweight[i];
                 tauyy[i] = tauyy[i]/totalweight[i];
@@ -3646 +3646 @@
 void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/
 
-        /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor.
+        /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. 
          * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/
 
@@ -3666 +3666 @@
         /*Get smoothed deviatoric stress tensor*/
         this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy);
-
+        
         /*Integrate the error over elements*/
    for(int i=0;i<this->elements->Size();i++){
@@ -3674 +3674 @@
       element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum);
       element->GetVerticesSidList(elem_vertices);
-
+                
                 /*Integrate*/
                 element->GetVerticesCoordinates(&xyz_list);
@@ -3689 +3689 @@
                                 ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n];
                         }
-                        error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2
-                }
-                /*Set the error in the global vector*/
+                        error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 
+                }
+                /*Set the error in the global vector*/  
       sid=element->Sid();
                 error = sqrt(error);//sqrt(e^2)
@@ -3705 +3705 @@
    /*Serialize and set output*/
    (*pelementerror)=velementerror->ToMPISerial();
-
+        
         /*Cleanup*/
         xDelete<IssmDouble>(smoothedtauxx);
@@ -3749 +3749 @@
       Tria* triaelement = xDynamicCast<Tria*>(element);
       weight            = triaelement->GetArea();//the tria area is a choice for the weight
-
+      
                 /*dH/dx*/
       vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL);
@@ -3817 +3817 @@
    /*Get smoothed deviatoric stress tensor*/
    this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy);
-
+   
         /*Integrate the error over elements*/
    for(int i=0;i<this->elements->Size();i++){
@@ -3903 +3903 @@
         IssmDouble* x                                   = NULL;
         IssmDouble* y                                   = NULL;
-        IssmDouble* z                                   = NULL;
+        IssmDouble* z                                   = NULL; 
         IssmDouble* xyz_list                    = NULL;
         IssmDouble x1,y1,x2,y2,x3,y3;
@@ -3912 +3912 @@
       //element->GetVerticesSidList(elem_vertices);
       int sid = element->Sid();
-                element->GetVerticesCoordinates(&xyz_list);
+                element->GetVerticesCoordinates(&xyz_list); 
                 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
                 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
@@ -3945 +3945 @@
                 _error_("level set type not implemented yet!");
         }
-
+        
         /*Outputs*/
         IssmDouble* zerolevelset_points                 = NULL;
         int npoints                                                                             = 0;
-
+        
         /*Intermediaries*/
         int elementswidth                       = this->GetElementsWidth();
@@ -3962 +3962 @@
         int count,sid;
         IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
-
+        
         /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/
    for(int i=0;i<this->elements->Size();i++){
@@ -3969 +3969 @@
                 element->GetVerticesSidList(elem_vertices);
                 sid= element->Sid();
-                element->GetVerticesCoordinates(&xyz_list);
+                element->GetVerticesCoordinates(&xyz_list); 
                 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
                 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
                 x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1];
                 xc      = NAN;
-                yc      = NAN;
+                yc      = NAN;  
         Tria* tria      = xDynamicCast<Tria*>(element);
                 if(tria->IsIceInElement()){/*verify if there is ice in the element*/
-                        if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||
+                        if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. ||  
                                 abs(levelset[0]*levelset[1])<DBL_EPSILON || abs(levelset[0]*levelset[2])<DBL_EPSILON) {
                                 xc=(x1+x2+x3)/3.;
@@ -4007 +4007 @@
                 }
         }
-
+        
         /*Assign outputs*/
         numberofpoints                          = npoints;
@@ -4047 +4047 @@
         responses_pointer=d_responses;
 
-        //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled.
+        //watch out, we have more d_numresponses than numresponsedescriptors, because the responses have been expanded if they were scaled. 
         //because we don't know the d_responses descriptors (the scaled ones) we can't key off them, so we will key off the responses_descriptors: */
 
@@ -4140 +4140 @@
 
         int         ns,nsmax;
-
+        
         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
         ns = elements->Size();
-
+        
         /*Figure out max of ns: */
         ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
@@ -4162 +4162 @@
                 }
         }
-
+        
         /*One last time: */
         pUp->Assemble();
@@ -4181 +4181 @@
 
         int         ns,nsmax;
-
+        
         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
         ns = elements->Size();
-
-        /*First, figure out the surface area of Earth: */
+        
+        /*First, figure out the surface area of Earth: */ 
         for(int i=0;i<ns;i++){
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
@@ -4209 +4209 @@
                 }
         }
-
+        
         /*One last time: */
         pUp->Assemble();
@@ -4226 +4226 @@
 #endif
 #ifdef _HAVE_SEALEVELRISE_
-void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/
+void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
 
         /*serialized vectors:*/
@@ -4262 +4262 @@
                 if(i<ns){
 
-                        if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%  ");
+                        if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "              convolution progress: " << (double)i/(double)ns*100 << "%  ");
 
                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
+                        element->SealevelriseEustatic(pRSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
                         eustatic_cpu+=eustatic_cpu_e;
                 }
-                if(i%100==0)pSgi->Assemble();
+                if(i%loop==0)pRSLgi->Assemble();
         }
         if(VerboseConvergence())_printf0_("\n");
-
+                
         /*One last time: */
-        pSgi->Assemble();
+        pRSLgi->Assemble();
 
         /*Sum all eustatic components from all cpus:*/
@@ -4285 +4285 @@
 }
 /*}}}*/
-void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution){/*{{{*/
+void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
 
         /*serialized vectors:*/
-        IssmDouble* Sg_old=NULL;
-
+        IssmDouble* RSLg_old=NULL;
+        
         IssmDouble  eartharea=0;
         IssmDouble  eartharea_cpu=0;
 
         int         ns,nsmax;
-
+        
         /*Serialize vectors from previous iteration:*/
-        Sg_old=pSg_old->ToMPISerial();
+        RSLg_old=pRSLg_old->ToMPISerial();
 
         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
@@ -4306 +4306 @@
                 eartharea_cpu += element->GetAreaSpherical();
         }
-
+        
         ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
@@ -4317 +4317 @@
         for(int i=0;i<nsmax;i++){
                 if(i<ns){
-                        if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%   ");
+                        if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "              convolution progress: " << (double)i/(double)ns*100 << "%   ");
                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea);
-                }
-                if(i%100==0)pSgo->Assemble();
-        }
-        if(verboseconvolution)if(VerboseConvergence())_printf_("\n");
-
+                        element->SealevelriseNonEustatic(pRSLgo,RSLg_old,latitude,longitude,radius,eartharea);
+                }
+                if(i%loop==0)pRSLgo->Assemble();
+        }
+        if(verboseconvolution)if(VerboseConvergence())_printf0_("\n");
+        
         /*Free ressources:*/
-        xDelete<IssmDouble>(Sg_old);
-}
-/*}}}*/
-void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
+        xDelete<IssmDouble>(RSLg_old);
+}
+/*}}}*/
+void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
 
         /*serialized vectors:*/
-        IssmDouble* Sg_old=NULL;
+        IssmDouble* RSLg_old=NULL;
         IssmDouble  eartharea=0;
         IssmDouble  eartharea_cpu=0;
@@ -4341 +4341 @@
 
         /*Serialize vectors from previous iteration:*/
-        Sg_old=pSg_old->ToMPISerial();
+        RSLg_old=pRSLg_old->ToMPISerial();
 
         /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
@@ -4355 +4355 @@
         for(int i=0;i<elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea);
+                element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,eartharea);
                 moi_list_cpu[0] += moi_list[0];
                 moi_list_cpu[1] += moi_list[1];
@@ -4399 +4399 @@
                                                 (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi)));
 
-                pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
+                pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
         }
 
         /*Assemble mesh velocity*/
-        pSgo_rot->Assemble();
+        pRSLgo_rot->Assemble();
 
         /*Assign output pointers:*/
-        *pIxz=moi_list[0];
-        *pIyz=moi_list[1];
-        *pIzz=moi_list[2];
+        if(pIxz)*pIxz=moi_list[0];
+        if(pIyz)*pIyz=moi_list[1];
+        if(pIzz)*pIzz=moi_list[2];
 
         /*Free ressources:*/
-        xDelete<IssmDouble>(Sg_old);
-
-}
-/*}}}*/
-void FemModel::SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
+        xDelete<IssmDouble>(RSLg_old);
+        
+        
+}
+/*}}}*/
+void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
 
         /*serialized vectors:*/
-        IssmDouble* Sg=NULL;
-
+        IssmDouble* RSLg=NULL;
+        
         IssmDouble  eartharea=0;
         IssmDouble  eartharea_cpu=0;
 
         int         ns,nsmax;
-
+        
         /*Serialize vectors from previous iteration:*/
-        Sg=pSg->ToMPISerial();
+        RSLg=pRSLg->ToMPISerial();
 
         /*Go through elements, and add contribution from each element to the deflection vector wg:*/
         ns = elements->Size();
-
+        
         /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
         for(int i=0;i<ns;i++){
@@ -4446 +4447 @@
         for(int i=0;i<nsmax;i++){
                 if(i<ns){
+                        if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "              convolution progress: " << (double)i/(double)ns*100 << "%  ");
                         Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea);
-                }
-                if(i%100==0){
+                        element->SealevelriseGeodetic(pUp,pNorth,pEast,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
+                }
+                if(i%loop==0){
                         pUp->Assemble();
-                        pNorth->Assemble();
-                        pEast->Assemble();
-                }
-        }
-
+                        if (horiz){
+                                pNorth->Assemble();
+                                pEast->Assemble();
+                        }
+                }
+        }
+        
         /*One last time: */
         pUp->Assemble();
-        pNorth->Assemble();
-        pEast->Assemble();
+        if(horiz){
+                pNorth->Assemble();
+                pEast->Assemble();
+        }
+        if(VerboseConvergence())_printf0_("\n");
 
         /*Free ressources:*/
-        xDelete<IssmDouble>(Sg);
-        xDelete<IssmDouble>(latitude);
-        xDelete<IssmDouble>(longitude);
-        xDelete<IssmDouble>(radius);
-        xDelete<IssmDouble>(xx);
-        xDelete<IssmDouble>(yy);
-        xDelete<IssmDouble>(zz);
-}
-/*}}}*/
-IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* Sg) { /*{{{*/
-
-        IssmDouble* Sg_serial=NULL;
+        xDelete<IssmDouble>(RSLg);
+}
+/*}}}*/
+IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* RSLg) { /*{{{*/
+
+        IssmDouble* RSLg_serial=NULL;
         IssmDouble  oceanvalue,oceanvalue_cpu;
         IssmDouble  oceanarea,oceanarea_cpu;
 
         /*Serialize vectors from previous iteration:*/
-        Sg_serial=Sg->ToMPISerial();
+        RSLg_serial=RSLg->ToMPISerial();
 
         /*Initialize:*/
@@ -4488 +4489 @@
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
                 oceanarea_cpu += element->OceanArea();
-                oceanvalue_cpu += element->OceanAverage(Sg_serial);
+                oceanvalue_cpu += element->OceanAverage(RSLg_serial);
         }
         ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-
+        
         ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
         ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
         /*Free ressources:*/
-        xDelete<IssmDouble>(Sg_serial);
-
+        xDelete<IssmDouble>(RSLg_serial);
+        
         return oceanvalue/oceanarea;
 }
@@ -4514 +4515 @@
         int*                eplzigzag_counter = NULL;
         int                 eplflip_lock;
-
+        
         HydrologyDCEfficientAnalysis* effanalysis =  new HydrologyDCEfficientAnalysis();
         HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
@@ -4521 +4522 @@
         mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
         recurence=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum));
-        this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum);
-        this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum);
+        this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); 
+        this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); 
         GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum);
-
+        
         for (int i=0;i<elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
@@ -4542 +4543 @@
         /*Assemble and serialize*/
         mask->Assemble();
-        serial_mask=mask->ToMPISerial();
-
+        serial_mask=mask->ToMPISerial();        
+        
         xDelete<int>(eplzigzag_counter);
         xDelete<IssmDouble>(serial_rec);
@@ -4585 +4586 @@
         int sum_counter;
         ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
         counter=sum_counter;
         *pEplcount = counter;
         if(VerboseSolution()) _printf0_("   Number of active nodes in EPL layer: "<< counter <<"\n");
-
-        /*Update dof indexings*/
-        this->UpdateConstraintsx();
-
-}
-/*}}}*/
-void FemModel::HydrologyIDSupdateDomainx(IssmDouble* pIDScount){ /*{{{*/
-
-        bool                isthermal;
-        Vector<IssmDouble>* mask                                = NULL;
-        Vector<IssmDouble>* active                              = NULL;
-        IssmDouble*         serial_mask = NULL;
-        IssmDouble*         serial_active       = NULL;
-
-        HydrologyDCInefficientAnalysis* inefanalysis =  new HydrologyDCInefficientAnalysis();
-        parameters->FindParam(&isthermal,TransientIsthermalEnum);
-
-        /*When solving a thermal model we update the thawed nodes*/
-        if(isthermal){
-                /*Step 1: update mask, the mask correspond to thawed nodes (that have a meltingrate)*/
-                mask=new Vector<IssmDouble>(this->nodes->NumberOfNodes(HydrologyDCInefficientAnalysisEnum));
-
-                for (int i=0;i<elements->Size();i++){
-                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        inefanalysis->HydrologyIDSGetMask(mask,element);
-                }
-                /*Assemble and serialize*/
-                mask->Assemble();
-                serial_mask=mask->ToMPISerial();
-                delete mask;
-        }
-        /*for other cases we just grab the mask from the initialisation value*/
-        else{
-                GetVectorFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
-        }
-        /*Update Mask and elementize*/
-        InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskThawedNodeEnum,NodeSIdEnum);
-        xDelete<IssmDouble>(serial_mask);
-        inefanalysis->ElementizeIdsMask(this);
-
-        /*get node mask coherent with element mask*/
-        active=new Vector<IssmDouble>(nodes->NumberOfNodes(HydrologyDCInefficientAnalysisEnum));
-        for (int i=0;i<elements->Size();i++){
-                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                inefanalysis->HydrologyIdsGetActive(active,element);
-        }
-
-        /*Assemble and serialize*/
-        active->Assemble();
-        serial_active=active->ToMPISerial();
-        delete active;
-
-        /*Update node activation accordingly*/
-        int counter =0;
-        for (int i=0;i<nodes->Size();i++){
-                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
-                if(node->InAnalysis(HydrologyDCInefficientAnalysisEnum)){
-                        if(serial_active[node->Sid()]==1.){
-                                node->Activate();
-                                if(!node->IsClone()) counter++;
-                        }
-                        else{
-                                node->Deactivate();
-                        }
-                }
-        }
-
-        xDelete<IssmDouble>(serial_active);
-        delete inefanalysis;
-        int sum_counter;
-        ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
-        counter=sum_counter;
-        *pIDScount = counter;
-        if(VerboseSolution()) _printf0_("   Number of active nodes in IDS layer: "<< counter <<"\n");
 
         /*Update dof indexings*/
@@ -4706 +4632 @@
         int sum_counter;
         ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm());                
         counter=sum_counter;
         *pL2count = counter;
@@ -4791 +4717 @@
 }
 /*}}}*/
-#ifdef _HAVE_JAVASCRIPT_
+#ifdef _HAVE_JAVASCRIPT_ 
 FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/
         /*configuration: */
@@ -4806 +4732 @@
         /*From command line arguments, retrieve different filenames needed to create the FemModel: */
         solution_type=StringToEnumx(solution);
-
+        
         /*Create femmodel from input files: */
         profiler->Start(MPROCESSOR);
         this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL);
         profiler->Stop(MPROCESSOR);
-
+        
         /*Save communicator in the parameters dataset: */
         this->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(incomm,FemModelCommEnum));
@@ -4826 +4752 @@
         size_t* poutputbuffersize;
 
-
+        
         /*Before we delete the profiler, report statistics for this run: */
         profiler->Stop(TOTAL);  //final tagging
@@ -4839 +4765 @@
                                 );
         _printf0_("\n");
-
+        
         /*Before we close the output file, recover the buffer and size:*/
         outputbufferparam = xDynamicCast<GenericParam<char**>*>(this->parameters->FindParamObject(OutputBufferPointerEnum));
@@ -4879 +4805 @@
 
         /*Open output file once for all and add output file descriptor to parameters*/
-        output_fid=open_memstream(&outputbuffer,&outputsize);
+        output_fid=open_memstream(&outputbuffer,&outputsize); 
         if(output_fid==NULL)_error_("could not initialize output stream");
         this->parameters->SetParam(output_fid,OutputFilePointerEnum);
@@ -4887 +4813 @@
 }/*}}}*/
 #endif
+
 
 #if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_)
@@ -4920 +4847 @@
                 /*Initialize hmaxvertices with NAN*/
                 hmaxvertices_serial=xNew<IssmDouble>(numberofvertices);
-                for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN;
+                for(int i=0;i<numberofvertices;i++) hmaxvertices_serial[i]=NAN; 
                 /*Fill hmaxvertices*/
                 if(this->amrbamg->groundingline_distance>0)             this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskGroundediceLevelsetEnum);
@@ -4927 +4854 @@
                 if(this->amrbamg->deviatoricerror_threshold>0)  this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum);
         }
-
+        
         if(my_rank==0){
                 this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist);
@@ -4947 +4874 @@
                 newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
         }
-        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
+        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());        
 
         /*Assign output pointers*/
@@ -4983 +4910 @@
         /*Create bamg data structures for bamg*/
         this->amrbamg = new AmrBamg();
-
+        
         /*Get amr parameters*/
         this->parameters->FindParam(&hmin,AmrHminEnum);
@@ -5007 +4934 @@
 
         /*Re-create original mesh and put it in bamg structure (only cpu 0)*/
-        if(my_rank==0){
+        if(my_rank==0){ 
                 this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements);
         }
@@ -5021 +4948 @@
 
         if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
-
+        
         /*Intermediaries*/
         int numberofvertices                     = this->vertices->NumberOfVertices();
@@ -5028 +4955 @@
 
         switch(levelset_type){
-                case MaskGroundediceLevelsetEnum:
+                case MaskGroundediceLevelsetEnum: 
                         threshold       = this->amrbamg->groundingline_distance;
                         resolution      = this->amrbamg->groundingline_resolution;
@@ -5042 +4969 @@
         this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type);
         if(!verticedistance) _error_("verticedistance is NULL!\n");
-
+        
         /*Fill hmaxVertices*/
         for(int i=0;i<numberofvertices;i++){
@@ -5056 +4983 @@
 /*}}}*/
 void FemModel::GethmaxVerticesFromEstimators(IssmDouble* hmaxvertices,int errorestimator_type){/*{{{*/
-
+   
         if(!hmaxvertices) _error_("hmaxvertices is NULL!\n");
 
@@ -5064 +4991 @@
         int numberofvertices                                            = this->vertices->NumberOfVertices();
         IssmDouble* maxlength                                   = xNew<IssmDouble>(numberofelements);
-        IssmDouble* error_vertices                              = xNewZeroInit<IssmDouble>(numberofvertices);
+        IssmDouble* error_vertices                              = xNewZeroInit<IssmDouble>(numberofvertices);   
         IssmDouble* error_elements                              = NULL;
         IssmDouble* x                                                           = NULL;
@@ -5077 +5004 @@
         /*Fill variables*/
         switch(errorestimator_type){
-                case ThicknessErrorEstimatorEnum:
+                case ThicknessErrorEstimatorEnum: 
                         threshold               = this->amrbamg->thicknesserror_threshold;
                         groupthreshold  = this->amrbamg->thicknesserror_groupthreshold;
@@ -5102 +5029 @@
         case ThicknessErrorEstimatorEnum:                       this->amrbamg->thicknesserror_maximum   = maxerror;break;
         case DeviatoricStressErrorEstimatorEnum:        this->amrbamg->deviatoricerror_maximum = maxerror;break;
-        }
+        }       
         }
 
         /*Get mesh*/
         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index);
-
+        
         /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/
         for(int i=0;i<numberofelements;i++){
@@ -5121 +5048 @@
                 error_vertices[v2]+=error_elements[i];
                 error_vertices[v3]+=error_elements[i];
-        }
+        }       
 
         /*Fill hmaxvertices with the criteria*/
@@ -5135 +5062 @@
                         }
                 }
-                /*Now, fill the hmaxvertices if requested*/
+                /*Now, fill the hmaxvertices if requested*/       
                 if(refine){
                         for(int j=0;j<elementswidth;j++){
@@ -5165 +5092 @@
         /*Output*/
         IssmDouble* verticedistance;
-
+        
         /*Intermediaries*/
    int numberofvertices       = this->vertices->NumberOfVertices();
@@ -5177 +5104 @@
         /*Get vertices coordinates*/
         VertexCoordinatesx(&x,&y,&z,this->vertices,false) ;
-
-        /*Get points which level set is zero (center of elements with zero level set)*/
+        
+        /*Get points which level set is zero (center of elements with zero level set)*/ 
         this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
 
@@ -5187 +5114 @@
                 for(int j=0;j<numberofpoints;j++){
                         distance=sqrt((x[i]-levelset_points[2*j])*(x[i]-levelset_points[2*j])+(y[i]-levelset_points[2*j+1])*(y[i]-levelset_points[2*j+1]));
-                        verticedistance[i]=min(distance,verticedistance[i]);
-                }
-        }
+                        verticedistance[i]=min(distance,verticedistance[i]);            
+                }
+        }       
 
         /*Assign the pointer*/
@@ -5223 +5150 @@
         if(this->amr->groundingline_distance>0)         this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskGroundediceLevelsetEnum);
    if(this->amr->icefront_distance>0)                           this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum);
-   if(this->amr->thicknesserror_threshold>0)            this->ThicknessZZErrorEstimator(&thicknesserror);
-        if(this->amr->deviatoricerror_threshold>0)      this->ZZErrorEstimator(&deviatoricerror);
-
+   if(this->amr->thicknesserror_threshold>0)            this->ThicknessZZErrorEstimator(&thicknesserror);       
+        if(this->amr->deviatoricerror_threshold>0)      this->ZZErrorEstimator(&deviatoricerror);       
+        
         if(my_rank==0){
                 this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror,
-                                                                                                &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist);
+                                                                                                &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); 
                 newz=xNewZeroInit<IssmDouble>(newnumberofvertices);
                 if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz.");
@@ -5242 +5169 @@
                 newelementslist=xNew<int>(newnumberofelements*this->GetElementsWidth());
         }
-        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());
-
-        /*Assign the pointers*/
+        ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+        ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+        ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); 
+        ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm());        
+
+        /*Assign the pointers*/ 
         (*pnewelementslist)     = newelementslist; //Matlab indexing
         (*pnewx)                                        = newx;
@@ -5263 +5190 @@
 /*}}}*/
 void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/
-
+        
         /*Define variables*/
         int my_rank                                                                             = IssmComm::GetRank();
@@ -5272 +5199 @@
         IssmDouble* z                                                                   = NULL;
         int* elements                                                                   = NULL;
-
+        
         /*Initialize field as NULL for now*/
         this->amr = NULL;
@@ -5280 +5207 @@
         this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements);
 
-        /*Create initial mesh (coarse mesh) in neopz data structure*/
+        /*Create initial mesh (coarse mesh) in neopz data structure*/ 
         /*Just CPU #0 should keep AMR object*/
    /*Initialize refinement pattern*/
         this->SetRefPatterns();
         this->amr = new AdaptiveMeshRefinement();
-        this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster)
+        this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) 
         /*Get amr parameters*/
         this->parameters->FindParam(&this->amr->level_max,AmrLevelMaxEnum);
@@ -5298 +5225 @@
         this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum);
         this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum);
-        if(my_rank==0){
+        if(my_rank==0){ 
                 this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements);
         }
@@ -5319 +5246 @@
         /*Output*/
         IssmDouble* elementdistance;
-
+        
         /*Intermediaries*/
    int numberofelements                                                 = this->elements->NumberOfElements();
@@ -5328 +5255 @@
         IssmDouble xc,yc,x1,y1,x2,y2,x3,y3;
         int numberofpoints;
-
-        /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/
+        
+        /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/      
         this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type);
 
@@ -5335 +5262 @@
       Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(i));
       int sid = element->Sid();
-                element->GetVerticesCoordinates(&xyz_list);
+                element->GetVerticesCoordinates(&xyz_list); 
                 x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1];
                 x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1];
@@ -5345 +5272 @@
                 for(int j=0;j<numberofpoints;j++){
                         distance =sqrt((xc-levelset_points[2*j])*(xc-levelset_points[2*j])+(yc-levelset_points[2*j+1])*(yc-levelset_points[2*j+1]));
-                        mindistance=min(distance,mindistance);
+                        mindistance=min(distance,mindistance);          
                 }
                 velementdistance->SetValue(sid,mindistance,INS_VAL);
                 xDelete<IssmDouble>(xyz_list);
-        }
+        }       
 
    /*Assemble*/

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -135 +135 @@
                 #endif
                 #ifdef _HAVE_ESA_
-                void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy);
-                void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
+                void EsaGeodetic2D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY, IssmDouble* xx, IssmDouble* yy); 
+                void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); 
                 #endif
                 #ifdef _HAVE_SEALEVELRISE_
-                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
-                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution);
-                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
-                void SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
+                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
+                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
+                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
+                void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); 
                 IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
                 #endif

issm/trunk-jpl/src/c/classes/Materials/Matpar.h

@@ -118 +118 @@
                 void       ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");};
                 void       ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,IssmDouble eps_eff){_error_("not supported");};
+                void       ViscosityBFS(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
+                void       ViscosityBHO(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
+                void       ViscosityBSSA(IssmDouble* pmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
                 /*}}}*/
                 /*Numerics: {{{*/

issm/trunk-jpl/src/c/cores/cores.h

@@ -54 +54 @@
 void damage_core(FemModel* femmodel);
 void sealevelrise_core(FemModel* femmodel);
+void geodetic_core(FemModel* femmodel);
+void steric_core(FemModel* femmodel);
 Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel);
-Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic);
+Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic);
+void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg);
+void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia,Vector<IssmDouble>** pN_gia,FemModel*  femmodel,Vector<IssmDouble>* RSLg);
+void sealevelrise_diagnostics(FemModel* femmodel,Vector<IssmDouble>* RSLg);
 IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);
 
@@ -68 +73 @@
 void TransferForcing(FemModel* femmodel,int forcingenum);
 void TransferSealevel(FemModel* femmodel,int forcingenum);
+void EarthMassTransport(FemModel* femmodel);
+void slrconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
 
 //solution configuration

issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp

@@ -10 +10 @@
 #include "../solutionsequences/solutionsequences.h"
 
+/*cores:*/
 void sealevelrise_core(FemModel* femmodel){ /*{{{*/
 
-        Vector<IssmDouble> *Sg    = NULL;
-        Vector<IssmDouble> *Sg_absolute  = NULL; 
-        Vector<IssmDouble> *Sg_eustatic  = NULL; 
-        Vector<IssmDouble> *slr_initial  = NULL; 
+
+        /*Parameters, variables:*/
+        bool save_results;
+        bool isslr=0;
+        int solution_type;
+                
+        /*Retrieve parameters:*/
+        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
+        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+        
+        /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
+        if(solution_type==SealevelriseSolutionEnum)isslr=1;
+
+        /*Should we be here?:*/
+        if(!isslr)return;
+
+        /*Verbose: */
+        if(VerboseSolution()) _printf0_("   computing sea level rise\n");
+
+        /*set SLR configuration: */
+        femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
+
+        /*Run geodetic:*/
+        geodetic_core(femmodel);
+
+        /*Run steric core for sure:*/
+        steric_core(femmodel);
+        
+        /*Save results: */
+        if(save_results){
+                int        numoutputs        = 0;
+                char     **requested_outputs = NULL;
+
+                if(VerboseSolution()) _printf0_("   saving results\n");
+                femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
+                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
+                if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
+        }
+
+        /*requested dependents: */
+        if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
+}
+/*}}}*/
+void geodetic_core(FemModel* femmodel){ /*{{{*/
+
+
+        /*variables:*/
+        Vector<IssmDouble> *RSLg    = NULL;
+        Vector<IssmDouble> *RSLg_rate    = NULL;
+        Vector<IssmDouble> *RSLg_eustatic  = NULL; 
+        Vector<IssmDouble> *U_esa  = NULL; 
+        Vector<IssmDouble> *U_esa_rate  = NULL; 
+        Vector<IssmDouble> *N_esa  = NULL; 
+        Vector<IssmDouble> *N_esa_rate  = NULL; 
+        Vector<IssmDouble> *U_north_esa   = NULL; 
+        Vector<IssmDouble> *U_east_esa    = NULL; 
+        Vector<IssmDouble> *N_gia= NULL; 
+        Vector<IssmDouble> *U_gia= NULL; 
+        Vector<IssmDouble> *N_gia_rate= NULL; 
+        Vector<IssmDouble> *U_gia_rate= NULL; 
+
+        /*parameters:*/
+        bool iscoupler;
+        int  solution_type;
+        int  configuration_type;
+        int  modelid,earthid;
+        bool istransientmasstransport;
+        int  frequency,count;
+        int  horiz;
+        int  geodetic=0;
+        IssmDouble          dt;
+
+        /*Should we even be here?:*/
+        femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return;
+        
+        /*Verbose: */
+        if(VerboseSolution()) _printf0_("         computing geodetic sea level rise\n");
+
+        /*retrieve more parameters:*/
+        femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
+        femmodel->parameters->FindParam(&frequency,SealevelriseGeodeticRunFrequencyEnum);
+        femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum);
+        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum);
+        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+
+        if(iscoupler){
+                femmodel->parameters->FindParam(&modelid,ModelIdEnum);
+                femmodel->parameters->FindParam(&earthid,EarthIdEnum);
+                femmodel->parameters->FindParam(&istransientmasstransport,TransientIsmasstransportEnum);
+        }
+        else{
+                /* we are here, we are not running in a coupler, so we will indeed compute SLR,
+                 * so make sure we are identified as being the Earth.:*/
+                modelid=1; earthid=1; 
+                /* in addition, if we are running solution_type SealevelriseSolutionEnum, make sure we 
+                 * run, irresepective of the time settings:*/
+                count=frequency;
+        }
+
+        /*If we are running in coupled mode, the Earth model needs to run its own mass transport (if 
+         * not already done by the mass trasnport module. For ice caps, they rely on the transient mass 
+         * transport module exclusively:*/
+        if(iscoupler) if(modelid==earthid) if(!istransientmasstransport) EarthMassTransport(femmodel);
+
+        /*increment counter, or call solution core if count==frequency:*/
+        if (count<frequency){
+                count++; femmodel->parameters->SetParam(count,SealevelriseRunCountEnum); 
+                return;
+        }
+
+        /*call sea-level rise sub cores:*/
+        if(iscoupler){
+                /*transfer cumulated deltathickness forcing from ice caps to earth model: */
+                TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum);
+
+                /*we have accumulated thicknesses, dump them in deltathcikness: */
+                if(modelid==earthid)InputDuplicatex(femmodel,SealevelriseCumDeltathicknessEnum,SealevelriseDeltathicknessEnum);
+        }
+
+        /*run cores:*/
+        if(modelid==earthid){
+
+                /*call eustatic core  (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
+                RSLg_eustatic=sealevelrise_core_eustatic(femmodel); 
+
+                /*call non-eustatic core (ocean loading tems  - 2nd and 5th terms on the RHS of Farrel and Clark) */
+                RSLg=sealevelrise_core_noneustatic(femmodel,RSLg_eustatic); 
+
+                /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */
+                sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg);
+
+                /*compute viscosus (GIA) geodetic signatures:*/
+                sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg);
+
+                /*compute sea-level rise (low-order spherical harmonics coefficients) diagnostics:*/
+                sealevelrise_diagnostics(femmodel,RSLg);
+
+                /*recover N_esa  = U_esa + RSLg:*/
+                N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1);
+                
+                /*transform these values into rates (as we only run this once each frequency turn:*/
+                N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency));
+                U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency));
+                N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency));
+                U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency));
+                RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency));
+                
+                /*get some results into elements:{{{*/
+                InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum); 
+                InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum); 
+                InputUpdateFromVectorx(femmodel,U_gia_rate,SealevelUGiaRateEnum,VertexSIdEnum);
+                InputUpdateFromVectorx(femmodel,N_gia_rate,SealevelNGiaRateEnum,VertexSIdEnum); 
+                InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum); 
+                InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum); 
+                InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum); 
+                InputUpdateFromVectorx(femmodel,U_gia,SealevelUGiaEnum,VertexSIdEnum); 
+                InputUpdateFromVectorx(femmodel,N_gia,SealevelNGiaEnum,VertexSIdEnum); 
+                InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum); 
+                InputUpdateFromVectorx(femmodel,RSLg_eustatic,SealevelRSLEustaticEnum,VertexSIdEnum); 
+
+                if (horiz){
+                        InputUpdateFromVectorx(femmodel,U_north_esa,SealevelUNorthEsaEnum,VertexSIdEnum);       // north motion 
+                        InputUpdateFromVectorx(femmodel,U_east_esa,SealevelUEastEsaEnum,VertexSIdEnum);         // east motion 
+                } /*}}}*/
+        }
+
+
+        if(iscoupler){
+                /*transfer sea level back to ice caps:*/
+                TransferSealevel(femmodel,SealevelNEsaRateEnum);
+                TransferSealevel(femmodel,SealevelNGiaRateEnum);
+                TransferSealevel(femmodel,SealevelUEsaRateEnum);
+                TransferSealevel(femmodel,SealevelUGiaRateEnum);
+
+                //reset cumdeltathickness  to 0: 
+                InputUpdateFromConstantx(femmodel,0.0,SealevelriseCumDeltathicknessEnum);
+        }
+
+        /*reset counter to 1:*/
+        femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter.
+
+        /*free ressources:{{{*/
+        delete RSLg;
+        delete RSLg_rate;
+        delete RSLg_eustatic;
+        delete U_esa;
+        delete U_esa_rate;
+        delete N_esa;
+        delete N_esa_rate;
+
+        if(horiz){
+                delete U_north_esa;
+                delete U_east_esa;
+        }
+        delete N_gia;
+        delete U_gia;
+        delete N_gia_rate;
+        delete U_gia_rate;
+        /*}}}*/
+
+} 
+/*}}}*/
+void steric_core(FemModel* femmodel){ /*{{{*/
+
+        /*variables:*/
+        Vector<IssmDouble> *bedrock  = NULL; 
+        Vector<IssmDouble> *SL  = NULL; 
         Vector<IssmDouble> *steric_rate_g  = NULL; 
-        Vector<IssmDouble> *U_radial  = NULL; 
-        Vector<IssmDouble> *U_north   = NULL; 
-        Vector<IssmDouble> *U_east    = NULL; 
-        bool save_results,isslr,iscoupler;
-        int configuration_type;
-        int solution_type;
-        int        numoutputs        = 0;
-        char     **requested_outputs = NULL;
-        
-        /*additional parameters: */
+        Vector<IssmDouble> *U_esa_rate= NULL;
+        Vector<IssmDouble> *N_esa_rate= NULL;
+        Vector<IssmDouble> *U_gia_rate= NULL;
+        Vector<IssmDouble> *N_gia_rate= NULL;
+                
+        /*parameters: */
+        bool isslr=0;
+        int  solution_type;
+        IssmDouble          dt;
+        int  geodetic=0;
+
+        /*Retrieve parameters:*/
+        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
+        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+        femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); 
+
+        /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
+        if(solution_type==SealevelriseSolutionEnum)isslr=1;
+        
+        /*Should we be here?:*/
+        if(!isslr)return;
+
+        /*Verbose: */
+        if(VerboseSolution()) _printf0_("         computing steric sea level rise\n");
+
+        /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/
+        GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
+        GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum);
+        GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexSIdEnum);
+        if(geodetic){
+                GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum);
+                GetVectorFromInputsx(&U_gia_rate,femmodel,SealevelUGiaRateEnum,VertexSIdEnum);
+                GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum);
+                GetVectorFromInputsx(&N_gia_rate,femmodel,SealevelNGiaRateEnum,VertexSIdEnum);
+        }
+
+        /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate)  * dt + steric_rate * dt*/
+        if(geodetic){
+                SL->AXPY(N_gia_rate,dt);
+                SL->AXPY(N_esa_rate,dt);
+        }
+        SL->AXPY(steric_rate_g,dt);
+
+        /*compute new bedrock position: */
+        if(geodetic){
+                bedrock->AXPY(U_esa_rate,dt);
+                bedrock->AXPY(U_gia_rate,dt);
+        }
+
+        /*update element inputs:*/
+        InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum); 
+        InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum); 
+
+        /*Free ressources:*/    
+        delete bedrock;
+        delete SL;
+        delete steric_rate_g;
+        if(geodetic){
+                delete U_esa_rate;
+                delete U_gia_rate;
+                delete N_esa_rate;
+                delete N_gia_rate;
+        }
+}
+/*}}}*/
+Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/
+ 
+        /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
+
+        Vector<IssmDouble> *RSLgi    = NULL;
+        IssmDouble          RSLgi_oceanaverage   = 0;
+
+        /*parameters: */
+        int  configuration_type;
         int  gsize;
         bool spherical=true;
+        IssmDouble *latitude  = NULL;
+        IssmDouble *longitude = NULL;
+        IssmDouble *radius    = NULL;
+        int         loop;
+
+        /*outputs:*/
+        IssmDouble eustatic;
+
+        /*recover parameters:*/
+        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
+
+        /*first, recover lat,long and radius vectors from vertices: */
+        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); 
+
+        /*Figure out size of g-set deflection vector and allocate solution vector: */
+        gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
+        
+        /*Initialize:*/
+        RSLgi = new Vector<IssmDouble>(gsize);
+
+        /*call the eustatic main module: */
+        femmodel->SealevelriseEustatic(RSLgi,&eustatic, latitude, longitude, radius,loop); //this computes 
+
+        /*we need to average RSLgi over the ocean: RHS term  4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
+        RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi);
+
+        /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the 
+         * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/
+        RSLgi->Shift(-eustatic-RSLgi_oceanaverage);
+        
+        /*save eustatic value for results: */
+        femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic));
+
+
+        /*clean up and return:*/
+        xDelete<IssmDouble>(latitude);
+        xDelete<IssmDouble>(longitude);
+        xDelete<IssmDouble>(radius);
+        return RSLgi;
+}/*}}}*/
+Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic){ /*{{{*/
+
+        /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
+          non eustatic core of the SLR solution */
+
+
+        Vector<IssmDouble> *RSLg    = NULL;
+        Vector<IssmDouble> *RSLg_old    = NULL;
+
+        Vector<IssmDouble> *RSLgo    = NULL; //ocean convolution of the perturbation to gravity potential.
+        Vector<IssmDouble> *RSLgo_rot= NULL; // rotational feedback 
+        IssmDouble          RSLgo_oceanaverage = 0;  //average of RSLgo over the ocean.
+
+        /*parameters: */
+        int count;
+        bool save_results;
+        int  gsize;
+        int  configuration_type;
+        bool spherical=true;
+        bool converged=true;
+        bool rotation=true;
+        bool verboseconvolution=true;
+        int max_nonlinear_iterations;
+        IssmDouble           eps_rel;
+        IssmDouble           eps_abs;
+        IssmDouble          *latitude    = NULL;
+        IssmDouble          *longitude    = NULL;
+        IssmDouble          *radius    = NULL;
+        IssmDouble           eustatic;
+        int                      loop;
+
+        /*Recover some parameters: */
+        femmodel->parameters->FindParam(&max_nonlinear_iterations,SealevelriseMaxiterEnum);
+        femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum);
+        femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum);
+        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
+        
+        /*computational flag: */
+        femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum);
+
+        /*first, recover lat,long and radius vectors from vertices: */
+        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); 
+
+        /*Figure out size of g-set deflection vector and allocate solution vector: */
+        gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
+        
+        /*Initialize:*/
+        RSLg = new Vector<IssmDouble>(gsize);
+        RSLg->Assemble();
+        RSLg_eustatic->Copy(RSLg);  //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic.
+
+        RSLg_old = new Vector<IssmDouble>(gsize);
+        RSLg_old->Assemble();
+
+        count=1;
+        converged=false;
+        
+        /*Start loop: */
+        for(;;){
+
+                //save pointer to old sea level rise
+                delete RSLg_old; RSLg_old=RSLg; 
+
+                /*Initialize solution vector: */
+                RSLg  = new Vector<IssmDouble>(gsize); RSLg->Assemble();
+                RSLgo = new Vector<IssmDouble>(gsize); RSLgo->Assemble();
+
+                /*call the non eustatic module: */
+                femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop);
+        
+                /*assemble solution vector: */
+                RSLgo->Assemble(); 
+
+                if(rotation){
+                        /*call rotational feedback  module: */
+                        RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
+                        femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,NULL,NULL,NULL,latitude,longitude,radius); 
+                        RSLgo_rot->Assemble(); 
+
+                        RSLgo->AXPY(RSLgo_rot,1); 
+                }
+                
+                /*we need to average RSLgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
+                RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo);
+        
+                /*RSLg is the sum of the eustatic term, and the ocean terms: */
+                RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1); 
+                RSLg->Shift(-RSLgo_oceanaverage);
+
+                /*convergence criterion:*/
+                slrconvergence(&converged,RSLg,RSLg_old,eps_rel,eps_abs);
+
+                /*free ressources: */
+                delete RSLgo;
+
+                /*Increase count: */
+                count++;
+                if(converged==true){
+                        break;
+                }
+                if(count>=max_nonlinear_iterations){
+                        _printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n"); 
+                        converged=true;
+                        break;
+                }       
+                
+                /*some minor verbosing adjustment:*/
+                if(count>1)verboseconvolution=false;
+                
+        }
+        if(VerboseConvergence()) _printf0_("\n              total number of iterations: " << count-1 << "\n");
+
+        xDelete<IssmDouble>(latitude);
+        xDelete<IssmDouble>(longitude);
+        xDelete<IssmDouble>(radius);
+        delete RSLg_old;
+
+        return RSLg;
+} /*}}}*/
+void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/
+
+        Vector<IssmDouble> *U_esa  = NULL; 
+        Vector<IssmDouble> *U_north_esa   = NULL; 
+        Vector<IssmDouble> *U_east_esa    = NULL; 
+                
+        /*parameters: */
+        int configuration_type;
+        int  gsize;
+        bool spherical=true;
+
         IssmDouble          *latitude   = NULL;
         IssmDouble          *longitude  = NULL;
@@ -35 +479 @@
         IssmDouble          *yy     = NULL;
         IssmDouble          *zz     = NULL;
+        int  loop;
+        int  horiz;
+
+        /*retrieve some parameters:*/
+        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
+        femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum);
+
+        /*find size of vectors:*/
+        gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
+
+        /*intialize vectors:*/
+        U_esa = new Vector<IssmDouble>(gsize);
+        if (horiz){
+                U_north_esa = new Vector<IssmDouble>(gsize);
+                U_east_esa = new Vector<IssmDouble>(gsize);
+        }
+
+        /*retrieve geometric information: */
+        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); 
+        VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); 
+        
+        /*call the elastic main modlule:*/ 
+        femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,latitude,longitude,radius,xx,yy,zz,loop,horiz);
+
+        /*Assign output pointers:*/
+        *pU_esa=U_esa;
+        if(horiz){
+                *pU_east_esa=U_east_esa;
+                *pU_north_esa=U_north_esa;
+        }
+
+        /*Free ressources: */
+        delete longitude; 
+        delete latitude; 
+        delete radius; 
+        delete xx; 
+        delete yy; 
+        delete zz; 
+}
+/*}}}*/
+void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia, Vector<IssmDouble>** pN_gia,FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/
+
+        /*variables:*/
+        Vector<IssmDouble> *U_gia  = NULL; 
+        Vector<IssmDouble> *N_gia  = NULL; 
+        
+        /*parameters:*/
+        int                                     frequency;
         IssmDouble          dt;
 
-        /*Recover some parameters: */
-        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
-        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
-        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
-        femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
-
-        /*first, recover lat,long and radius vectors from vertices: */
-        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); 
-
-        /*recover x,y,z vectors from vertices: */
-        VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); 
-
-        /*Figure out size of g-set deflection vector and allocate solution vector: */
-        gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
-
-        /*several cases here, depending on value of iscoupler and isslr: 
-        solution_type == SealevelriseSolutionEnum)       we are running sea level rise core (no coupler)
-        ( !iscoupler & !isslr)       we are not interested in being here :) 
-        ( !iscoupler & isslr)        we are running in uncoupled mode
-        ( iscoupler & isslr)         we are running in coupled mode, and better be earth
-        ( iscoupler & !isslr)        we are running in coupled mode, and better be an ice cap
-        */
-
-        if(solution_type==SealevelriseSolutionEnum){
-                isslr=1;
-                iscoupler=0;
-        }
-
-        /*early return: */
-        if( !iscoupler & !isslr) return;  //we are not interested in being here :) 
-
-        /*In what follows we assume we are all running slr, either in coupled, or uncoupled mode:*/
-        if(VerboseSolution()) _printf0_("   computing sea level rise\n");
-
-        /*set configuration: */
-        if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
-
-        /*transfer deltathickness forcing from ice caps to earth model: */
-        if(iscoupler) TransferForcing(femmodel,SealevelriseDeltathicknessEnum);
-
-        /*call sea-level rise sub cores:*/
-        if(isslr){
-                Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS.
-        
-                /* The following is for reproducing Farrell&Clark1976 Fig. 1, and aimed for the workshop! 
-                 * md.slr.sealevel is considered as the "initial sea-level" where 1 m slr is distributed uniformly over the ocean 
-                 * similar strategy may work for computin SAL due to "internal mass distribution" associated with ocean dynamics 
-                 * if it breaks the code, it will be reverted. And, we will strategize how we want to accomodate */
-                GetVectorFromInputsx(&slr_initial,femmodel,SealevelEnum,VertexSIdEnum);
-                Sg_eustatic->AXPY(slr_initial,1); 
-
-                Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //ocean loading tems  (2nd and 5th terms on the RHS of Farrel and Clark)
-
-                /*compute other geodetic signatures, such as absolute sea level chagne, components of 3-D crustal motion: */
-                /*Initialize:*/
-                U_radial = new Vector<IssmDouble>(gsize);
-                U_north = new Vector<IssmDouble>(gsize);
-                U_east = new Vector<IssmDouble>(gsize);
-                Sg_absolute = new Vector<IssmDouble>(gsize); 
-                
-                /*call the geodetic main modlule:*/ 
-                femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz); 
-
-                /*Now deal with steric ocean expansion by just shifting Sg by a spatial rate pattern : */
-                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-                GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexSIdEnum);
-                Sg->AXPY(steric_rate_g,dt);
-                
-                /*compute: absolute sea level change = relative sea level change + vertical motion*/
-                Sg->Copy(Sg_absolute); Sg_absolute->AXPY(U_radial,1); 
-                
-                /*get results into elements:*/
-                InputUpdateFromVectorx(femmodel,U_radial,SealevelUmotionEnum,VertexSIdEnum);    // radial displacement 
-                InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum);     // north motion 
-                InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum);              // east motion 
-                InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelAbsoluteEnum,VertexSIdEnum); //absolute sea level
-                InputUpdateFromVectorx(femmodel,Sg,SealevelEnum,VertexSIdEnum); //relative sea level
-                
-                if(save_results){
-                        if(VerboseSolution()) _printf0_("   saving results\n");
-                        femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
-                        femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
-                }
-                        
-                if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
-
-                /*Free ressources:*/    
-                delete Sg;
-                delete Sg_eustatic;
-                delete U_radial;
-                delete U_north;
-                delete U_east;
-                delete Sg_absolute;
-                delete steric_rate_g;
-                if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
-        }
-
-        /*transfer sea-level back to ice caps: */
-        if(iscoupler)TransferSealevel(femmodel,SealevelEnum);
-} 
+        /*retrieve some parameters:*/
+        femmodel->parameters->FindParam(&frequency,SealevelriseGeodeticRunFrequencyEnum);
+        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
+        /*recover GIA rates:*/
+        GetVectorFromInputsx(&U_gia,femmodel,SealevelUGiaRateEnum,VertexSIdEnum);
+        GetVectorFromInputsx(&N_gia,femmodel,SealevelNGiaRateEnum,VertexSIdEnum);
+
+        /*we just loaded rates, that's not what's being asked, scale by time:*/ 
+        U_gia->Scale(frequency*dt);
+        N_gia->Scale(frequency*dt);
+
+        /*Assign output pointers:*/
+        *pU_gia=U_gia;
+        *pN_gia=N_gia;
+
+}
 /*}}}*/
+void sealevelrise_diagnostics(FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/
+
+        /*compute spherical harmonics deg 1 and deg 2 coefficeints:*/
+
+}
+/*}}}*/
+
+/*support routines:*/
 void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/
 
@@ -211 +628 @@
                  *First, build the global delta thickness vector in the earth model: */
                 nv=femmodel->vertices->NumberOfVertices();
-                forcingglobal= new Vector<IssmDouble>(nv);
+                GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
 
                 /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/
@@ -232 +649 @@
                                 /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular 
                                  * ice cap: */
-                                forcingglobal->SetValues(M,index,forcingfromcap,INS_VAL);
+                                forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL);
                                 xDelete<int>(index);
                         }
@@ -372 +789 @@
 
 } /*}}}*/
+void EarthMassTransport(FemModel* femmodel){ /*{{{*/
+
+        IssmDouble time,dt;
+        Vector<IssmDouble> *oldthickness    = NULL;
+        Vector<IssmDouble> *newthickness    = NULL;
+        Vector<IssmDouble> *deltathickness    = NULL;
+        Vector<IssmDouble> *cumdeltathickness    = NULL;
+        int nv;
+        
+        if(VerboseSolution()) _printf0_("              computing earth mass transport\n");
+
+        /*This mass transport module for the Earth is because we might have thickness variations as spcs 
+         * specified in the md.slr class, outside of what we will get from the icecaps. That's why we get t
+         * the thickness variations from SealevelriseSpcthicknessEnum.*/
+
+        /*No mass transport module was called, so we are just going to retrieve the geometry thickness 
+         * at this time step, at prior time step, and plug the difference as deltathickness: */
+        femmodel->parameters->FindParam(&time,TimeEnum);
+        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+        nv=femmodel->vertices->NumberOfVertices();
+
+        GetVectorFromInputsx(&newthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum);
+        GetVectorFromInputsx(&oldthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum,time-dt);
+
+        /*compute deltathickness: */
+        deltathickness = new Vector<IssmDouble>(nv); 
+        newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1); 
+
+        /*plug into elements:*/
+        InputUpdateFromVectorx(femmodel,deltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum);
+
+        /*add to cumulated delta thickness: */
+        GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
+        cumdeltathickness->AXPY(deltathickness,1); 
+        InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
+
+        /*free ressources:*/
+        delete oldthickness;
+        delete newthickness;
+        delete deltathickness;
+        delete cumdeltathickness;
+
+} /*}}}*/
+void slrconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
+        
+        bool converged=true;
+        IssmDouble ndS,nS; 
+        Vector<IssmDouble> *dRSLg    = NULL;
+
+        //compute norm(du) and norm(u) if requested
+        dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
+        ndS=dRSLg->Norm(NORM_TWO); 
+        
+        if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");
+        
+        if(!xIsNan<IssmDouble>(eps_rel)){
+                nS=RSLg_old->Norm(NORM_TWO);
+                if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
+        }
+
+
+        //clean up
+        delete dRSLg;
+
+        //print
+        if(!xIsNan<IssmDouble>(eps_rel)){
+                if((ndS/nS)<eps_rel){
+                        if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " < " << eps_rel*100 << " %\n");
+                }
+                else{ 
+                        if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " > " << eps_rel*100 << " %\n");
+                        converged=false;
+                }
+        }
+        if(!xIsNan<IssmDouble>(eps_abs)){
+                if(ndS<eps_abs){
+                        if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)" << ndS << " < " << eps_abs << " \n");
+                }
+                else{ 
+                        if(VerboseConvergence()) _printf0_(setw(50) << left << "              convergence criterion: norm(dS)" << ndS << " > " << eps_abs << " \n");
+                        converged=false;
+                }
+        }
+
+        /*assign output*/
+        *pconverged=converged;
+
+} /*}}}*/

issm/trunk-jpl/src/c/cores/transient_core.cpp

@@ -358 +358 @@
 
                 /*Sea level rise: */
-                if(isslr || iscoupler) sealevelrise_core(femmodel);
+                if(isslr) sealevelrise_core(femmodel);
 
                 /*unload results*/

issm/trunk-jpl/src/m/classes/maskpsl.m

@@ -10 +10 @@
                 ocean_levelset = NaN;
                 land_levelset = NaN;
+                glacier_levelset = NaN;
         end
         methods (Static)
@@ -68 +69 @@
                         fielddisplay(self,'ocean_levelset','is the vertex on the ocean ? yes if = 1, no if = 0');
                         fielddisplay(self,'land_levelset','is the vertex on the land ? yes if = 1, no if = 0');
+                        fielddisplay(self,'glacier_levelset','is the vertex on a glacier ? yes if = 1, no if = 0');
                 end % }}}
                 function marshall(self,prefix,md,fid) % {{{
+                        WriteData(fid,prefix,'name','md.mask.type','data',class(md.mask),'format','String');
                         WriteData(fid,prefix,'object',self,'class','mask','fieldname','groundedice_levelset','format','DoubleMat','mattype',1);
                         WriteData(fid,prefix,'object',self,'class','mask','fieldname','ice_levelset','format','DoubleMat','mattype',1);

issm/trunk-jpl/src/m/classes/mesh3dsurface.m

@@ -103 +103 @@
                 end % }}}
                 function disp(obj) % {{{
-                        disp(sprintf('   2D tria Mesh (horizontal):')); 
+                        disp(sprintf('   3D tria Mesh (surface):')); 
 
                         disp(sprintf('\n      Elements and vertices:'));

issm/trunk-jpl/src/m/classes/model.js

@@ -83 +83 @@
                         this.levelset         = new levelset();
                         this.calving          = new calving();
+                        this.gia              = new gia();
                         this.love             = new fourierlove();
-                        this.gia              = new giaivins();
-
                         this.esa              = new esa();
                         this.autodiff         = new autodiff();
                         this.inversion        = new inversion();
                         this.qmu              = new qmu();
-                        this.amr              = new amr();
-
+                        this.amr                                         = new amr();
+                        this.radaroverlay     = new radaroverlay();
                         this.results          = {};
                         this.outputdefinition = new outputdefinition();
@@ -764 +763 @@
                 this.inversion        = 0;
                 this.qmu              = 0;
-                this.amr              = 0;
+                this.amr                                         = 0;
 
                 this.results          = 0;

issm/trunk-jpl/src/m/classes/model.m

@@ -38 +38 @@
                 steadystate      = 0;
                 transient        = 0;
-                levelset                          = 0;
+                levelse                  = 0;
                 calving          = 0;
-                love                         = 0;
-                gia                               = 0;
+                gia                              = 0;
                 esa              = 0;
-
+                love                     = 0;
                 autodiff         = 0;
                 inversion        = 0;
                 qmu              = 0;
-                amr                               = 0;
-
+                amr                              = 0;
                 results          = 0;
                 outputdefinition = 0;
@@ -128 +126 @@
                         %2016 October 11
                         if isa(md.esa,'double'); md.esa=esa(); end
-                        %2017 Dec 21st (needs to be here)
-                        if isempty(md.settings)
-                                disp('Warning: md.settings had to be reset, make sure to adjust md.settings.output_frequency and other fields');
-                                md.settings = issmsettings();
-                        end
                         %2017 February 10th
                         if md.settings.solver_residue_threshold==0,
@@ -147 +140 @@
                                 disp('Warning: calvingdev is now calvingvonmises');
                                 md.calving=calvingvonmises(md.calving); 
+                        end
+                        %2017 Dec 21st (needs to be here)
+                        if isempty(md.settings)
+                                disp('Warning: md.settings had to be reset, make sure to adjust md.settings.output_frequency and other fields');
+                                md.settings = issmsettings();
                         end
 
@@ -328 +326 @@
                         end
 
+                        %lat long
+                        if numel(md.mesh.lat) ==md.mesh.numberofvertices,  md.mesh.lat=project2d(md,md.mesh.lat,1); end
+                        if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end
+
                         %outputdefinitions
                         for i=1:length(md.outputdefinition.definitions)
@@ -361 +363 @@
 
                 end % }}}
-                function md2 = extract(md,area) % {{{
+                function md2 = extract(md,area,varargin) % {{{
                         %extract - extract a model according to an Argus contour or flag list
                         %
@@ -383 +385 @@
                         md1=md;
 
+                        %recover optoins: 
+                        options=pairoptions(varargin{:});
+
                         %some checks
-                        if ((nargin~=2) | (nargout~=1)),
+                        if ((nargin<2) | (nargout~=1)),
                                 help extract
                                 error('extract error message: bad usage');
@@ -396 +401 @@
 
                         %kick out all elements with 3 dirichlets
-                        spc_elem=find(~flag_elem);
-                        spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
-                        flag=ones(md1.mesh.numberofvertices,1);
-                        flag(spc_node)=0;
-                        pos=find(sum(flag(md1.mesh.elements),2)==0);
-                        flag_elem(pos)=0;
+                        if getfieldvalue(options,'spccheck',1)
+                                spc_elem=find(~flag_elem);
+                                spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
+                                flag=ones(md1.mesh.numberofvertices,1);
+                                flag(spc_node)=0;
+                                pos=find(sum(flag(md1.mesh.elements),2)==0);
+                                flag_elem(pos)=0;
+                        end
 
                         %extracted elements and nodes lists
@@ -827 +834 @@
                         md.calving=extrude(md.calving,md);
                         md.hydrology = extrude(md.hydrology,md);
+                        md.slr = extrude(md.slr,md);
 
                         %connectivity
@@ -1167 +1175 @@
                         md.levelset         = levelset();
                         md.calving          = calving();
-                        md.gia              = giaivins();
+                        md.gia                            = giaivins();
+                        md.esa              = esa();
                         md.love             = fourierlove();
                         md.esa              = esa();
@@ -1173 +1182 @@
                         md.inversion        = inversion();
                         md.qmu              = qmu();
-                        md.amr              = amr();
+                        md.amr                            = amr();
                         md.radaroverlay     = radaroverlay();
                         md.results          = struct();
@@ -1341 +1350 @@
                         disp(sprintf('%19s: %-22s -- %s','levelset'        ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)'));
                         disp(sprintf('%19s: %-22s -- %s','calving'         ,['[1x1 ' class(self.calving) ']'],'parameters for calving'));
-                        disp(sprintf('%19s: %-22s -- %s','gia'        ,['[1x1 ' class(self.gia) ']'],'parameters for gia solution'));
-                        disp(sprintf('%19s: %-22s -- %s','love'        ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
+                        disp(sprintf('%19s: %-22s -- %s','gia'             ,['[1x1 ' class(self.gia) ']'],'parameters for gia solution'));
                         disp(sprintf('%19s: %-22s -- %s','esa'             ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
+                        disp(sprintf('%19s: %-22s -- %s','love'            ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
                         disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
                         disp(sprintf('%19s: %-22s -- %s','inversion'       ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
@@ -1443 +1452 @@
 
                 end % }}}
-                function savemodeljs(md,modelname,websiteroot) % {{{
-
+                function savemodeljs(md,modelname,websiteroot,varargin) % {{{
+                
                         %the goal of this routine is to save the model as a javascript array that can be included in any html 
                         %file: 
 
+                        options=pairoptions(varargin{:});
+                        optimization=getfieldvalue(options,'optimize',0);
+
+                        
                         %disp: 
                         disp(['saving model ''' modelname ''' in file ' websiteroot '/js/' modelname '.js']);
@@ -1466 +1479 @@
                                 end
 
+                                %some optimization: 
+                                if optimization==1,
+                                        %optimize for plotting only:
+                                        if ~ismember(field,{'geometry','mesh','mask'}),
+                                                continue;
+                                        end
+                                end
+
                                 %Check that current field is an object
                                 if ~isobject(md.(field))

issm/trunk-jpl/src/m/classes/model.py

@@ -111 +111 @@
                 self.levelset         = levelset()
                 self.calving          = calving()
+                self.gia              = giaivins()
                 self.love             = fourierlove()
-                self.gia              = giaivins()
                 self.esa              = esa()
-
                 self.autodiff         = autodiff()
                 self.inversion        = inversion()
@@ -128 +127 @@
         def properties(self):    # {{{
                 # ordered list of properties since vars(self) is random
-                return ['mesh',
+                return ['mesh',
                         'mask',
                         'geometry',

issm/trunk-jpl/src/m/classes/sealevelmodel.m

@@ -18 +18 @@
                 settings         = 0;
                 private          = 0;
+                range            = 0;
+                mergedcaps     = 0;
                 %}}}
         end
@@ -43 +45 @@
                         for i=1:length(slm.icecaps),
                                 if slm.icecaps{i}.transient.iscoupler==0,
-                                        error(sprintf('sealevelmodel checkconsistenty error:  icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name));
+                                        warning(sprintf('sealevelmodel checkconsistenty error:  icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name));
                                 end
                         end
                                 
                         if slm.earth.transient.iscoupler==0,
-                                error('sealevelmodel checkconsistenty error:  earth model should have the transient coupler option turned on!');
+                                warning('sealevelmodel checkconsistenty error:  earth model should have the transient coupler option turned on!');
                         end
 
@@ -54 +56 @@
                         for i=1:length(slm.icecaps),
                                 if slm.icecaps{i}.mesh.numberofvertices ~= length(slm.earth.slr.transitions{i}),
-                                        error('sealevelmodel checkconsistenty issue with size of transition vectors!');
+                                        error(['sealevelmodel checkconsistenty issue with size of transition vector for ice cap: ' num2str(i) ' name: ' slm.icecaps{i}.miscellaneous.name]);
+                                end
+                        end
+                        
+                        %check that run_frequency is the same everywhere: 
+                        for i=1:length(slm.icecaps),
+                                if slm.icecaps{i}.slr.geodetic_run_frequency~=slm.earth.slr.geodetic_run_frequency,
+                                        error(sprintf('sealevelmodel checkconsistenty error:  icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name));
                                 end
                         end
 
+                        %make sure steric_rate is the same everywhere: 
+                        for i=1:length(slm.icecaps),
+                                md= slm.icecaps{i}; 
+                                if ~isempty(find(md.slr.steric_rate - slm.earth.slr.steric_rate(slm.earth.slr.transitions{i}))),
+                                        error(sprintf('steric rate on ice cap %s is not the same as for the earth\n',md.miscellaneous.name));
+                                end
+                        end
 
                 end
@@ -70 +86 @@
                         slm.private           = private();
                         slm.cluster           = generic();
+                        slm.range             = {};
                 end
                 %}}}
@@ -78 +95 @@
                         disp(sprintf('%19s: %-22s -- %s','cluster'         ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of cpus...)'));
                         disp(sprintf('%19s: %-22s -- %s','miscellaneous'   ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
+                        disp(sprintf('%19s: %-22s -- %s','range'   ,['[1x1 ' class(self.range) ']'],'ranges'));
+                end % }}}
+                function self=mergeresults(self) % {{{
+                        champs=fields(self.icecaps{1}.results.TransientSolution);
+                        for i=1:length(self.mergedcaps)/2,
+                                md=self.mergedcaps{2*(i-1)+1}; trans=self.mergedcaps{2*(i-1)+2};
+                                icecaps=self.icecaps(self.range{2*(i-1)+2});
+                                for j=1:length(self.icecaps{1}.results.TransientSolution),
+                                        for k=1:length(champs),
+                                                if strcmpi(class(icecaps{1}.results.TransientSolution(j).(champs{k})),'double'),
+                                                        %vertex or element? 
+                                                        if length(icecaps{1}.results.TransientSolution(j).(champs{k}))==icecaps{1}.mesh.numberofvertices,
+                                                                md.results.TransientSolution(j).(champs{k})=zeros(md.mesh.numberofvertices,1);
+                                                                for l=1:length(trans),
+                                                                        resultcap=icecaps{l}.results.TransientSolution(j).(champs{k});
+                                                                        md.results.TransientSolution(j).(champs{k})(trans{l})=resultcap;
+                                                                end
+                                                        else 
+                                                                if strcmpi(champs{k},'IceVolume') | strcmpi(champs{k},'IceVolumeAboveFloatation') ,
+                                                                        md.results.TransientSolution(j).(champs{k})=0;
+                                                                        for l=1:length(trans),
+                                                                                resultcap=icecaps{l}.results.TransientSolution(j).(champs{k});
+                                                                                md.results.TransientSolution(j).(champs{k})= md.results.TransientSolution(j).(champs{k})+resultcap;
+                                                                        end
+                                                                elseif strcmpi(champs{k},'time'),
+                                                                        md.results.TransientSolution(j).(champs{k})= icecaps{1}.results.TransientSolution(j).(champs{k});
+                                                                else
+                                                                        continue;
+                                                                end
+                                                        end
+                                                else
+                                                        continue;
+                                                end
+                                                                                        end
+                                end
+                                self.mergedcaps{2*(i-1)+1}=md;
+                        end
+                end % }}}
+                function listcaps(self) % {{{
+                        for  i=1:length(self.icecaps),
+                                disp(sprintf('%i: %s',i,self.icecaps{i}.miscellaneous.name));
+                        end
+                end % }}}
+                function viscousiterations(self) % {{{
+                        for  i=1:length(self.icecaps),
+                                ic=self.icecaps{i}; 
+                                mvi=ic.results.TransientSolution(1).StressbalanceConvergenceNumSteps;
+                                for j=2:length(ic.results.TransientSolution)-1,
+                                        mvi=max(mvi,ic.results.TransientSolution(j).StressbalanceConvergenceNumSteps);
+                                end
+                                disp(sprintf('%i, %s: %i',i,self.icecaps{i}.miscellaneous.name,mvi));
+                        end
+                end % }}}
+                function maxtimestep(self) % {{{
+                        for  i=1:length(self.icecaps),
+                                ic=self.icecaps{i}; 
+                                mvi=length(ic.results.TransientSolution);
+                                timei=ic.results.TransientSolution(end).time;
+                                disp(sprintf('%i, %s: %i/%g',i,self.icecaps{i}.miscellaneous.name,mvi,timei));
+                        end
+                        mvi=length(self.earth.results.TransientSolution);
+                        timei=self.earth.results.TransientSolution(end).time;
+                        disp(sprintf('Earth: %i/%g',mvi,timei));
+                end % }}}
+                function self=homogeneize(self,noearth) % {{{
+                        if nargin==1,
+                                noearth=0;
+                        end
+                        mintimestep=Inf;
+                        for  i=1:length(self.icecaps),
+                                ic=self.icecaps{i}; 
+                                mintimestep=min(mintimestep, length(ic.results.TransientSolution));
+                        end
+                        if ~noearth,
+                                mintimestep=min(mintimestep, length(self.earth.results.TransientSolution));
+                        end
+                        
+                        for  i=1:length(self.icecaps),
+                                ic=self.icecaps{i}; 
+                                ic.results.TransientSolution=ic.results.TransientSolution(1:mintimestep);
+                                self.icecaps{i}=ic;
+                        end
+                        ic=self.earth;
+                        if ~noearth,
+                                ic.results.TransientSolution=ic.results.TransientSolution(1:mintimestep);
+                        end
+                        self.earth=ic;
                 end % }}}
         end

issm/trunk-jpl/src/m/classes/slr.js

@@ -88 +88 @@
                         console.log(sprintf('   Sealevelrise solution parameters:'));
 
-                        fielddisplay(this,'deltathickness','thickness change: ice height equivalent [m]');
-                        fielddisplay(this,'sealevel','current sea level (prior to computation) [m]');
-                        fielddisplay(this,'reltol','sea level rise relative convergence criterion, (NaN: not applied)');
-                        fielddisplay(this,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)');
-                        fielddisplay(this,'maxiter','maximum number of nonlinear iterations');
-                        fielddisplay(this,'love_h','load Love number for radial displacement');
-                        fielddisplay(this,'love_k','load Love number for gravitational potential perturbation');
-                        fielddisplay(this,'love_l','load Love number for horizontal displacements'); 
-                        fielddisplay(this,'tide_love_h','tidal love number (degree 2)');
-                        fielddisplay(this,'tide_love_k','tidal love number (degree 2)');
-                        fielddisplay(this,'fluid_love','secular fluid Love number');
-                        fielddisplay(this,'equatorial_moi','mean equatorial moment of inertia [kg m^2]');
-                        fielddisplay(this,'polar_moi','polar moment of inertia [kg m^2]');
-                        fielddisplay(this,'angular_velocity','mean rotational velocity of earth [per second]'); 
-                        fielddisplay(this,'rigid','rigid earth graviational potential perturbation');
-                        fielddisplay(this,'elastic','elastic earth graviational potential perturbation');
-                        fielddisplay(this,'rotation','rotational earth potential perturbation');
-                        fielddisplay(this,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); 
-                        fielddisplay(this,'steric_rate','rate of steric ocean expansion [mm/yr]'); 
-                        fielddisplay(this,'degacc',"accuracy (default .01 deg) for numerical discretization of the Green's functions");
-                        fielddisplay(this,'transitions','indices into parts of the mesh that will be icecaps');
-                        fielddisplay(this,'requested_outputs','additional outputs requested');
+                fielddisplay(this,'deltathickness','thickness change (main loading of the slr solution core [m]');
+                fielddisplay(this,'sealevel','current sea level (prior to computation) [m]');
+                fielddisplay(this,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)');
+                fielddisplay(this,'abstol','sea level rise absolute convergence criterion, NaN: not applied');
+                fielddisplay(this,'maxiter','maximum number of nonlinear iterations');
+                fielddisplay(this,'love_h','load Love number for radial displacement');
+                fielddisplay(this,'love_k','load Love number for gravitational potential perturbation');
+                fielddisplay(this,'love_l','load Love number for horizontal displacements'); 
+                fielddisplay(this,'tide_love_h','tidal love number (degree 2)');
+                fielddisplay(this,'tide_love_k','tidal love number (degree 2)');
+                fielddisplay(this,'fluid_love','secular fluid Love number');
+                fielddisplay(this,'equatorial_moi','mean equatorial moment of inertia [kg m^2]');
+                fielddisplay(this,'polar_moi','polar moment of inertia [kg m^2]');
+                fielddisplay(this,'angular_velocity','mean rotational velocity of earth [per second]'); 
+                fielddisplay(this,'rigid','rigid earth graviational potential perturbation');
+                fielddisplay(this,'elastic','elastic earth graviational potential perturbation');
+                fielddisplay(this,'rotation','rotational earth potential perturbation');
+                fielddisplay(this,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); 
+                fielddisplay(this,'degacc',"accuracy (default .01 deg) for numerical discretization of the Green's functions");
+                fielddisplay(this,'transitions','indices into parts of the mesh that will be icecaps');
+                fielddisplay(this,'requested_outputs','additional outputs requested');
+                fielddisplay(this,'steric_rate','rate of steric ocean expansion [mm/yr]'); 
                 } //}}}
                 this.marshall=function(md,prefix,fid) { //{{{

issm/trunk-jpl/src/m/classes/slr.m

@@ -8 +8 @@
                 deltathickness = NaN;
                 sealevel       = NaN; 
+                spcthickness   = NaN; 
                 maxiter        = 0;
                 reltol         = 0;
@@ -25 +26 @@
                 ocean_area_scaling = 0;
                 steric_rate    = 0; %rate of ocean expansion from steric effects. 
+                geodetic_run_frequency  = 1; %how many time steps we skip before we run the geodetic part of the solver during transient
+                geodetic       = 0; %compute geodetic SLR? (in addition to steric?)
                 degacc         = 0;
+                loop_increment = 0;
+                horiz = 0;
+                Ngia  = 0;
+                Ugia  = 0;
                 requested_outputs      = {};
                 transitions    = {};
@@ -46 +53 @@
                 %maximum of non-linear iterations.
                 self.maxiter=5;
+                self.loop_increment=200;
 
                 %computational flags: 
+                self.geodetic=0;
                 self.rigid=1;
                 self.elastic=1;
-                self.rotation=0;
                 self.ocean_area_scaling=0;
+                self.rotation=1;
 
                 %tidal love numbers: 
@@ -73 +82 @@
                 self.steric_rate=0;
         
+                %how many time steps we skip before we run SLR solver during transient
+                self.geodetic_run_frequency=1;
         
                 %output default:
@@ -79 +90 @@
                 %transitions should be a cell array of vectors: 
                 self.transitions={};
+
+                %horizontal displacement?  (not by default)
+                self.horiz=0;
                 
                 end % }}}
@@ -86 +100 @@
                         md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
                         md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
+                        md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1);
                         md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
                         md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1);
@@ -98 +113 @@
                         md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
                         md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
+                        md = checkfield(md,'fieldname','slr.geodetic_run_frequency','size',[1 1],'>=',1);
                         md = checkfield(md,'fieldname','slr.steric_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
                         md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
                         md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1);
+                        md = checkfield(md,'fieldname','slr.loop_increment','NaN',1,'Inf',1,'>=',1);
+                        md = checkfield(md,'fieldname','slr.horiz','NaN',1,'Inf',1,'values',[0 1]);
+                        md = checkfield(md,'fieldname','slr.Ngia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
+                        md = checkfield(md,'fieldname','slr.Ugia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
 
                         %check that love numbers are provided at the same level of accuracy: 
@@ -112 +132 @@
                         [els,vertices]=find(maskpos>0);
                         if length(els),
-                                error('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
-                        end
+                                warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
+                        end
+
+                        %check that  if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not, 
+                        %a coupler to a planet model is provided. 
+                        if self.geodetic,
+                                if md.transient.iscoupler, 
+                                        %we are good; 
+                                else
+                                        if strcmpi(class(md.mesh),'mesh3dsurface'),
+                                                %we are good
+                                        else
+                                                error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!');
+                                        end
+                                end
+                        end
+
 
                 end % }}}
@@ -124 +159 @@
                         fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]');
                         fielddisplay(self,'sealevel','current sea level (prior to computation) [m]');
-                        fielddisplay(self,'reltol','sea level rise relative convergence criterion, (NaN: not applied)');
-                        fielddisplay(self,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)');
+                        fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]');
+                        fielddisplay(self,'reltol','sea level rise relative convergence criterion, (default, NaN: not applied)');
+                        fielddisplay(self,'abstol','sea level rise absolute convergence criterion, NaN: not applied');
                         fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
                         fielddisplay(self,'love_h','load Love number for radial displacement');
@@ -136 +172 @@
                         fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]');
                         fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]'); 
+                        fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'); 
+                        fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)'); 
+                        fielddisplay(self,'Ngia','rate of viscous (GIA) geoid expansion (in mm/yr)'); 
+                        fielddisplay(self,'Ugia','rate of viscous (GIA) bedrock uplift (in mm/yr)'); 
+                        fielddisplay(self,'loop_increment','vector assembly (in the convolution) framentation');
+                        fielddisplay(self,'geodetic','compute geodetic SLR? ( in addition to steric?) default 0');
+                        fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)');
                         fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
                         fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
                         fielddisplay(self,'rotation','earth rotational potential perturbation');
-                        fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'); 
-                        fielddisplay(self,'steric_rate','rate of steric ocean expansion [mm/yr]'); 
                         fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
                         fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
@@ -149 +190 @@
                         %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2);
                         WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
+                        %WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofelements+1);
                         WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
+                        WriteData(fid,prefix,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
                         WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double');
                         WriteData(fid,prefix,'object',self,'fieldname','abstol','format','Double');
@@ -166 +209 @@
                         WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean');
                         WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean');
+                        WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer');
                         WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
+                        WriteData(fid,prefix,'object',self,'fieldname','Ngia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
+                        WriteData(fid,prefix,'object',self,'fieldname','Ugia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
                         WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
                         WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
+                        WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer');
+                        WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer');
+                        WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer');
                         
                         %process requested outputs
@@ -184 +233 @@
                         writejs1Darray(fid,[modelname '.slr.deltathickness'],self.deltathickness);
                         writejs1Darray(fid,[modelname '.slr.sealevel'],self.sealevel);
+                        writejs1Darray(fid,[modelname '.slr.spcthickness'],self.spcthickness);
                         writejsdouble(fid,[modelname '.slr.maxiter'],self.maxiter);
                         writejsdouble(fid,[modelname '.slr.reltol'],self.reltol);
@@ -200 +250 @@
                         writejsdouble(fid,[modelname '.slr.rotation'],self.rotation);
                         writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling);
+                        writejsdouble(fid,[modelname '.slr.geodetic_run_frequency'],self.geodetic_run_frequency);
                         writejs1Darray(fid,[modelname '.slr.steric_rate'],self.steric_rate);
                         writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
@@ -205 +256 @@
                         writejscellarray(fid,[modelname '.slr.transitions'],self.transitions);
                 end % }}}
+                function self = extrude(self,md) % {{{
+                        self.sealevel=project3d(md,'vector',self.sealevel,'type','node');
+                end % }}}
         end
 end

issm/trunk-jpl/src/m/classes/slr.py

@@ -60 +60 @@
                         string="%s\n%s"%(string,fielddisplay(self,'elastic','elastic earth graviational potential perturbation'))
                         string="%s\n%s"%(string,fielddisplay(self,'rotation','earth rotational potential perturbation'))
-                        string="%s\n%s"%(string,fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)'))
+                        string="%s\n%s"%(string,fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'))
                         string="%s\n%s"%(string,fielddisplay(self,'steric_rate','rate of steric ocean expansion [mm/yr]'))
                         string="%s\n%s"%(string,fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'))
@@ -71 +71 @@
                 
                 #Convergence criterion: absolute, relative and residual
-                self.reltol=0.01 # 1 per cent 
-                self.abstol=float('NaN') #default
+                self.reltol=float('NaN') #default
+                self.abstol=0.001 #1 mm of sea level rise
 
                 #maximum of non-linear iterations.
@@ -87 +87 @@
                 self.tide_love_k=0.3055; #degree 2
                 
-      #secular fluid love number: 
+                #secular fluid love number: 
                 self.fluid_love=0.942; 
                 
@@ -93 +93 @@
                 self.equatorial_moi=8.0077*10**37; # [kg m^2] 
                 self.polar_moi     =8.0345*10**37; # [kg m^2] 
-
+                
                 #mean rotational velocity of earth 
                 self.angular_velocity=7.2921*10**-5; # [s^-1]