source: issm/oecreview/Archive/25834-26739/ISSM-26109-26110.diff@ 26740

../trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h

@@ -51 +51 @@
                 virtual void PointwiseDivide(IssmAbsVec* x,IssmAbsVec* y)=0;
                 virtual void PointwiseMult(IssmAbsVec* x,IssmAbsVec* y)=0;
                 virtual void Pow(doubletype scale_factor)=0;
+                virtual void Sum(doubletype* pvalue)=0;
 };
 
 #endif //#ifndef _ISSM_ABS_VEC_H_

../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h

@@ -593 +593 @@
 
                 }
                 /*}}}*/
+                void Sum(doubletype* pvalue){/*{{{*/
+                        _error_("not support yet!");
+                }
+                /*}}}*/
                 void BucketsBuildScatterBuffers(int** pnumvalues_forcpu,int** prow_indices_forcpu,doubletype** pvalues_forcpu,int** pmodes_forcpu,DataSet** bucketsforcpu,int num_procs){/*{{{*/
 
                         /*intermediary: */

../trunk-jpl/src/c/toolkits/issm/IssmVec.h

@@ -219 +219 @@
                         vector->Pow(scale_factor);
                 }
                 /*}}}*/
+                void Sum(doubletype*  pvalue){/*{{{*/
+                        vector->Sum(pvalue);
+                }
+                /*}}}*/
 };
 
 #endif //#ifndef _ISSMVEC_H_

../trunk-jpl/src/c/toolkits/objects/Vector.h

@@ -405 +405 @@
                         else this->ivector->Pow(scale_factor);
                 }
                 /*}}}*/
-};
+void Sum(doubletype* pvalue){ /*{{{*/
+        _assert_(this);/*{{{*/
+
+        if(type==PetscVecType){
+                #ifdef _HAVE_PETSC_
+                this->pvector->Sum(pvalue);
+                #endif
+        }
+        else this->ivector->Sum(pvalue);
+}
+/*}}}*/
+}; /*}}}*/
 #endif //#ifndef _VECTOR_H_

../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h

@@ -57 +57 @@
                 IssmDouble  Max(void);
                 void        Scale(IssmDouble scale_factor);
                 void        Pow(IssmDouble scale_factor);
+                void        Sum(IssmDouble* pvalue);
                 void        PointwiseDivide(PetscVec* x,PetscVec* y);
                 void        PointwiseMult(PetscVec* x,PetscVec* y);
                 IssmDouble  Dot(PetscVec* vector);

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

@@ -96 +96 @@
         ./classes/Vertex.cpp \
         ./classes/Hook.cpp \
         ./classes/Radar.cpp \
+        ./classes/BarystaticContributions.cpp \
         ./classes/ExternalResults/Results.cpp \
         ./classes/Elements/Element.cpp \
         ./classes/Elements/Elements.cpp \

../trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp

@@ -182 +182 @@
                 parameters->AddObject(new DoubleMatParam(CumBslcOceanPartitionEnum,bslcocean_partition,npartocean,1));
                 xDelete<IssmDouble>(partitionocean);
         }
-
+        /*New optimized code:*/
+        BarystaticContributions* barystaticcontributions=new BarystaticContributions(iomodel);
+        parameters->AddObject(new GenericParam<BarystaticContributions*>(barystaticcontributions,BarystaticContributionsEnum));
         
         /*Deal with external multi-model ensembles: {{{*/
         if(isexternal){

../trunk-jpl/src/c/classes/BarystaticContributions.cpp

@@ +1 @@
+/*!\file BarystaticContributions.c
+ * \brief: implementation of the BarystaticContributions object
+ */
+
+/*Include files: {{{*/
+#ifdef HAVE_CONFIG_H
+        #include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "./BarystaticContributions.h" 
+#include "../toolkits/toolkits.h"
+#include "./classes.h"
+/*}}}*/
+
+/*Constructors and destructors:*/
+BarystaticContributions::BarystaticContributions(IoModel* iomodel ){ /*{{{*/
+
+        int nel;
+
+        iomodel->FetchData(&nice,"md.solidearth.npartice");
+        if(nice){
+                iomodel->FetchData(&pice,&nel,NULL,"md.solidearth.partitionice");
+                ice=new Vector<IssmDouble>(nice);
+                cumice=new Vector<IssmDouble>(nice); cumice->Set(0); cumice->Assemble();
+        }
+
+        iomodel->FetchData(&nhydro,"md.solidearth.nparthydro");
+        if(nhydro){
+                iomodel->FetchData(&phydro,&nel,NULL,"md.solidearth.partitionhydro");
+                hydro=new Vector<IssmDouble>(nhydro);
+                cumhydro=new Vector<IssmDouble>(nhydro); cumhydro->Set(0); cumhydro->Assemble();
+        }
+        iomodel->FetchData(&nocean,"md.solidearth.npartocean");
+        if(nocean){
+                iomodel->FetchData(&pocean,&nel,NULL,"md.solidearth.partitionocean");
+                ocean=new Vector<IssmDouble>(nocean);
+                cumocean=new Vector<IssmDouble>(nocean); cumocean->Set(0); cumocean->Assemble();
+        }
+
+} /*}}}*/
+BarystaticContributions::~BarystaticContributions(){ /*{{{*/
+        delete ice;   delete cumice;
+        delete hydro; delete cumhydro;
+        delete ocean; delete cumocean;
+        if(nice)xDelete<IssmDouble>(pice);
+        if(nhydro)xDelete<IssmDouble>(phydro);
+        if(nocean)xDelete<IssmDouble>(pocean);
+}; /*}}}*/
+
+/*Support routines:*/
+IssmDouble BarystaticContributions::Total(){ /*{{{*/
+
+        IssmDouble  sumice,sumhydro,sumocean;
+        
+        ice->Assemble();
+        hydro->Assemble();
+        ocean->Assemble();
+
+        ice->Sum(&sumice);
+        hydro->Sum(&sumhydro);
+        ocean->Sum(&sumocean);
+
+        return sumice+sumhydro+sumocean;
+
+} /*}}}*/
+IssmDouble BarystaticContributions::CumTotal(){ /*{{{*/
+
+        IssmDouble sumice,sumhydro,sumocean;
+
+        cumice->Assemble();
+        cumhydro->Assemble();
+        cumocean->Assemble();
+
+        cumice->Sum(&sumice);
+        cumhydro->Sum(&sumhydro);
+        cumocean->Sum(&sumocean);
+
+
+        return sumice+sumhydro+sumocean;
+
+} /*}}}*/
+void BarystaticContributions::Cumulate(Parameters* parameters){ /*{{{*/
+
+        cumice->AXPY(ice,1);
+        cumocean->AXPY(ocean,1);
+        cumhydro->AXPY(hydro,1);
+
+
+} /*}}}*/
+void BarystaticContributions::Save(Results* results, Parameters* parameters, IssmDouble oceanarea){ /*{{{*/
+
+        int        step;
+        IssmDouble time;
+        IssmDouble rho_water;
+
+        IssmDouble* cumice_serial=NULL;
+        IssmDouble* cumhydro_serial=NULL;
+        IssmDouble* cumocean_serial=NULL;
+
+        IssmDouble sumice,sumhydro,sumocean;
+
+        parameters->FindParam(&step,StepEnum);
+        parameters->FindParam(&time,TimeEnum);
+        parameters->FindParam(&rho_water,TimeEnum);
+
+        ice->Sum(&sumice); hydro->Sum(&sumhydro); ocean->Sum(&sumocean);
+        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcEnum,-this->Total()/oceanarea/rho_water,step,time));
+        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcIceEnum,-sumice/oceanarea/rho_water,step,time));
+        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcHydroEnum,-sumice/oceanarea/rho_water,step,time));
+        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,BslcOceanEnum,-sumocean/oceanarea/rho_water,step,time));
+
+        cumice->Sum(&sumice); cumhydro->Sum(&sumhydro); cumocean->Sum(&sumocean);
+        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcEnum,this->CumTotal()/oceanarea/rho_water,step,time));
+        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcIceEnum,sumice/oceanarea/rho_water,step,time));
+        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcHydroEnum,sumhydro/oceanarea/rho_water,step,time));
+        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcOceanEnum,sumocean/oceanarea/rho_water,step,time));
+
+        cumice_serial=this->cumice->ToMPISerial0(); for (int i=0;i<nice;i++)cumice_serial[i]=cumice_serial[i]/oceanarea/rho_water;
+        cumhydro_serial=this->cumhydro->ToMPISerial0(); for (int i=0;i<nhydro;i++)cumhydro_serial[i]=cumhydro_serial[i]/oceanarea/rho_water;
+        cumocean_serial=this->cumocean->ToMPISerial0(); for (int i=0;i<nocean;i++)cumocean_serial[i]=cumocean_serial[i]/oceanarea/rho_water;
+        
+        results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcIcePartitionEnum,cumice_serial,nice,1,step,time));
+        results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcHydroPartitionEnum,cumhydro_serial,nhydro,1,step,time));
+        results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcOceanPartitionEnum,cumocean_serial,nocean,1,step,time));
+
+        if(IssmComm::GetRank()==0){
+                xDelete<IssmDouble>(cumice_serial);
+                xDelete<IssmDouble>(cumhydro_serial);
+                xDelete<IssmDouble>(cumocean_serial);
+        }
+        return;
+
+} /*}}}*/
+void BarystaticContributions::Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue){ /*{{{*/
+        
+        int id;
+                
+        id=reCast<int>(pice[eid]);
+        ice->SetValue(id,icevalue,ADD_VAL);
+
+        id=reCast<int>(phydro[eid]);
+        hydro->SetValue(id,hydrovalue,ADD_VAL);
+
+        id=reCast<int>(pocean[eid]);
+        ocean->SetValue(id,oceanvalue,ADD_VAL);
+
+} /*}}}*/

../trunk-jpl/src/c/classes/BarystaticContributions.h

@@ +1 @@
+/*!\file BarystaticContributions.h
+ * \brief: header file for barystatic contribution object
+ */
+
+#ifndef _BARYSTATICCONTRIBUTIONS_H_
+#define _BARYSTATICCONTRIBUTIONS_H_
+
+/*Headers:*/
+class IoModel;
+class Parameters;
+class Results;
+template <class doubletype> class Vector;
+#include "../shared/shared.h"
+
+class BarystaticContributions {
+
+        public: 
+
+                Vector<IssmDouble>* ice;  //contributions to every ice partition
+                Vector<IssmDouble>* cumice;  //cumulated contributions to every ice partition
+                int                 nice; //number of ice partitions 
+                IssmDouble*         pice; //ice partition
+
+                Vector<IssmDouble>* hydro;  //contributions to every hydro partition
+                Vector<IssmDouble>* cumhydro;  //cumulated contributions to every hydro partition
+                int                 nhydro; //number of hydro partitions 
+                IssmDouble*         phydro; //hydro partition
+
+                Vector<IssmDouble>* ocean;  //contributions to every ocean partition
+                Vector<IssmDouble>* cumocean;  //cumulated contributions to every ocean partition
+                int                 nocean; //number of ocean partitions 
+                IssmDouble*         pocean; //ocean partition
+
+                /*BarystaticContributions constructors, destructors :*/
+                BarystaticContributions(IoModel* iomodel );
+                ~BarystaticContributions();
+
+                /*routines:*/
+                IssmDouble Total();
+                IssmDouble CumTotal();
+                void Cumulate(Parameters* parameters);
+                void Save(Results* results, Parameters* parameters, IssmDouble oceanarea);
+                void Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue);
+
+};
+#endif  /* _BARYSTATICCONTRIBUTIONS_H_ */

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

@@ -37 +37 @@
 template <class doubletype> class Vector;
 class ElementMatrix;
 class ElementVector;
+class BarystaticContributions;
 /*}}}*/
 
 class Element: public Object{
@@ -388 +389 @@
                 virtual void          SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* mask)=0;
                 virtual void          DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0;
                 virtual void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0;
+
+                virtual void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0;
+                virtual void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0;
+                virtual void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks)=0;
+                virtual void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks)=0;
+                virtual void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0;
                 #endif
 
 };

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

@@ -226 +226 @@
                 void    SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
                 void    DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
                 void           GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
+
+                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
+                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
+                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
                 #endif
 
                 /*}}}*/

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

@@ -181 +181 @@
                 void    DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
                 IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
                 void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
+
+                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
+                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
+                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
+
 #endif
 
 #ifdef _HAVE_DAKOTA_

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

@@ -188 +188 @@
                 void    DeformationFromSurfaceLoads(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
                 IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
                 void        GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
+                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
+                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
+                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
+                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
+
 #endif
 
 #ifdef _HAVE_DAKOTA_

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

@@ -5431 +5431 @@
 /*}}}*/
 #endif
 #ifdef _HAVE_SEALEVELCHANGE_
+//old code
 IssmDouble Tria::OceanAverage(IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
 
         if(masks->isoceanin[this->lid]){
@@ -5450 +5451 @@
 
 }
 /*}}}*/
-void    Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
+void          Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
         /*early return if we are not on an ice cap OR ocean:*/
         if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){
                 dI_list[0] = 0.0; // this is important!!!
@@ -5544 +5545 @@
 
         return;
 }/*}}}*/
-void    Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
+void       Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
 
         masks->isiceonly[this->lid]=this->IsIceOnlyInElement();
         masks->isoceanin[this->lid]=this->IsOceanInElement();
@@ -5557 +5558 @@
         /*are we not fully grounded: */
         if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true;
         else masks->notfullygrounded[this->lid]=false;
-
 }
 /*}}}*/
-void    Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
+void       Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
 
         /*diverse:*/
         int gsize;
@@ -5950 +5950 @@
         return bslchydro;
 }
 /*}}}*/
-void    Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/
+void       Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/
 
         /*diverse:*/
         int gsize;
@@ -6000 +6000 @@
         return;
 }
 /*}}}*/
-void    Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/
+void       Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/
 
         /*diverse:*/
         int gsize,dummy;
@@ -6042 +6042 @@
         return;
 }
 /*}}}*/
-void    Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
+void       Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
 
         /*diverse:*/
         int gsize;
@@ -6131 +6131 @@
         return;
 }
 /*}}}*/
-void    Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/
+void       Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/
 
         IssmDouble xyz_list[NUMVERTICES][3];
 
@@ -6237 +6237 @@
 }
 /*}}}*/
 
-//new logic 
-void    Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
+//new code 
+void       Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/
 
         /*diverse:*/
         int nel;
@@ -6427 +6427 @@
         return;
 }
 /*}}}*/
-void    Tria::SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix<IssmDouble>* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/
+void       Tria::SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){ /*{{{*/
 
         /*diverse:*/
         IssmDouble area;
@@ -6435 +6435 @@
         IssmDouble phi_water=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
         IssmDouble I,W,BP;  //change in ice thickness, water column or bottom pressure (Farrel and Clarke, Equ. 4)
         bool notfullygrounded=false;
-        bool scaleoceanarea= false;
         bool computerigid= false;
         int  glfraction=1;
         int  npartice,nparthydro,npartocean;
@@ -6468 +6467 @@
                 #ifdef _ISSM_DEBUG_
                         constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
                 #endif
-                /*barystatic_contribution[0]+=0;
-                barystatic_contribution[1]+=0;
-                barystatic_contribution[2]+=0;*/
                 return;
                 }
         }
@@ -6482 +6478 @@
                         #ifdef _ISSM_DEBUG_
                         this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
                         #endif
-                        /*barystatic_contribution[0]+=0;
-                        barystatic_contribution[1]+=0;
-                        barystatic_contribution[2]+=0;*/
                         return;
                 }
         }
@@ -6493 +6486 @@
          * hydrology:*/
         if(!isice  && !ishydro){
                 if(!masks->isoceanin[this->lid]){
-                        /*barystatic_contribution[0]+=0;
-                        barystatic_contribution[1]+=0;
-                        barystatic_contribution[2]+=0;*/
                         return;
                 }
 
@@ -6508 +6498 @@
         rho_ice=FindParam(MaterialsRhoIceEnum);
         rho_water=FindParam(MaterialsRhoSeawaterEnum);
         this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
-        this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
         this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
         this->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
         this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
@@ -6517 +6506 @@
         this->Element::GetInputValue(&area,AreaEnum);
 
         /*Deal with ice loads if we are on grounded ice:*/
-        if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){ /*{{{*/
+        if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){ 
 
-                /*Compute fraction of the element that is grounded: */
+                /*Compute fraction of the element that is grounded: {{{*/
                 if(notfullygrounded){
                         IssmDouble xyz_list[NUMVERTICES][3];
                         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
@@ -6531 +6520 @@
 #endif
                 }
                 else phi_ice=1.0;
+                /*}}}*/
 
                 /*Inform mask: */
                 constant+=100; //1 for ice.
@@ -6568 +6558 @@
                 }
                 /*}}}*/
 
-                /*Compute barystatic contribution:*/
-                _assert_(oceanarea>0.);
-                if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
-                bslcice = rho_ice*area*I/(oceanarea*rho_water);
+                /*Compute barystatic contribution in kg:*/
+                bslcice = rho_ice*area*I;
                 _assert_(!xIsNan<IssmDouble>(bslcice));
 
                 /*Transfer thickness change into kg/m^2:*/
                 I=I*rho_ice*phi_ice; 
-        } /*}}}*/
+        } 
 
         /*Deal with water loads if we are on ground:*/
         if(!masks->isfullyfloating[this->lid]){
@@ -6586 +6574 @@
                 if (!deltathickness_input)_error_("DeltaTwsEnum input needed to compute sea level change!");
                 deltathickness_input->GetInputAverage(&W);
 
-                /*Compute barystatic component:*/
-                _assert_(oceanarea>0.);
-                if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
-                bslchydro = rho_freshwater*area*phi_water*W/(oceanarea*rho_water);
+                /*Compute barystatic component in kg:*/
+                bslchydro = rho_freshwater*area*phi_water*W;
                 _assert_(!xIsNan<IssmDouble>(bslchydro));
 
                 /*convert from m to kg/m^2:*/
@@ -6604 +6590 @@
                 if (!bottompressure_change_input)_error_("DeltaBottomPressureEnum pressure input needed to compute sea level change fingerprint!");
                 bottompressure_change_input->GetInputAverage(&BP);
                 
-                /*Compute barystatic component:*/
-                bslcbp = rho_water*area*BP/(oceanarea*rho_water);
+                /*Compute barystatic component in kg:*/
+                bslcbp = rho_water*area*BP;
 
                 /*convert from m to kg/m^2:*/
                 BP=BP*rho_water;
@@ -6612 +6598 @@
         }
 
         /*Plug all loads into total load vector:*/
-        localloads[this->lid]+=I+W+BP;
+        loads->SetValue(this->sid,I+W+BP,INS_VAL);
         
-        /*Plug bslcice into barystatic contribution vector:*/
-        if(barystatic_contribution_onpartition){
-                int idi=reCast<int>(partition[this->Sid()])+1;
-                int idj=0;
-                idj=0;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcice,ADD_VAL);
-                idj=1;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslchydro,ADD_VAL);
-                idj=2;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcbp,ADD_VAL);
-        }
+        /*Keep track of barystatic contributions:*/
+        barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
 
-        barystatic_contribution[0]+=bslcice;
-        barystatic_contribution[1]+=bslchydro;
-        barystatic_contribution[2]+=bslcbp;
-
 }
 /*}}}*/
-IssmDouble Tria::SealevelchangeConvolution(IssmDouble* loads){ /*{{{*/
+void       Tria::SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){ /*{{{*/
 
         /*sal green function:*/
         IssmDouble* G=NULL;
-        IssmDouble Sealevel[NUMVERTICES]={0,0,0};
+        IssmDouble SealevelGRD[NUMVERTICES]={0,0,0};
+        IssmDouble oceanaverage,oceanarea=0;
 
         bool sal = false;
         int  size;
@@ -6646 +6623 @@
                 
                 for(int i=0;i<NUMVERTICES;i++) {
                         for (int e=0;e<nel;e++){
-                                Sealevel[i]+=G[i*nel+e]*loads[e]; 
+                                SealevelGRD[i]+=G[i*nel+e]*loads[e]; 
                         }
                 }
+        }
 
-                this->AddInput(SealevelRSLEnum,Sealevel,P1Enum);
+        /*compute ocean average over element:*/
+        OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
+        
+        /*add ocean average in the global sealevelloads vector:*/
+        sealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
+        oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
+        
+        return;
+} /*}}}*/
+void       Tria::SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){ /*{{{*/
+
+        bool converged=false;
+        IssmDouble SealevelGRD[3]={0,0,0};
+        IssmDouble oceanaverage,oceanarea=0;
+        int nel;
+        bool sal = false;
+        IssmDouble* G=NULL;
+        int size;
+        
+        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
+        this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
+                
+        if(sal){
+                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
+
+                for(int i=0;i<NUMVERTICES;i++) {
+                        for (int e=0;e<nel;e++){
+                                SealevelGRD[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]);
+                        }
+                }
         }
+        OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
+        newsealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
 
 } /*}}}*/
+void       Tria::OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/
 
+        IssmDouble phi=1.0;
+        bool iscoastline=false;
+        IssmDouble area;
+        IssmDouble Sg_avg=0; //output
+        
+        /*Do we have an ocean?:*/
+        if(!masks->isoceanin[this->lid]){
+                *poceanarea=0;
+                *poceanaverage=0;
+        }
+
+        /*Do we have a coastline?:*/
+        if(!masks->isfullyfloating[this->lid])iscoastline=true;
+
+        if(iscoastline){
+                IssmDouble xyz_list[NUMVERTICES][3];
+                ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+                phi=1.0-this->GetGroundedPortion(&xyz_list[0][0]); 
+        }
+
+        /*Get area of element:*/
+        this->Element::GetInputValue(&area,AreaEnum);
+
+        /*Average over ocean if there is no coastline, area of the ocean 
+         *is the areaa of the element:*/
+        if(!iscoastline){
+
+                /*Average Sg over vertices:*/
+                for(int i=0;i<NUMVERTICES;i++) Sg_avg+=Sg[i]/NUMVERTICES;
+
+                *poceanaverage=Sg_avg;
+                *poceanarea=area;
+                return;
+        }
+
+        /*Average over  the ocean only if there is a coastline. Area of the 
+         * ocean will be the fraction times the area of the element:*/
+        area=phi*area;
+   
+        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);
+        //!mainlyfloating so that the integration is done on the ocean (and not the grounded) part
+        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;
+        Sg_avg=0;
+        while(gauss->next()){
+                IssmDouble Sg_gauss=0;
+                TriaRef::GetInputValue(&Sg_gauss, Sg, gauss,P1Enum);
+                Sg_avg+=Sg_gauss*gauss->weight;
+                total_weight+=gauss->weight;
+        }
+        Sg_avg=Sg_avg/total_weight;
+        delete gauss;
+
+        *poceanaverage=Sg_avg;
+        *poceanarea=area;
+        return;
+
+}
+/*}}}*/
 #endif
 
 #ifdef _HAVE_DAKOTA_

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

@@ -161 +161 @@
                 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_SEALEVELCHANGE_
+                void       SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
+                void       SealevelchangeGeometry(IssmDouble* lat, IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
+                IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea);
+                IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea);
+                void       SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks);
+                void       SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks);
+                void       DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks);
+                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
                 IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks);
-                void    SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
-                void    SetSealevelMasks(SealevelMasks* masks);
-                void    SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
-                void    SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
-
-                IssmDouble    SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea);
-                IssmDouble    SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea);
-                void    SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks);
-                void    SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks);
-                void    DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks);
-                void    GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
-                void    SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix<IssmDouble>* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea);
-                IssmDouble SealevelchangeConvolution(IssmDouble* loads);
+                void       SetSealevelMasks(SealevelMasks* masks);
+                
+                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
+                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks);
+                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks);
+                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks);
+                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks);
                 #endif
                 /*}}}*/
                 /*Tria specific routines:{{{*/

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

@@ -4748 +4748 @@
 /*}}}*/
 #endif
 #ifdef _HAVE_SEALEVELCHANGE_
-void FemModel::SealevelchangeBarystatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition,SealevelMasks* masks) { /*{{{*/
-
-        /*serialized vectors:*/
-        IssmDouble  bslcice       = 0.;
-        IssmDouble  bslcice_cpu   = 0.;
-        IssmDouble  bslchydro       = 0.;
-        IssmDouble  bslchydro_cpu   = 0.;
-        IssmDouble  area      = 0.;
-        IssmDouble  oceanarea      = 0.;
-        IssmDouble  oceanarea_cpu  = 0.;
-        int bp_compute_fingerprints= 0;
-        bool isoceantransport=false;
-
-        Vector<IssmDouble>* bslcice_partition=NULL;
-        IssmDouble* bslcice_partition_serial=NULL;
-        IssmDouble* partitionice=NULL;
-        int npartice,nel;
-
-        Vector<IssmDouble>* bslchydro_partition=NULL;
-        IssmDouble* bslchydro_partition_serial=NULL;
-        IssmDouble* partitionhydro=NULL;
-        bool istws=0;
-        int nparthydro;
-                
-        int npartocean;
-        Vector<IssmDouble>* bslcocean_partition=NULL;
-        IssmDouble* bslcocean_partition_serial=NULL;
-        IssmDouble* partitionocean=NULL;
-
-   /*Initialize temporary vector that will be used to sum barystatic components
-    * on all local elements, prior to assembly:*/
-        int gsize = this->nodes->NumberOfDofs(GsetEnum);
-        IssmDouble* RSLgi=xNewZeroInit<IssmDouble>(gsize);
-        int* indices=xNew<int>(gsize);
-   for(int i=0;i<gsize;i++) indices[i]=i;
-
-        /*First, figure out the area of the ocean, which is needed to compute the barystatic component: */
-        int i = -1;
-        for(Object* & object : this->elements->objects){
-                i +=1;
-                Element* element = xDynamicCast<Element*>(object);
-                element->GetInputValue(&area,AreaEnum);
-                if (masks->isoceanin[i]) oceanarea_cpu += area;
-        }
-        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());
-        _assert_(oceanarea>0.);
-
-        /*Initialize partition vectors to retrieve barystatic contributions: */
-        this->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
-        if(npartice){
-                this->parameters->FindParam(&partitionice,&nel,NULL,SolidearthPartitionIceEnum);
-                bslcice_partition= new Vector<IssmDouble>(npartice);
-        }
-
-        this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
-        if(nparthydro){
-                this->parameters->FindParam(&partitionhydro,&nel,NULL,SolidearthPartitionHydroEnum);
-                bslchydro_partition= new Vector<IssmDouble>(nparthydro);
-        }
-
-        this->parameters->FindParam(&npartocean,SolidearthNpartOceanEnum);
-        if(npartocean){
-                this->parameters->FindParam(&partitionocean,&nel,NULL,SolidearthPartitionOceanEnum);
-                bslchydro_partition= new Vector<IssmDouble>(npartocean);
-        }
-        /*For later:
-        npartbarystatic=npartice;
-        if(nparthydro>npartbarystatic)npartbarystatic=nparthydro;
-        if(npartocean>npartbarystatic)npartbarystatic=npartocean;
-        bslc_partition=new Matrix(IssmDouble>(npartbarystatic,3);
-
-        bslc_cpu[0]=0; bslc_cpu[1]=0; bslc_cpu[2]=0;
-        for(Object* & object : this->elements->objects){
-                Element* element = xDynamicCast<Element*>(object);
-                element->SealevelchangeBarystaticLoads(&bslc_cpu[0], localloads,masks, bslcice_partition,partitionice,oceanarea);
-        }
-        MPI Bcast localloads -> loads
-
-        for(Object* & object : this->elements->objects){
-                Element* element = xDynamicCast<Element*>(object);
-                element->SealevelchangeConvolution(loads);
-        }
-        */
-
-
-
-
-
-        /*Call the barystatic sea level change core for ice : */
-        bslcice_cpu=0;
-        for(Object* & object : this->elements->objects){
-                Element* element = xDynamicCast<Element*>(object);
-                bslcice_cpu+=element->SealevelchangeBarystaticIce(RSLgi,masks, bslcice_partition,partitionice,oceanarea);
-        }
-
-        /*Call the barystatic sea level change core for hydro: */
-        bslchydro_cpu=0; //make sure to initialize this, so we have a total barystatic contribution computed at 0.
-        this->parameters->FindParam(&istws,TransientIshydrologyEnum);
-        if(istws){
-                for(int i=0;i<elements->Size();i++){
-                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        bslchydro_cpu+=element->SealevelchangeBarystaticHydro(RSLgi,masks, bslchydro_partition,partitionhydro,oceanarea);
-                }
-        }
-
-        /*Call the barystatic sea level change core for bottom pressures: */
-        this->parameters->FindParam(&bp_compute_fingerprints,SolidearthSettingsComputeBpGrdEnum);
-        this->parameters->FindParam(&isoceantransport,TransientIsoceantransportEnum);
-        if(bp_compute_fingerprints && isoceantransport){
-                for(int i=0;i<elements->Size();i++){
-                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
-                        element->SealevelchangeBarystaticBottomPressure(RSLgi,masks);
-                }
-        }
-
-        /*Plug values once and assemble: */
-        pRSLgi->SetValues(gsize,indices,RSLgi,ADD_VAL);
-        pRSLgi->Assemble();
-
-        /*Sum all barystatic components from all cpus:*/
-        ISSM_MPI_Reduce (&bslcice_cpu,&bslcice,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&bslcice,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        _assert_(!xIsNan<IssmDouble>(bslcice));
-
-        ISSM_MPI_Reduce (&bslchydro_cpu,&bslchydro,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
-        ISSM_MPI_Bcast(&bslchydro,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
-        _assert_(!xIsNan<IssmDouble>(bslchydro));
-
-        /*Take care of partition vectors:*/
-        if(bslcice_partition){
-                bslcice_partition->Assemble();
-                bslcice_partition_serial=bslcice_partition->ToMPISerial();
-        }
-        if(bslchydro_partition){
-                bslchydro_partition->Assemble();
-                bslchydro_partition_serial=bslchydro_partition->ToMPISerial();
-        }
-
-
-        /*Free ressources:*/
-        xDelete<int>(indices);
-        xDelete<IssmDouble>(RSLgi);
-        if(bslchydro_partition)delete bslchydro_partition;
-        if(bslcice_partition)delete bslcice_partition;
-        if(partitionhydro)xDelete<IssmDouble>(partitionhydro);
-        if(partitionice)xDelete<IssmDouble>(partitionice);
-
-        /*Assign output pointers:*/
-        *poceanarea = oceanarea;
-        *pbslcice  = bslcice;
-        *pbslchydro  = bslchydro;
-        *pbslc=bslchydro+bslcice;
-        *pbslcice_partition=bslcice_partition_serial;
-        *pbslchydro_partition=bslchydro_partition_serial;
-
-}
-/*}}}*/
 void FemModel::SealevelchangeSal(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old,  SealevelMasks* masks, bool verboseconvolution){/*{{{*/
 
         /*serialized vectors:*/

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

@@ -18 +18 @@
 #include "./Massfluxatgate.h"
 #include "./Misfit.h"
 #include "./SealevelMasks.h"
+#include "./BarystaticContributions.h"
 #include "./Nodalvalue.h"
 #include "./Numberedcostfunction.h"
 #include "./Cfsurfacesquare.h"

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

@@ -80 +80 @@
 void WriteLockFile(char* filename);
 void ResetBoundaryConditions(FemModel* femmodel, int analysis_type);
 void PrintBanner(void);
-void TransferForcing(FemModel* femmodel,int forcingenum);
-void TransferSealevel(FemModel* femmodel,int forcingenum);
 void EarthMassTransport(FemModel* femmodel);
-void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
 
 //solution configuration
 void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype);

../trunk-jpl/src/c/cores/sealevelchange_core.cpp

@@ -11 +11 @@
 #include "../shared/shared.h"
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
+/*support routines local definitions:{{{*/
+void TransferForcing(FemModel* femmodel,int forcingenum);
+void TransferSealevel(FemModel* femmodel,int forcingenum);
+void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
+IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea);
+/*}}}*/
 
 /*main cores:*/
 void sealevelchange_core(FemModel* femmodel){ /*{{{*/
@@ -383 +389 @@
         }
 }; /*}}}*/
 
+void grd_core_optim(FemModel* femmodel,SealevelMasks* masks) { /*{{{*/
+
+        /*variables:{{{*/
+        int nel;
+        BarystaticContributions* barycontrib=NULL;
+        GenericParam<BarystaticContributions*>* barycontribparam=NULL;
+        IssmDouble               barystatic;
+        
+        Vector<IssmDouble>*    loads=NULL;
+        IssmDouble*            allloads=NULL; 
+        Vector<IssmDouble>*    sealevelloads=NULL;
+        Vector<IssmDouble>*    oldsealevelloads=NULL;
+        IssmDouble             sealevelloadsaverage;
+        IssmDouble*            allsealevelloads=NULL;
+        Vector<IssmDouble>*    oceanareas=NULL;
+        IssmDouble             oceanarea;
+        bool                   scaleoceanarea=false;
+        IssmDouble             rho_water;
+
+        IssmDouble           eps_rel;
+        IssmDouble           eps_abs;
+        int                  step;
+        IssmDouble           time; 
+        
+        IssmDouble cumbslc;
+        IssmDouble cumbslcice;
+        IssmDouble cumbslchydro;
+        /*}}}*/
+
+        /*retrieve parameters: */
+        femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
+        barycontribparam = xDynamicCast<GenericParam<BarystaticContributions*>*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum));
+        barycontrib=barycontribparam->GetParameterValue();
+        femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
+        femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
+        femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
+        
+        /*initialize matrices and vectors:*/
+        femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
+        loads=new Vector<IssmDouble>(nel);
+        sealevelloads=new Vector<IssmDouble>(nel);
+        oceanareas=new Vector<IssmDouble>(nel);
+        
+        /*buildup loads: */
+        for(Object* & object : femmodel->elements->objects){
+                Element* element = xDynamicCast<Element*>(object);
+                element->SealevelchangeBarystaticLoads(loads, barycontrib,masks); 
+        }
+
+        //Communicate loads from local to global: 
+        loads->Assemble(); allloads=loads->ToMPISerial();
+
+        /*convolve loads:*/
+        for(Object* & object : femmodel->elements->objects){
+                Element* element = xDynamicCast<Element*>(object);
+                element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,masks);
+        }
+
+        //Get ocean area: 
+        oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
+        if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
+
+        //Get sea level loads ocean average: 
+        sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
+        
+        //substract ocean average and barystatic contributionfrom sea level loads: 
+        barystatic=barycontrib->Total()/oceanarea/rho_water;
+        sealevelloads->Shift(-sealevelloadsaverage+barystatic);
+        allsealevelloads=sealevelloads->ToMPISerial();
+
+        bool converged=false;
+        for(;;){
+                        
+                oldsealevelloads=sealevelloads->Duplicate();
+
+                /*convolve load and sealevel loads on oceans:*/
+                for(Object* & object : femmodel->elements->objects){
+                        Element* element = xDynamicCast<Element*>(object);
+                        element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,masks);
+                }
+                sealevelloads->Assemble();
+                
+                //substract ocean average and barystatic contribution
+                sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
+                sealevelloads->Shift(-sealevelloadsaverage+barystatic);
+
+                //broadcast loads 
+                allsealevelloads=sealevelloads->ToMPISerial();
+
+                //convergence?
+                slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs);
+                if (converged)break;
+        }
+
+        /*cumulate barystatic contributions and save to results: */
+        barycontrib->Cumulate(femmodel->parameters);
+        barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
+
+}
+/*}}}*/
+
 //Geometry:
 SealevelMasks* sealevel_masks(FemModel* femmodel) {  /*{{{*/
 
@@ -1061 +1168 @@
         *pconverged=converged;
 
 } /*}}}*/
+IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea){ /*{{{*/
+        IssmDouble sealevelloadsaverage;        
+        Vector<IssmDouble>* sealevelloadsvolume=sealevelloads->Duplicate();
+        sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas);
+        sealevelloadsvolume->Sum(&sealevelloadsaverage);
+        delete sealevelloadsvolume; 
+        return sealevelloadsaverage/oceanarea;
+} /*}}}*/

../trunk-jpl/src/c/shared/Enum/Enum.vim

@@ -746 +746 @@
 syn keyword cConstant BslcEnum
 syn keyword cConstant BslcIceEnum
 syn keyword cConstant BslcHydroEnum
+syn keyword cConstant BslcOceanEnum
 syn keyword cConstant BslcRateEnum
 syn keyword cConstant GmtslcEnum
 syn keyword cConstant SealevelRSLBarystaticEnum
@@ -1070 +1071 @@
 syn keyword cConstant BalancevelocitySolutionEnum
 syn keyword cConstant BasalforcingsIsmip6Enum
 syn keyword cConstant BasalforcingsPicoEnum
+syn keyword cConstant BarystaticContributionsEnum
 syn keyword cConstant BeckmannGoosseFloatingMeltRateEnum
 syn keyword cConstant BedSlopeSolutionEnum
 syn keyword cConstant BoolExternalResultEnum
@@ -1432 +1434 @@
 syn keyword cType AmrBamg
 syn keyword cType AmrNeopz
 syn keyword cType ArrayInput
+syn keyword cType BarystaticContributions
 syn keyword cType BoolInput
 syn keyword cType BoolParam
 syn keyword cType Cfdragcoeffabsgrad

../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -742 +742 @@
         BslcEnum,
         BslcIceEnum,
         BslcHydroEnum,
+        BslcOceanEnum,
         BslcRateEnum,
         GmtslcEnum,
         SealevelRSLBarystaticEnum,
@@ -1069 +1070 @@
         BalancevelocitySolutionEnum,
         BasalforcingsIsmip6Enum,
         BasalforcingsPicoEnum,
+        BarystaticContributionsEnum,
         BeckmannGoosseFloatingMeltRateEnum,
         BedSlopeSolutionEnum,
         BoolExternalResultEnum,

../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -748 +748 @@
                 case BslcEnum : return "Bslc";
                 case BslcIceEnum : return "BslcIce";
                 case BslcHydroEnum : return "BslcHydro";
+                case BslcOceanEnum : return "BslcOcean";
                 case BslcRateEnum : return "BslcRate";
                 case GmtslcEnum : return "Gmtslc";
                 case SealevelRSLBarystaticEnum : return "SealevelRSLBarystatic";
@@ -1072 +1073 @@
                 case BalancevelocitySolutionEnum : return "BalancevelocitySolution";
                 case BasalforcingsIsmip6Enum : return "BasalforcingsIsmip6";
                 case BasalforcingsPicoEnum : return "BasalforcingsPico";
+                case BarystaticContributionsEnum : return "BarystaticContributions";
                 case BeckmannGoosseFloatingMeltRateEnum : return "BeckmannGoosseFloatingMeltRate";
                 case BedSlopeSolutionEnum : return "BedSlopeSolution";
                 case BoolExternalResultEnum : return "BoolExternalResult";

../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -766 +766 @@
               else if (strcmp(name,"Bslc")==0) return BslcEnum;
               else if (strcmp(name,"BslcIce")==0) return BslcIceEnum;
               else if (strcmp(name,"BslcHydro")==0) return BslcHydroEnum;
+              else if (strcmp(name,"BslcOcean")==0) return BslcOceanEnum;
               else if (strcmp(name,"BslcRate")==0) return BslcRateEnum;
               else if (strcmp(name,"Gmtslc")==0) return GmtslcEnum;
               else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum;
@@ -873 +874 @@
               else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
               else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
               else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
-              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
+              if (strcmp(name,"SmbT")==0) return SmbTEnum;
+              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
               else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
               else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
               else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
@@ -996 +997 @@
               else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
               else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
               else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
-              else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
+              if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
+              else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
               else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum;
               else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum;
               else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum;
@@ -1096 +1097 @@
               else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
               else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum;
               else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum;
+              else if (strcmp(name,"BarystaticContributions")==0) return BarystaticContributionsEnum;
               else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
               else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
               else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
@@ -1118 +1120 @@
               else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum;
               else if (strcmp(name,"Closed")==0) return ClosedEnum;
               else if (strcmp(name,"Colinear")==0) return ColinearEnum;
-              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
-              else if (strcmp(name,"Contact")==0) return ContactEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"Contour")==0) return ContourEnum;
+              if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
+              else if (strcmp(name,"Contact")==0) return ContactEnum;
+              else if (strcmp(name,"Contour")==0) return ContourEnum;
               else if (strcmp(name,"Contours")==0) return ContoursEnum;
               else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
               else if (strcmp(name,"ControlInputGrad")==0) return ControlInputGradEnum;
@@ -1241 +1243 @@
               else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum;
               else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
               else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
-              else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
-              else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
+              if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
+              else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
+              else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
               else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
               else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum;
               else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
@@ -1364 +1366 @@
               else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
               else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum;
               else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;
-              else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
-              else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
          else stage=12;
    }
    if(stage==12){
-              if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
+              if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
+              else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
+              else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
               else if (strcmp(name,"Scaled")==0) return ScaledEnum;
               else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
               else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;

../trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h

@@ -323 +323 @@
 
                 }
                 /*}}}*/
+                void Sum(doubletype* pvalue){/*{{{*/
+
+                        doubletype value=0;
+                        int i;
+                        for(i=0;i<this->M;i++)value+=this->vector[i];
+                        *pvalue=value;
+                }
+                /*}}}*/
                 
 };
 #endif //#ifndef _ISSM_SEQ_VEC_H_

../trunk-jpl/src/c/toolkits/objects/Matrix.h

@@ -281 +281 @@
                         return output;
                 }
                 /*}}}*/
+                doubletype* ToMPISerial(void){/*{{{*/
+
+                        doubletype* output=NULL;
+
+                        if(type==PetscMatType){
+                                #ifdef _HAVE_PETSC_
+                                output=this->pmatrix->ToMPISerial();
+                                #endif
+                        }
+                        else{
+                                _error_("not implemented yet!");
+                        }
+
+                        return output;
+                }
+                /*}}}*/
                 void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){/*{{{*/
 
                         if(type==PetscMatType){
@@ -306 +322 @@
 
                 }
                 /*}}}*/
-                /*
-                * sets all values to 0 but keeps the structure of a sparse matrix
-                */
                 void SetZero(void) {/*{{{*/
+                        // sets all values to 0 but keeps the structure of a sparse matrix
                         if(type==PetscMatType){
                                 #ifdef _HAVE_PETSC_
                                 this->pmatrix->SetZero();

../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp

@@ -248 +248 @@
 
 }
 /*}}}*/
+void PetscVec::Sum(IssmDouble* pvalue){/*{{{*/
+
+        _assert_(this->vector);
+        VecSum(this->vector,pvalue);
+
+}
+/*}}}*/
 IssmDouble PetscVec::Dot(PetscVec* input){/*{{{*/
 
         IssmDouble dot;

../trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h

@@ -31 +31 @@
 void PetscOptionsDetermineSolverType(int* psolver_type);
 void MatMultPatch(Mat A,Vec X, Vec AX,ISSM_MPI_Comm comm);
 void MatToSerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
+void MatToMPISerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm);
 Vec  SerialToVec(double* vector,int vector_size);
 InsertMode ISSMToPetscInsertMode(InsMode mode);
 NormType ISSMToPetscNormMode(NormMode mode);

../trunk-jpl/test/NightlyRun/test2004.m

@@ -133 +133 @@
 
                 disp('      reading bedrock');
                 md.geometry.bed=-ones(md.mesh.numberofvertices,1);
+                md.geometry.base=md.geometry.bed;
+                md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1);
+                md.geometry.surface=md.geometry.bed+md.geometry.thickness;
+
         end % }}}
         %Slc: {{{
         if bas.iscontinentany('antarctica'),
@@ -166 +170 @@
                         md.masstransport.spcthickness=mean(delHAIS(md.mesh.elements),2)/100;
                 end
 
-                md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1);
+                md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
 
-                md.dsl.global_average_thermosteric_sea_level_change=[0;0];
-                md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1);
-                md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
+                md.dsl.global_average_thermosteric_sea_level=[0;0];
+                md.dsl.sea_surface_height_above_geoid=zeros(md.mesh.numberofvertices+1,1);
+                md.dsl.sea_water_pressure_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
 
         end %}}}
         % material properties: {{{
@@ -312 +316 @@
         %geometry:  {{{
         di=md.materials.rho_ice/md.materials.rho_water;
         md.geometry.bed=-ones(md.mesh.numberofvertices,1);
+        md.geometry.base=md.geometry.bed;
+        md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1);
+        md.geometry.surface=md.geometry.bed+md.geometry.thickness;
         % }}}
         %materials:  {{{
         md.materials=materials('hydro');
@@ -345 +352 @@
 sl.transfer('mask.ice_levelset');
 sl.transfer('mask.ocean_levelset');
 sl.transfer('geometry.bed');
+sl.transfer('geometry.surface');
+sl.transfer('geometry.thickness');
+sl.transfer('geometry.base');
 sl.transfer('mesh.lat');
 sl.transfer('mesh.long');
 sl.transfer('masstransport.spcthickness');
@@ -399 +409 @@
 md.transient.ismasstransport=1;
 md.transient.isslc=1;
 
+%Initializations: 
+md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
+md.initialization.bottompressure=zeros(md.mesh.numberofvertices,1);
+md.initialization.dsl=zeros(md.mesh.numberofvertices,1);
+md.initialization.str=0;
+md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
 
 %max number of iterations reverted back to 10 (i.e. the original default value)
 md.solidearth.settings.maxiter=10;