Changeset 26222


Ignore:
Timestamp:
04/29/21 11:55:48 (4 years ago)
Author:
Eric.Larour
Message:

CHG: Mathieu, here it is, just one year too late :) or one year early depending on how you see it.

Dakota: pre core capability, for things called only once prior to sampling, which could
not be called in the ModelProcessorx.

SealevelGeometry now used to optimize structures passed betweens sea level cores.

SealevelMasks: gone.

Mask logic now moved to the geometry core in the sea level solution. We now compute
barycentres of loads, for elements that can be interesected by up to two levelsets.
Allows for loading differentailly for hydro, ice, bp and ocean loads. Essentially
a sub-element parameterization of sea level solver.

Location:
issm/trunk-jpl
Files:
5 added
1 deleted
25 edited

Legend:

Unmodified
Added
Removed

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

    r26156 r26222  
    249249        ./cores/ResetBoundaryConditions.cpp \
    250250        ./cores/WrapperCorePointerFromSolutionEnum.cpp \
     251        ./cores/WrapperPreCorePointerFromSolutionEnum.cpp \
    251252        ./cores/CorePointerFromSolutionEnum.cpp \
    252253        ./cores/ad_core.cpp \
    … …  
    557558issm_sources += \
    558559        ./cores/sealevelchange_core.cpp \
    559         ./analyses/SealevelchangeAnalysis.cpp
     560        ./analyses/SealevelchangeAnalysis.cpp\
     561        ./classes/SealevelGeometry.cpp
    560562
    561563#gia ivins physics (only if have fortran)

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

    r26165 r26222  
    143143        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.compute_bp_grd",SolidearthSettingsComputeBpGrdEnum));
    144144        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum));
    145         parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum));
    146145        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.cross_section_shape",SolidearthSettingsCrossSectionShapeEnum));
    147146        parameters->AddObject(new DoubleParam(CumBslcEnum,0.0));
    … …  
    429428
    430429        /*run sea level change core geometry only once, after the Model Processor is done:*/
    431         sealevelchange_geometry(femmodel);
     430        sealevelchange_initialgeometry(femmodel);
    432431
    433432}/*}}}*/

  • issm/trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp

    r25947 r26222  
    5555                bool       control_analysis         = false;
    5656                void     (*solutioncore)(FemModel*) = NULL;
     57                void     (*solutionprecore)(FemModel*) = NULL;
    5758                bool       nodakotacore             = true;
    5859
    … …  
    9495                responses=xNewZeroInit<IssmDouble>(numFns);
    9596
    96                 /*Hack:*/
    97                 for(int i=0;i<femmodel_init->nummodels;i++) if(femmodel_init->analysis_type_list[i]==SealevelchangeAnalysisEnum) sealevelchange_geometry(femmodel_init);
     97                /*Launch cores that are not used during the uncertainty quantification: */
     98                WrapperPreCorePointerFromSolutionEnum(&solutionprecore,femmodel_init->parameters,solution_type);
     99                if(solutionprecore)solutionprecore(femmodel_init);
    98100
    99101                /*Make a copy of femmodel, so we start this new evaluation run for this specific sample with a brand

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

    r26194 r26222  
    20342034bool       Element::IsIceOnlyInElement(){/*{{{*/
    20352035        Input* input=this->GetInput(MaskIceLevelsetEnum); _assert_(input);
    2036         return (input->GetInputMax()<0.);
     2036        return (input->GetInputMax()<=0.);
    20372037}
    20382038/*}}}*/
    … …  
    20452045        Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input);
    20462046        return (input->GetInputMin()<0.);
     2047}
     2048/*}}}*/
     2049bool       Element::IsOceanOnlyInElement(){/*{{{*/
     2050        Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input);
     2051        return (input->GetInputMax()<=0.);
    20472052}
    20482053/*}}}*/

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

    r26167 r26222  
    3232class IoModel;
    3333class SealevelMasks;
     34class SealevelGeometry;
    3435class Gauss;
    3536class ElementVector;
    … …  
    146147                bool               IsIceOnlyInElement();
    147148                bool               IsOceanInElement();
     149                bool               IsOceanOnlyInElement();
    148150                bool               IsLandInElement();
    149151                void               Ismip6FloatingiceMeltingRate();
    … …  
    240242                virtual void       ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
    241243                virtual void       ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze, Vector<IssmDouble>* vareae, bool spherical=false)=0;
     244                virtual void       ElementCoordinates(Vector<IssmDouble>* vlonge,Vector<IssmDouble>* vlate,Vector<IssmDouble>* vareae)=0;
    242245                virtual int        FiniteElement(void)=0;
    243246                virtual IssmDouble FloatingArea(bool scaled)=0;
    … …  
    248251                virtual void       GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0;
    249252                virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0;
    250                 virtual void        GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl)=0;
    251                 virtual IssmDouble  GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl)=0;
     253                virtual void        GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl)=0;
     254                virtual void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelsetenum)=0;
     255                virtual void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum)=0;
     256                virtual void        GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius)=0;
     257
    252258                virtual IssmDouble GetIcefrontArea(){_error_("not implemented");};
    253259                virtual void       GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0;
    … …  
    379385                #endif
    380386                #ifdef _HAVE_SEALEVELCHANGE_
    381                 virtual void          LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){_error_("not implemented");};
    382                 virtual void          SetSealevelMasks(SealevelMasks* masks)=0;
    383387                virtual IssmDouble    GetArea3D(void)=0;
    384388                virtual IssmDouble    GetAreaSpherical(void)=0;
    385                 virtual void          SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads)=0;
     389                virtual IssmDouble    GetTriangleAreaSpherical(IssmDouble xyz_list[3][3])=0;
    386390                virtual void          GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0;
    387391
    388                 virtual void          SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
    389                 virtual void          SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0;
    390                 virtual void          SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector)=0;
    391                 virtual void          SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector)=0;
    392                 virtual void          SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks)=0;
     392                virtual void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom)=0;
     393                virtual void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom)=0;
     394                virtual void       SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom)=0;
     395                virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
     396                virtual void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
     397                virtual void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0;
     398                virtual void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom)=0;
     399                virtual void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
     400                virtual void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
     401                virtual void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom)=0;
    393402                #endif
    394403

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

    r26167 r26222  
    7171                void           ElementResponse(IssmDouble* presponse,int response_enum);
    7272                void           ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,Vector<IssmDouble>* vareae,bool spherical=false){_error_("not implemented yet");};
     73                void           ElementCoordinates(Vector<IssmDouble>* vlonge,Vector<IssmDouble>* vlate,Vector<IssmDouble>* vareae){_error_("not implemented yet");};
    7374                void           ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    7475                int            FiniteElement(void);
    … …  
    7879                IssmDouble     GetAreaSpherical(void){_error_("not implemented yet!");};
    7980                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
     81                IssmDouble     GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]){_error_("not implemented yet");};
    8082                Element*       GetBasalElement(void);
    8183                Penta*         GetBasalPenta(void);
    … …  
    8385                void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
    8486                IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
    85                 void           GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");};
    86                 IssmDouble     GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");};
     87                void           GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){_error_("not implemented yet");};
     88                void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelsetenum){_error_("not implemented yet");};
     89                void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum){_error_("not implemented yet");};
     90                void        GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius){_error_("not implemented yet");};
    8791                IssmDouble              GetIcefrontArea();
    8892                void           GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
    … …  
    218222                #endif
    219223                #ifdef _HAVE_SEALEVELCHANGE_
    220                 void       SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    221                 void       SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
    222224                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    223 
    224                 void       SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    225                 void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    226                 void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");};
    227                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");};
    228                 void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
     225                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
     226                void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     227                void       SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
     228                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
     229                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
     230                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
     231                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
     232                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     233                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     234                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    229235                #endif
    230236

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

    r26167 r26222  
    5555                void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
    5656                void        ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,Vector<IssmDouble>* vareae,bool spherical=false){_error_("not implemented yet");};
     57                void        ElementCoordinates(Vector<IssmDouble>* vlonge,Vector<IssmDouble>* vlate,Vector<IssmDouble>* vareae){_error_("not implemented yet");};
    5758                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
    5859                int         FiniteElement(void);
    … …  
    6364                void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
    6465                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
    65                 void        GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");};
    66                 IssmDouble  GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");};
    67                 void               GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
     66                void        GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){_error_("not implemented yet");};
     67                void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelsetenum){_error_("not implemented yet");};
     68                void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelset1enum, int levelset2enum){_error_("not implemented yet");};
     69                void        GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius){_error_("not implemented yet");};
     70                void            GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
    6871                Input*     GetInput(int enumtype);
    6972                Input*     GetInput(int enumtype,IssmDouble time);
    … …  
    167170                IssmDouble     GetArea3D(void){_error_("not implemented yet!");};
    168171                IssmDouble     GetAreaSpherical(void){_error_("not implemented yet!");};
     172                IssmDouble GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]){_error_("not implemented yet");};
    169173
    170174#ifdef _HAVE_ESA_
    … …  
    173177#endif
    174178#ifdef _HAVE_SEALEVELCHANGE_
    175                 void       SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
    176                 void       SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    177179                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    178                 void       SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    179                 void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    180                 void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");};
    181                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");};
    182                 void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
    183 
     180                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
     181                void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     182                void       SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
     183                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
     184                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
     185                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
     186                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
     187                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     188                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     189                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    184190#endif
    185191

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

    r26167 r26222  
    5454                void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
    5555                void        ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze,Vector<IssmDouble>* vareae,bool spherical=false){_error_("not implemented yet");};
     56                void        ElementCoordinates(Vector<IssmDouble>* vlonge,Vector<IssmDouble>* vlate,Vector<IssmDouble>* vareae){_error_("not implemented yet");};
    5657                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    5758                void        FaceOnBaseIndices(int* pindex1,int* pindex2,int* pindex3);
    … …  
    6364                IssmDouble  GetArea3D(void){_error_("not implemented yet!");};
    6465                IssmDouble  GetAreaSpherical(void){_error_("not implemented yet!");};
     66                IssmDouble  GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]){_error_("not implemented yet");};
    6567                Element*    GetBasalElement(void){_error_("not implemented yet");};
    6668                int         GetElementType(void);
    6769                void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
    6870                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
    69                 void        GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){_error_("not implemented yet");};
    70                 IssmDouble  GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){_error_("not implemented yet");};
     71                void        GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){_error_("not implemented yet");};
     72                void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelsetenum){_error_("not implemented yet");};
     73                void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelset1enum, int levelset2enum){_error_("not implemented yet");};
     74                void        GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius){_error_("not implemented yet");};
    7175                void               GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
    7276                Input*     GetInput(int enumtype);
    … …  
    180184#endif
    181185#ifdef _HAVE_SEALEVELCHANGE_
    182                 void       SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    183                 void       SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
    184186                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    185                 void       SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    186                 void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    187                 void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector){_error_("not implemented yet");};
    188                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){_error_("not implemented yet");};
    189                 void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){_error_("not implemented yet");};
    190 
     187                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
     188                void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
     189                void       SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
     190                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
     191                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
     192                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};
     193                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){_error_("not implemented yet");};
     194                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     195                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
     196                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    191197#endif
    192198

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

    r26165 r26222  
    13121312}
    13131313/*}}}*/
     1314void       Tria::ElementCoordinates(Vector<IssmDouble>* vlonge,Vector<IssmDouble>* vlate,Vector<IssmDouble>* vareae){ /*{{{*/
     1315
     1316        IssmDouble planetradius;
     1317
     1318        /*Look for x,y,z coordinates:*/
     1319        IssmDouble xyz_list[NUMVERTICES][3];
     1320        ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
     1321        IssmDouble area,xe,ye,ze,late,longe;
     1322
     1323        /*Find centroid:*/
     1324        xe=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
     1325        ye=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
     1326        ze=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
     1327
     1328        late= asin(ze/sqrt(pow(xe,2.0)+pow(ye,2.0)+pow(ze,2.0)));
     1329        longe= atan2(ye,xe);
     1330        area=this->GetAreaSpherical();
     1331
     1332        vlonge->SetValue(this->sid,longe,INS_VAL);
     1333        vlate->SetValue(this->sid,late,INS_VAL);
     1334        vareae->SetValue(this->sid,area,INS_VAL);
     1335               
     1336        return;
     1337}
     1338/*}}}*/
    13141339void       Tria::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
    13151340
    … …  
    17231748}
    17241749/*}}}*/
    1725 IssmDouble Tria::GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl){/*{{{*/
    1726        
    1727         /*Computeportion of the element that is on a levelset: */
    1728         bool              mainlyfloating = true;
    1729         int               domaintype,index1,index2;
    1730         const IssmPDouble epsilon        = 1.e-15;
    1731         IssmDouble        phi,s1,s2;
    1732 
    1733         /*Recover parameters and values*/
    1734         parameters->FindParam(&domaintype,DomainTypeEnum);
    1735 
    1736         /*Be sure that values are not zero*/
    1737         if(gl[0]==0.) gl[0]=gl[0]+epsilon;
    1738         if(gl[1]==0.) gl[1]=gl[1]+epsilon;
    1739         if(gl[2]==0.) gl[2]=gl[2]+epsilon;
    1740 
    1741         if(domaintype==Domain2DverticalEnum){
    1742                 this->EdgeOnBaseIndices(&index1,&index2);
    1743                 if(gl[index1]>0 && gl[index2]>0) phi=1; // All grounded
    1744                 else if(gl[index1]<0 && gl[index2]<0) phi=0; // All floating
    1745                 else if(gl[index1]<0 && gl[index2]>0){ //index2 grounded
    1746                         phi=1./(1.-gl[index1]/gl[index2]);
    1747                 }
    1748                 else if(gl[index2]<0 && gl[index1]>0){ //index1 grounded
    1749                         phi=1./(1.-gl[index2]/gl[index1]);
    1750                 }
    1751 
    1752         }
    1753         else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
    1754                 /*Check that not all nodes are grounded or floating*/
    1755                 if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
    1756                         phi=1;
    1757                 }
    1758                 else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
    1759                         phi=0;
    1760                 }
    1761                 else{
    1762                         /*Figure out if two nodes are floating or grounded*/
    1763                         if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=false;
    1764 
    1765                         if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
    1766                                 s1=gl[2]/(gl[2]-gl[1]);
    1767                                 s2=gl[2]/(gl[2]-gl[0]);
    1768                         }
    1769                         else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
    1770                                 s1=gl[0]/(gl[0]-gl[1]);
    1771                                 s2=gl[0]/(gl[0]-gl[2]);
    1772                         }
    1773                         else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
    1774                                 s1=gl[1]/(gl[1]-gl[0]);
    1775                                 s2=gl[1]/(gl[1]-gl[2]);
    1776                         }
    1777                         else _error_("case not possible");
    1778                         if(mainlyfloating){
    1779                                 phi = (1-s1*s2);
    1780                         }
    1781                         else{
    1782                                 phi = s1*s2;
    1783                         }
    1784                 }
    1785         }
    1786         else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
    1787 
    1788         _assert_(phi<=1. && phi>=0.);
    1789         return 1-phi;
    1790 }
    1791 /*}}}*/
    1792 void       Tria::GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl){/*{{{*/
     1750void       Tria::GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){/*{{{*/
    17931751       
    17941752        /*Computeportion of the element that is grounded*/
    1795         bool               negative=false;
     1753        bool               trapezeisnegative=true; //default value
    17961754        int                point;
    17971755        const IssmPDouble  epsilon= 1.e-15;
    1798         IssmDouble         f1,f2;
     1756        IssmDouble         f1,f2,phi;
     1757
     1758        /*Weights: */
     1759        Gauss* gauss=NULL;
     1760        IssmDouble loadweights_g[NUMVERTICES];
     1761        IssmDouble total_weight=0;
     1762       
     1763        _assert_(!xIsNan<IssmDouble>(gl[0]));
     1764        _assert_(!xIsNan<IssmDouble>(gl[1]));
     1765        _assert_(!xIsNan<IssmDouble>(gl[2]));
    17991766
    18001767        /*Be sure that values are not zero*/
    … …  
    18151782        }
    18161783        else{
    1817                 if(gl[0]*gl[1]*gl[2]<0) negative=true;
     1784                if(gl[0]*gl[1]*gl[2]<0) trapezeisnegative=false; //no matter what configuration, there has to be two positive vertices, which means the trapeze is positive.
    18181785
    18191786                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
    … …  
    18341801                else _error_("case not possible");
    18351802        }
    1836         *point1=point;
    1837         *fraction1=f1;
    1838         *fraction2=f2;
    1839         *pmainlynegative=negative;
    1840 }
    1841 /*}}}*/
     1803        if(trapezeisnegative) phi=1-f1*f2;
     1804        else phi=f1*f2;
     1805       
     1806
     1807        /*Compute weights:*/
     1808        gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
     1809       
     1810        total_weight=0;
     1811        for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
     1812        while(gauss->next()){
     1813                TriaRef::GetNodalFunctions(&loadweights_g[0], gauss,P1Enum);
     1814                for(int i=0;i<NUMVERTICES;i++)weights[i]+=loadweights_g[i]*gauss->weight;
     1815                total_weight+=gauss->weight;
     1816        }
     1817        //normalize to phi.
     1818        if(total_weight)for(int i=0;i<NUMVERTICES;i++)weights[i]/=total_weight/phi;
     1819        else for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
     1820
     1821        /*free ressources:*/
     1822        delete gauss;
     1823
     1824        /*Assign output pointers:*/
     1825        *pphi=phi;
     1826        *ppoint1=point;
     1827        *pfraction1=f1;
     1828        *pfraction2=f2;
     1829        *ptrapezeisnegative=trapezeisnegative;
     1830}
     1831/*}}}*/
     1832IssmDouble Tria::GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]){/*{{{*/
     1833       
     1834        IssmDouble x1,y1,z1,x2,y2,z2,x3,y3,z3;
     1835        IssmDouble arc12,arc23,arc31,semi_peri,excess;
     1836        IssmDouble lat1,lat2,lat3;
     1837        IssmDouble long1,long2,long3;
     1838        IssmDouble r1,r2,r3;
     1839
     1840        /*retrieve x,y,z coordinates: */
     1841        x1=xyz_list[0][0]; y1=xyz_list[0][1]; z1=xyz_list[0][2];
     1842        x2=xyz_list[1][0]; y2=xyz_list[1][1]; z2=xyz_list[1][2];
     1843        x3=xyz_list[2][0]; y3=xyz_list[2][1]; z3=xyz_list[2][2];
     1844
     1845        /*Build lat,long, r:*/
     1846        r1=sqrt(pow(x1,2.0)+pow(y1,2.0)+pow(z1,2.0));
     1847        r2=sqrt(pow(x2,2.0)+pow(y2,2.0)+pow(z2,2.0));
     1848        r3=sqrt(pow(x3,2.0)+pow(y3,2.0)+pow(z3,2.0));
     1849
     1850        lat1=asin(z1/r1); long1=atan2(y1,x1);
     1851        lat2=asin(z2/r2); long2=atan2(y2,x2);
     1852        lat3=asin(z3/r3); long3=atan2(y3,x3);
     1853
     1854        /*compute great circle distance between vertices */
     1855        arc12=2.*asin(sqrt(pow(sin(0.5*(lat2-lat1)),2)+cos(lat1)*cos(lat2)*pow(sin(0.5*(long2-long1)),2)));
     1856        arc23=2.*asin(sqrt(pow(sin(0.5*(lat3-lat2)),2)+cos(lat2)*cos(lat3)*pow(sin(0.5*(long3-long2)),2)));
     1857        arc31=2.*asin(sqrt(pow(sin(0.5*(lat1-lat3)),2)+cos(lat3)*cos(lat1)*pow(sin(0.5*(long1-long3)),2)));
     1858
     1859        /*semi parameter */
     1860        semi_peri=(arc12+arc23+arc31)/2;
     1861
     1862        /*spherical excess */
     1863        excess=4.*atan(sqrt(tan(semi_peri/2)*tan((semi_peri-arc12)/2)*tan((semi_peri-arc23)/2)*tan((semi_peri-arc31)/2)));
     1864
     1865        /*area = excess*radius^2 */
     1866        return excess*pow((r1+r2+r3)/3,2);
     1867}
     1868/*}}}*/
     1869void       Tria:: GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius){ /*{{{*/
     1870
     1871        int i0,i1,i2;
     1872
     1873        IssmDouble xyz0[3][3];
     1874        IssmDouble barycenter[3]={0};
     1875        IssmDouble centroid[3]={0};
     1876
     1877        ::GetVerticesCoordinates(&xyz0[0][0],vertices,NUMVERTICES); // initial triangle
     1878
     1879        i0=point1;
     1880        i1=(point1+1)%3;
     1881        i2=(point1+2)%3;
     1882
     1883        //Barycenter of the subelement triangle:
     1884        for (int i=0;i<3;i++) barycenter[i]=xyz0[i0][i]*(3.0-fraction1-fraction2)/3.0 + fraction1/3.0*xyz0[i1][i] + fraction2/3.0*xyz0[i2][i];
     1885
     1886        if (istrapeze1){
     1887                centroid[0]=planetradius*cos(late*M_PI/180.0) * cos(longe*M_PI/180.0); //x
     1888                centroid[1]=planetradius*cos(late*M_PI/180.0) * sin(longe*M_PI/180.0);  //y
     1889                centroid[2]=planetradius*sin(late*M_PI/180.0);                                  //z
     1890
     1891                // centroid_el *area_el= barycenter_triangle * area_triangle + barycenter_trapeze * area_trapeze
     1892                // and phi_trapeze = area_trapeze/area_el = (1 - area_triangle/area_el)
     1893                // => barycenter_trapeze = (centroid_el - barycenter_triangle * (1-phi_trapeze) )/phi_trapeze
     1894                for (int i=0;i<3;i++) barycenter[i] =(centroid[i] -barycenter[i]*(1.0-phi))/phi;
     1895
     1896        }
     1897
     1898        //recompute planetradius from the barycenter onwards:
     1899        planetradius=sqrt( pow(barycenter[0],2.0)+ pow(barycenter[1],2.0)+ pow(barycenter[2],2.0));
     1900
     1901        *platbar=asin(barycenter[2]/planetradius)*180.0/M_PI;
     1902        *plongbar=atan2(barycenter[1],barycenter[0])*180.0/M_PI;
     1903
     1904} /*}}}*/
     1905void       Tria::GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelsetenum){ /*{{{*/
     1906
     1907        IssmDouble phi;
     1908        IssmDouble fraction1,fraction2;
     1909        bool istrapeze1; 
     1910        bool flip1=false;
     1911        int  point1;
     1912        IssmDouble levelset[NUMVERTICES];
     1913        IssmDouble planetradius;
     1914       
     1915        this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
     1916
     1917        //figure out if we are flipping the levelsets:
     1918        if(levelsetenum<0){
     1919                levelsetenum=-levelsetenum;
     1920                flip1=true;
     1921        }
     1922        //figure out area where we have loads
     1923        Element::GetInputListOnVertices(&levelset[0],levelsetenum);
     1924        if(flip1)for(int i=0;i<NUMVERTICES;i++)levelset[i]=-levelset[i];
     1925       
     1926        //compute sea level load weights
     1927        this->GetFractionGeometry(loadweights,&phi,&point1,&fraction1,&fraction2,&istrapeze1,levelset);
     1928
     1929        this->GetBarycenterFromLevelset(platbar,plongbar, phi, fraction1, fraction2, late, longe, point1,istrapeze1,planetradius);
     1930
     1931        /*assign output pointers:*/
     1932        *ploadarea=phi*area;
     1933
     1934} /*}}}*/
     1935void       Tria::GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelset1enum, int levelset2enum){ /*{{{*/
     1936
     1937
     1938        bool istrapeze1, istrapeze2;
     1939        IssmDouble phi1,phi2, f11,f12,f21,f22;
     1940        int point1, point2,  i0,i1,i2,j0,j1,j2;
     1941        IssmDouble weights1[3],weights2[3];
     1942        IssmDouble levelset1[3];
     1943        IssmDouble levelset2[3];
     1944
     1945        bool flip1=false;
     1946        bool flip2=false;
     1947
     1948        IssmDouble xyz0[3][3];
     1949        IssmDouble xyz1[3][3]={0};
     1950        IssmDouble xyz2[3][3]={0};
     1951        IssmDouble xyz3[3][3]={0};
     1952        IssmDouble xyz[8][3]={0};
     1953        IssmDouble f1o;
     1954        IssmDouble w[8][NUMVERTICES]={0};
     1955        IssmDouble areasub=0;
     1956        IssmDouble area1=0;
     1957        IssmDouble area2=0;
     1958        IssmDouble area3=0;
     1959
     1960        int tria0[3]={0,1,2};
     1961        int tria1[3]={-1};
     1962        int tria2[3]={-1};
     1963        int tria3[3]={-1};
     1964
     1965        IssmDouble w1[3][3]={0};
     1966        IssmDouble w2[3][3]={0};
     1967        IssmDouble w3[3][3]={0};
     1968
     1969        IssmDouble barycenter[3]={0};
     1970        IssmDouble planetradius;
     1971
     1972        //figure out if we are flipping the levelsets:
     1973        if(levelset1enum<0){
     1974                levelset1enum=-levelset1enum;
     1975                flip1=true;
     1976        }
     1977        if(levelset2enum<0){
     1978                levelset2enum=-levelset2enum;
     1979                flip2=true;
     1980        }
     1981
     1982        //recover levelsets:
     1983        Element::GetInputListOnVertices(&levelset1[0],levelset1enum);
     1984        if(flip1)for(int i=0;i<NUMVERTICES;i++)levelset1[i]=-levelset1[i];
     1985        Element::GetInputListOnVertices(&levelset2[0],levelset2enum);
     1986        if(flip2)for(int i=0;i<NUMVERTICES;i++)levelset2[i]=-levelset2[i];
     1987
     1988        //We want the fraction of the element where both levelsets are negative.
     1989        //Early return if either of them is >=0 on all vertices
     1990        if (levelset1[0]>=0 && levelset1[1]>=0 && levelset1[2]>=0) {
     1991                for (int i=0;i<NUMVERTICES;i++) loadweights[i]=0.0;
     1992                *ploadarea= 0.0;
     1993                *platbar=late; //just default to centroid of triangle
     1994                *plongbar=longe;
     1995                return;
     1996        }
     1997        if (levelset2[0]>=0 && levelset2[1]>=0 && levelset2[2]>=0) {
     1998                for (int i=0;i<NUMVERTICES;i++) loadweights[i]=0.0;
     1999                *ploadarea= 0.0;
     2000                *platbar=late; //just default to centroid of triangle
     2001                *plongbar=longe;
     2002                return;
     2003        }
     2004
     2005
     2006        //If everyone is negative, no need to calculate any fraction
     2007        if (levelset1[0]<=0 && levelset1[1]<=0 && levelset1[2]<=0 && levelset2[0]<=0 && levelset2[1]<=0 && levelset2[2]<=0) {
     2008                for (int i=0;i<NUMVERTICES;i++) loadweights[i]=1.0/NUMVERTICES;
     2009                *ploadarea= area;
     2010                *platbar=late;
     2011                *plongbar=longe;
     2012                return;
     2013        }
     2014       
     2015        /*recover planet radius:*/
     2016        this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
     2017
     2018
     2019        //If just one levelset is all negative, just take the partitioning of the other, no interaction between them
     2020        if (levelset1[0]<=0 && levelset1[1]<=0 && levelset1[2]<=0) {
     2021                this->GetFractionGeometry(loadweights,&phi2,&point2,&f21,&f22,&istrapeze2,levelset2);
     2022                this->GetBarycenterFromLevelset(platbar,plongbar, phi2, f21, f22, late, longe, point2,istrapeze2,planetradius);
     2023                *ploadarea=area*phi2;
     2024                return;
     2025        }
     2026        if (levelset2[0]<=0 && levelset2[1]<=0 && levelset2[2]<=0) {
     2027                this->GetFractionGeometry(loadweights,&phi1,&point1,&f11,&f12,&istrapeze1,levelset1);
     2028                this->GetBarycenterFromLevelset(platbar,plongbar, phi1, f11, f12, late, longe, point1,istrapeze1,planetradius);
     2029                *ploadarea=area*phi1;
     2030                return;
     2031        }
     2032
     2033
     2034        this->GetFractionGeometry(&weights1[0],&phi1,&point1,&f11,&f12,&istrapeze1,levelset1);
     2035        this->GetFractionGeometry(&weights2[0],&phi2,&point2,&f21,&f22,&istrapeze2,levelset2);
     2036
     2037
     2038        //Early return if levelsets are not independent
     2039        if (istrapeze1==istrapeze2 && point1==point2 && phi1==phi2){
     2040                //the two levelsets are redundant: levelset1 = positivescalar * levelset2
     2041                this->GetBarycenterFromLevelset(platbar,plongbar, phi1, f11, f12, late, longe, point1,istrapeze1,planetradius);
     2042                *ploadarea=area*phi1;
     2043                for (int i=0;i<NUMVERTICES;i++) loadweights[i]=weights1[i];
     2044                return;
     2045        }
     2046        if (istrapeze1!=istrapeze2 && point1==point2 && phi1==(1.0-phi2)){
     2047                //the two levelsets are incompatible: levelset1 = negativescalar * levelset2
     2048                *plongbar=longe;
     2049                *platbar=late;
     2050                *ploadarea=0.0;
     2051                for (int i=0;i<NUMVERTICES;i++) loadweights[i]=0.0;
     2052                return;
     2053        }
     2054
     2055       
     2056        ::GetVerticesCoordinates(&xyz0[0][0],vertices,NUMVERTICES); // initial triangle
     2057
     2058        i0=point1;
     2059        i1=(point1+1)%3;
     2060        i2=(point1+2)%3;
     2061
     2062        j0=point2;
     2063        j1=(point2+1)%3;
     2064        j2=(point2+2)%3;
     2065
     2066        //f1o: fraction of segment {point_f11 -> point_f12} where the levelsets intersect (only used when f1o>=0 and f1o<=1)
     2067        if(point2==i0) f1o= f22*(f11-f21)/(f11*f22-f12*f21);
     2068        else if(point2==i1) f1o=f21*(1.0-f22-f11)/(f21*(f12-f11)-f12*f22);
     2069        else f1o= (f22*(1.0-f21-f11)+f21*f11)/(f22*(f12-f11) +f21*f11);
     2070
     2071
     2072        //interpolant weights of each point. Any field F[0,1,2] provided at the original vertices [0,1,2] will be equal on point k to sum_i (F[i] * w[k][i])
     2073        w[0][0]=1;
     2074        w[1][1]=1;
     2075        w[2][2]=1;
     2076        w[3][i0]=1.0-f11; w[3][i1]=f11;
     2077        w[4][i0]=1.0-f12; w[4][i2]=f12;
     2078        w[5][j0]=1.0-f21; w[5][j1]=f21;
     2079        w[6][j0]=1.0-f22; w[6][j2]=f22;
     2080        for (int j=0;j<3;j++) w[7][j]=w[3][j]*(1.0-f1o) + w[4][j]*f1o; //we interpolate the intersection point between point_f11 and point_f12 at fraction f1o
     2081
     2082        for (int k=0;k<8;k++){
     2083                for (int i=0;i<NUMVERTICES;i++) {
     2084                        for (int j=0;j<3;j++) xyz[k][j]+=xyz0[i][j]*w[k][i];
     2085                }
     2086        }
     2087
     2088                //point2 can be either i0,i1 or i2. We start the search with i1 and i2 as they have less computational cost in ifs
     2089                if(point2==i2){ /*{{{*/
     2090                        if (f12>1.0-f21){ /*{{{*/
     2091                                if (!istrapeze1 && !istrapeze2){
     2092                                        tria1[0]=5; tria1[1]= 7; tria1[2]= 4;
     2093                                }
     2094                                else if (!istrapeze1 && istrapeze2){
     2095                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 7;
     2096                                        tria2[0]=i0; tria2[1]= 7; tria2[2]= 5;
     2097                                }
     2098                                else if (istrapeze1 && !istrapeze2){
     2099                                        tria1[0]=7; tria1[1]= 6; tria1[2]= 4;
     2100                                        tria2[0]=4; tria2[1]= 6; tria2[2]= i2;
     2101                                }
     2102                                else { //istrapeze1 && istrapeze2
     2103                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 7;
     2104                                        tria2[0]=7; tria2[1]= i1; tria2[2]= 6;
     2105                                }  /*}}}*/
     2106                        }
     2107                        else if (f12<=1.0-f21){ /*{{{*/
     2108                                if (!istrapeze1 && !istrapeze2){
     2109                                }
     2110                                else if (!istrapeze1 && istrapeze2){
     2111                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 4;
     2112                                }
     2113                                else if (istrapeze1 && !istrapeze2){
     2114                                        tria1[0]=5; tria1[1]= 6; tria1[2]= i2;
     2115                                }
     2116                                else { //istrapeze1 && istrapeze2
     2117                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 4;
     2118                                        tria2[0]=4; tria2[1]= i1; tria2[2]= 5;
     2119                                        tria3[0]=5; tria3[1]= i1; tria3[2]= 6;
     2120                                } /*}}}*/
     2121                        }
     2122                }/*}}}*/   
     2123                else if(point2==i1){ /*{{{*/
     2124                        if (f11>1.0-f22){ /*{{{*/
     2125                                if (!istrapeze1 && !istrapeze2){
     2126                                        tria1[0]=6; tria1[1]= 3; tria1[2]= 7;
     2127                                }
     2128                                else if (!istrapeze1 && istrapeze2){
     2129                                        tria1[0]=i0; tria1[1]= 6; tria1[2]= 7;
     2130                                        tria2[0]=i0; tria2[1]= 7; tria2[2]= 4;
     2131                                }
     2132                                else if (istrapeze1 && !istrapeze2){
     2133                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 7;
     2134                                        tria2[0]=7; tria2[1]= i1; tria2[2]= 5;
     2135                                }
     2136                                else { //istrapeze1 && istrapeze2
     2137                                        tria1[0]=7; tria1[1]= 5; tria1[2]= 4;
     2138                                        tria2[0]=4; tria2[1]= 5; tria2[2]= i2;
     2139                                }  /*}}}*/
     2140                        }
     2141                        else if (f11<=1.0-f22){ /*{{{*/
     2142                                if (!istrapeze1 && !istrapeze2){
     2143                                }
     2144                                else if (!istrapeze1 && istrapeze2){
     2145                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 4;
     2146                                }
     2147                                else if (istrapeze1 && !istrapeze2){
     2148                                        tria1[0]=6; tria1[1]= i1; tria1[2]= 5;
     2149                                }
     2150                                else { //istrapeze1 && istrapeze2
     2151                                        tria1[0]=3; tria1[1]= 6; tria1[2]= 4;
     2152                                        tria2[0]=4; tria2[1]= 6; tria2[2]= 5;
     2153                                        tria3[0]=4; tria3[1]= 5; tria3[2]= i2;
     2154                                } /*}}}*/
     2155                        }
     2156                       
     2157                }/*}}}*/
     2158                else{ /*{{{*/
     2159                        if (f11<=f21 && f12>=f22){  /*{{{*/
     2160                                if (!istrapeze1 && !istrapeze2){
     2161                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 7;
     2162                                        tria2[0]=i0; tria2[1]= 7; tria2[2]= 6;
     2163                                }
     2164                                else if (!istrapeze1 && istrapeze2){
     2165                                        tria1[0]=6; tria1[1]= 7; tria1[2]= 4;
     2166                                }
     2167                                else if (istrapeze1 && !istrapeze2){
     2168                                        tria1[0]=3; tria1[1]= 5; tria1[2]= 7;
     2169                                }
     2170                                else { //istrapeze1 && istrapeze2
     2171                                        tria1[0]=7; tria1[1]= 5; tria1[2]= i1;
     2172                                        tria2[0]=7; tria2[1]= i1; tria2[2]= 4;
     2173                                        tria3[0]=4; tria3[1]= i1; tria3[2]= i2;
     2174                                } /*}}}*/
     2175                        }
     2176                        else if (f11>=f21 && f12<=f22){ /*{{{*/
     2177                                if (!istrapeze1 && !istrapeze2){
     2178                                        tria1[0]=i0; tria1[1]= 5; tria1[2]= 7;
     2179                                        tria2[0]=i0; tria2[1]= 4; tria2[2]= 7;
     2180                                }else if (!istrapeze1 && istrapeze2){
     2181                                        tria1[0]=5; tria1[1]= 3; tria1[2]= 7;
     2182                                }else if (istrapeze1 && !istrapeze2){
     2183                                        tria1[0]=4; tria1[1]= 7; tria1[2]= 6;
     2184                                }else { //istrapeze1 && istrapeze2
     2185                                        tria1[0]=3; tria1[1]= i1; tria1[2]= 7;
     2186                                        tria2[0]=7; tria2[1]= i1; tria2[2]= 6;
     2187                                        tria3[0]=6; tria3[1]= i1; tria3[2]= i2;
     2188                                }  /*}}}*/
     2189                        }
     2190                        else if (f11<=f21 && f12<=f22){ /*{{{*/
     2191                                if (!istrapeze1 && !istrapeze2){
     2192                                        tria1[0]=i0; tria1[1]= 3; tria1[2]= 4;
     2193                                }else if (!istrapeze1 && istrapeze2){
     2194                                }else if (istrapeze1 && !istrapeze2){
     2195                                        tria1[0]=3; tria1[1]= 5; tria1[2]= 4;
     2196                                        tria2[0]=4; tria2[1]= 5; tria2[2]= 6;
     2197                                }else { //istrapeze1 && istrapeze2
     2198                                        tria1[0]=6; tria1[1]= 5; tria1[2]= i1;
     2199                                        tria2[0]=6; tria2[1]= i1; tria2[2]= i2;
     2200                                }  /*}}}*/
     2201                        }
     2202                        else if (f11>=f21 && f12>=f22){ /*{{{*/
     2203                                if (!istrapeze1 && !istrapeze2){
     2204                                        tria1[0]=i0; tria1[1]= 5; tria1[2]= 6;
     2205                                }else if (!istrapeze1 && istrapeze2){
     2206                                        tria1[0]=5; tria1[1]= 3; tria1[2]= 6;
     2207                                        tria2[0]=6; tria2[1]= 3; tria2[2]= 4;
     2208                                }else if (istrapeze1 && !istrapeze2){
     2209                                }else { //istrapeze1 && istrapeze2
     2210                                        tria1[0]=4; tria1[1]= 3; tria1[2]= i1;
     2211                                        tria2[0]=i1; tria2[1]= i2; tria2[2]= 4;
     2212                                }  /*}}}*/
     2213                        }
     2214                } /*}}}*/
     2215
     2216        if(tria1[0]>-1){
     2217                for (int i=0;i<NUMVERTICES;i++){
     2218                        for (int j=0;j<3;j++) {
     2219                                xyz1[i][j]=xyz[tria1[i]][j];
     2220                                w1[i][j]=w[tria1[i]][j];
     2221                        }
     2222                }
     2223                area1= GetTriangleAreaSpherical(xyz1);
     2224        }
     2225        if(tria2[0]>-1){
     2226                for (int i=0;i<NUMVERTICES;i++){
     2227                        for (int j=0;j<3;j++) {
     2228                                xyz2[i][j]=xyz[tria2[i]][j];
     2229                                w2[i][j]=w[tria2[i]][j];
     2230                        }
     2231                }
     2232                area2= GetTriangleAreaSpherical(xyz2);
     2233        }
     2234        if(tria3[0]>-1){
     2235                for (int i=0;i<NUMVERTICES;i++){
     2236                        for (int j=0;j<3;j++) {
     2237                                xyz3[i][j]=xyz[tria3[i]][j];
     2238                                w3[i][j]=w[tria3[i]][j];
     2239                        }
     2240                }
     2241                area3= GetTriangleAreaSpherical(xyz3);
     2242        }
     2243
     2244        areasub=area1+area2+area3;
     2245
     2246        if (areasub>0){
     2247                for (int j=0;j<3;j++){
     2248                        for (int i=0;i<NUMVERTICES;i++) {
     2249                                loadweights[j]=w1[i][j]*area1 + w2[i][j]*area2 + w3[i][j]*area3;
     2250                                barycenter[j]+=xyz1[i][j]*area1+xyz2[i][j]*area2+xyz3[i][j]*area3;
     2251                        }
     2252                        loadweights[j]/=area;
     2253                        barycenter[j]/=areasub *3.0;
     2254                }
     2255                *platbar=asin(barycenter[2]/sqrt(pow(barycenter[0],2.0)+pow(barycenter[1],2.0)+pow(barycenter[2],2.0)))*180.0/M_PI;
     2256                *plongbar=atan2(barycenter[1],barycenter[0])*180.0/M_PI;
     2257        }
     2258        else {
     2259                for(int j=0;j<3;j++)loadweights[j]=0.0;
     2260                *platbar=late;
     2261                *plongbar=longe;
     2262        }
     2263        *ploadarea=areasub;
     2264
     2265} /*}}}*/
    18422266IssmDouble Tria::GetIcefrontArea(){/*{{{*/
    18432267
    … …  
    55555979#endif
    55565980#ifdef _HAVE_SEALEVELCHANGE_
    5557 //old code
    5558 void       Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/
    5559 
    5560         masks->isiceonly[this->lid]=this->IsIceOnlyInElement();
    5561         masks->isoceanin[this->lid]=this->IsOceanInElement();
    5562 
    5563         /*are we fully floating:*/
    5564         Input* gr_input=this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
    5565         if (gr_input->GetInputMax()<=0)masks->isfullyfloating[this->lid]=true;
    5566         else masks->isfullyfloating[this->lid]=false;
    5567 
    5568         /*are we not fully grounded: */
    5569         if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true;
    5570         else masks->notfullygrounded[this->lid]=false;
    5571 }
    5572 /*}}}*/
    55735981void       Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/
    55745982
    … …  
    56766084}
    56776085/*}}}*/
    5678 void       Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/
    5679 
    5680         /*diverse:*/
     6086void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/
     6087
     6088        /*Declarations:{{{*/
    56816089        int nel;
    56826090        IssmDouble area,planetarea,planetradius;
    … …  
    56846092        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
    56856093        IssmDouble rho_earth;
     6094        IssmDouble NewtonG;
     6095        IssmDouble g;
    56866096        IssmDouble lati,longi;
     6097        IssmDouble latitude[NUMVERTICES];
     6098        IssmDouble longitude[NUMVERTICES];
    56876099        IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
    5688         int sidlist[NUMVERTICES];
    5689         int sid;
     6100        IssmDouble xyz_list[NUMVERTICES][3];
    56906101
    56916102        #ifdef _HAVE_RESTRICT_
    … …  
    56986109        IssmDouble* __restrict__ U_elastic_precomputed=NULL;
    56996110        IssmDouble* __restrict__ H_elastic_precomputed=NULL;
    5700         IssmDouble* __restrict__ lates=NULL;
    5701         IssmDouble* __restrict__ longes=NULL;
    57026111        #else
    57036112        IssmDouble* G=NULL;
    … …  
    57096118        IssmDouble* U_elastic_precomputed=NULL;
    57106119        IssmDouble* H_elastic_precomputed=NULL;
    5711         IssmDouble* lates=NULL;
    5712         IssmDouble* longes=NULL;
    57136120        #endif
    57146121       
    … …  
    57286135        IssmDouble* load_love_k  = NULL;
    57296136        IssmDouble  tide_love_k2secular;
    5730         IssmDouble  moi_e, moi_p, omega, g;
     6137        IssmDouble  moi_e, moi_p, omega;
    57316138        IssmDouble Grotm1[3];
    57326139        IssmDouble Grotm2[3];
    57336140        IssmDouble Grotm3[3];
    57346141        IssmDouble pre;
    5735 
    5736 
    5737         /*recover parameters: */
     6142        /*}}}*/
     6143        /*Recover parameters:{{{ */
    57386144        rho_earth=FindParam(MaterialsEarthDensityEnum);
    57396145        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    … …  
    57436149        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    57446150        this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
     6151        this->parameters->FindParam(&NewtonG,ConstantsNewtonGravityEnum);
    57456152        this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    57466153
    … …  
    57546161                parameters->FindParam(&omega,RotationalAngularVelocityEnum);
    57556162        }
     6163        /*}}}*/
    57566164
    57576165        /*early return:*/
    57586166        if(!computerigid)return;
    57596167
    5760         /*recover precomputed green function kernels:*/
     6168        /*Recover precomputed green function kernels:{{{*/
    57616169        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter);
    57626170        parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
    5763 
    5764         /*allocate arrays:*/
    5765         G=xNewZeroInit<IssmDouble>(3*nel);
    5766         longes=xNewZeroInit<IssmDouble>(nel);
    5767         lates=xNewZeroInit<IssmDouble>(nel);
    57686171
    57696172        if(computeelastic){
    … …  
    57776180                parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
    57786181
     6182               
     6183        }
     6184        /*}}}*/
     6185        /*Compute lat long of all vertices in the element:{{{*/
     6186        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6187        for(int i=0;i<NUMVERTICES;i++){
     6188                latitude[i]= asin(xyz_list[i][2]/planetradius);
     6189                if((xyz_list[i][2]/planetradius)==1.0)latitude[i]=M_PI/2;
     6190                longitude[i]= atan2(xyz_list[i][1],xyz_list[i][0]);
     6191        }
     6192        /*}}}*/
     6193        /*Compute green functions:{{{ */
     6194        constant=3/rho_earth/planetarea;
     6195       
     6196        //Allocate:
     6197        G=xNewZeroInit<IssmDouble>(3*nel);
     6198        if(computeelastic){
    57796199                GU=xNewZeroInit<IssmDouble>(3*nel);
    57806200                if(horiz){
    … …  
    57836203                }
    57846204        }
    5785 
    5786         /*compute centroids of all elements:*/
    5787         for(int i=0;i<nel;i++){
    5788                 lates[i]= asin(zze[i]/planetradius);
    5789                 longes[i]= atan2(yye[i],xxe[i]);
    5790         }
    5791        
    5792        
    5793         constant=3/rho_earth/planetarea;
    5794 
    5795         /*Recover vertex absolute id: */
    5796         this->GetVerticesSidList(&sidlist[0]);
    57976205
    57986206        for(int e=0;e<nel;e++){
    … …  
    58026210                        IssmDouble delPhi,delLambda;
    58036211                        /*recover info for this element and vertex:*/
    5804                         IssmDouble late=lates[e];
    5805                         IssmDouble longe=longes[e];
    5806                         sid=sidlist[i];
     6212                        IssmDouble late= asin(zze[e]/sqrt( pow(xxe[e],2.0)+ pow(yye[e],2.0)+ pow(zze[e],2.0)));
     6213                        IssmDouble longe= atan2(yye[e],xxe[e]);
     6214                       
     6215                        lati=latitude[i];
     6216                        longi=longitude[i];
    58076217
    58086218                        /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    5809                         lati=latitude[sid]/180.*M_PI; longi=longitude[sid]/180.*M_PI;
    58106219                        delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    58116220                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    58126221                        index=reCast<int,IssmDouble>( alpha/M_PI*reCast<IssmDouble,int>(M-1) );
     6222                        _assert_(index>0 && index<M);
    58136223
    58146224                        /*Rigid earth gravitational perturbation: */
    … …  
    58256235                                if(horiz){
    58266236                                        /*Compute azimuths, both north and east components: */
    5827                                         x = xx[sid]; y = yy[sid]; z = zz[sid];
    5828                                         if(latitude[sid]==90){
     6237                                        x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
     6238                                        if(lati==M_PI/2){
    58296239                                                x=1e-12; y=1e-12;
    58306240                                        }
    5831                                         if(latitude[sid]==-90){
     6241                                        if(lati==-M_PI/2){
    58326242                                                x=1e-12; y=1e-12;
    58336243                                        }
    … …  
    58416251                        }
    58426252                }
    5843         }
    5844         if(computerotation){
    5845 
    5846                 for(int i=0;i<3;i++){
    5847                         sid=sidlist[i];
    5848                         lati=latitude[sid]/180.*M_PI;
    5849                         longi=longitude[sid]/180.*M_PI;
    5850 
    5851                         pre=((1.0+tide_love_k[2]-tide_love_h[2])/9.81)*pow(omega*planetradius,2.0);
    5852                         Grotm1[i]= - pre* 0.5*sin(2.*lati)*cos(longi);
    5853                         Grotm2[i]= - pre* 0.5*sin(2.*lati)*sin(longi);
    5854                         Grotm3[i]= - pre* (1.0/6.0 - 0.5*cos(2.0*lati));
    5855                 }
    5856                 this->AddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum);
    5857                 this->AddInput(SealevelGrotm2Enum,&Grotm2[0],P1Enum);
    5858                 this->AddInput(SealevelGrotm3Enum,&Grotm3[0],P1Enum);
    58596253        }
    58606254
    … …  
    58696263                }
    58706264        }
    5871 
    5872         /*Free allocations:*/
     6265        /*}}}*/
     6266        /*Compute rotation kernel:{{{*/
     6267        if(computerotation){
     6268
     6269                /*What is the gravity at planet's surface: */
     6270                g=4.0/3.0*M_PI*rho_earth*NewtonG*planetradius;
     6271
     6272                for(int i=0;i<3;i++){
     6273                        lati=latitude[i];
     6274                        longi=longitude[i];
     6275
     6276                        pre=((1.0+tide_love_k[2]-tide_love_h[2])/g)*pow(omega*planetradius,2.0);
     6277                        Grotm1[i]= - pre* 0.5*sin(2.*lati)*cos(longi);
     6278                        Grotm2[i]= - pre* 0.5*sin(2.*lati)*sin(longi);
     6279                        Grotm3[i]= - pre* (1.0/6.0 - 0.5*cos(2.0*lati));
     6280                }
     6281                this->AddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum);
     6282                this->AddInput(SealevelGrotm2Enum,&Grotm2[0],P1Enum);
     6283                this->AddInput(SealevelGrotm3Enum,&Grotm3[0],P1Enum);
     6284        }
     6285        /*}}}*/
     6286        /*Free allocations:{{{*/
    58736287        #ifdef _HAVE_RESTRICT_
    58746288        delete G;
    … …  
    58806294                }
    58816295        }
    5882         delete lates;
    5883         delete longes;
    58846296        #else
    58856297        xDelete(G);
    … …  
    58916303                }
    58926304        }
    5893         delete lates;
    5894         delete longes;
    58956305        #endif
    5896 
     6306        /*}}}*/
     6307
     6308       
    58976309        return;
    5898 }
    5899 /*}}}*/
    5900 void       Tria::SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){ /*{{{*/
    5901 
    5902         //Compute sea level barystatic loads (ice, water column or bottom pressure, see Farrel and Clarke, Equ. 4)
    5903 
    5904         /*diverse:*/
     6310
     6311}
     6312/*}}}*/
     6313void       Tria::SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/
     6314
     6315        /* Classic buildup of load weights, centroids and areas *for elements which are fully inside a mask.
     6316         * At the same time, we'll tag the elements that are fractionally only inside a mask*/
     6317       
     6318        IssmDouble loadweights[3]={0};
    59056319        IssmDouble area;
    5906         IssmDouble phi_ice=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
    5907         IssmDouble phi_water=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
    5908         IssmDouble phi_bp=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic
    5909         IssmDouble I=0; //Do not change this!
    5910         IssmDouble W=0; //Do not change this!
    5911         IssmDouble BP=0; //Do not change this!
    5912 
    5913                 bool notfullygrounded=false;
    5914         bool computerigid= false;
    5915         int  glfraction=1;
    5916         int  npartice,nparthydro,npartocean;
     6320        IssmDouble loadweightsocean[3]; //to keep memory of these loads, no need to recompute for bottom pressure.
     6321        IssmDouble xyz_list[NUMVERTICES][3];
     6322        IssmDouble planetradius;
     6323        IssmDouble late,longe;
    59176324
    59186325        /*flags:*/
    59196326        bool isocean=false;
    5920         bool ishydro=false;
     6327        bool isoceanonly=false;
    59216328        bool isice=false;
     6329        bool isiceonly=false;
     6330        bool  computeice=false;
     6331        bool  computebp=false;
     6332        bool  computehydro=false;
     6333
     6334        /*constants:*/
     6335        IssmDouble constant=0;
     6336       
     6337        /*recover parameters:*/
     6338        this->parameters->FindParam(&computeice,TransientIsmasstransportEnum);
     6339        this->parameters->FindParam(&computebp,TransientIsoceantransportEnum);
     6340        this->parameters->FindParam(&computehydro,TransientIshydrologyEnum);
     6341        this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
     6342
     6343        /*get vertex information:*/
     6344        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6345
     6346        /*answer mask questions:*/
     6347        isiceonly=this->IsIceOnlyInElement();
     6348        isice=this->IsIceInElement();
     6349        isoceanonly=this->IsOceanOnlyInElement();
     6350        isocean=this->IsOceanInElement();
     6351        slgeom->isoceanin[this->lid]=isocean; //keep track for later.
     6352        area=areae[this->sid];
     6353
     6354        /*set barycentre for all elements, to be updated for fractional loads in the next routine: */
     6355        //late= asin(zze[this->sid]/planetradius)*180.0/M_PI;
     6356        late= asin(zze[this->sid]/sqrt( pow(xxe[this->sid],2.0)+ pow(yye[this->sid],2.0)+ pow(zze[this->sid],2.0)))*180.0/M_PI;
     6357        longe= atan2(yye[this->sid],xxe[this->sid])*180.0/M_PI;
     6358        slgeom->longe[this->lid]=longe;
     6359        slgeom->late[this->lid]=late;
     6360
     6361        /*compute fractional areas and load weights for ocean:*/
     6362        if(isoceanonly){
     6363                slgeom->LoadArea[SLGEOM_OCEAN][this->lid]=area;
     6364                for(int i=0;i<NUMVERTICES;i++) slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]=1.0/3.0;
     6365
     6366                #ifdef _ISSM_DEBUG_ /*{{{*/
     6367                /*Inform mask: */
     6368                constant=1.0;
     6369                for(int i=0;i<NUMVERTICES;i++) loadweightsocean[i]=1.0/3.0;
     6370                this->AddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum);
     6371                this->AddInput(SealevelBarystaticOceanWeightsEnum,loadweightsocean,P1DGEnum);
     6372                this->AddInput(SealevelBarystaticOceanAreaEnum,&area,P0Enum);
     6373                #endif /*}}}*/
     6374        }
     6375        else if(!isocean){
     6376                slgeom->LoadArea[SLGEOM_OCEAN][this->lid]=0;
     6377                for(int i=0;i<NUMVERTICES;i++) slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]=0.0;
     6378                #ifdef _ISSM_DEBUG_ /*{{{*/
     6379                /*Inform mask: */
     6380                constant=0.0;
     6381                for(int i=0;i<NUMVERTICES;i++) loadweightsocean[i]=0.0;
     6382                this->AddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum);
     6383                this->AddInput(SealevelBarystaticOceanWeightsEnum,loadweightsocean,P1DGEnum);
     6384                this->AddInput(SealevelBarystaticOceanAreaEnum,&constant,P0Enum);
     6385                #endif /*}}}*/
     6386        }
     6387        else{
     6388                slgeom->issubelement[SLGEOM_OCEAN][this->lid]=true;
     6389                slgeom->nsubel[SLGEOM_OCEAN]++;
     6390        }
     6391               
     6392        /*early return if we are not on an ice sheet , and we are not requesting
     6393         *hydrology or bottom pressure loads :*/
     6394        if(!computebp && !computehydro){
     6395                if(!isice  || isoceanonly) {
     6396                        #ifdef _ISSM_DEBUG_
     6397                        constant=0;
     6398                        this->AddInput(SealevelBarystaticIceMaskEnum,&constant,P0Enum);
     6399                        this->AddInput(SealevelBarystaticIceAreaEnum,&constant,P0Enum);
     6400                        this->AddInput(SealevelBarystaticIceWeightsEnum,loadweights,P1DGEnum);
     6401                        this->AddInput(SealevelBarystaticHydroMaskEnum,&constant,P0Enum);
     6402                        this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum);
     6403                        this->AddInput(SealevelBarystaticHydroAreaEnum,&constant,P0Enum);
     6404                        this->AddInput(SealevelBarystaticBpMaskEnum,&constant,P0Enum);
     6405                        this->AddInput(SealevelBarystaticBpWeightsEnum,loadweights,P1DGEnum);
     6406                        this->AddInput(SealevelBarystaticBpAreaEnum,&constant,P0Enum);
     6407                        #endif
     6408                        for(int i=0;i<NUMVERTICES;i++){
     6409                                slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid]=0;
     6410                                slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]=0;
     6411                        }
     6412                        slgeom->LoadArea[SLGEOM_ICE][this->lid]=0;
     6413                        slgeom->LoadArea[SLGEOM_WATER][this->lid]=0;
     6414                        return;
     6415                }
     6416        }
     6417
     6418        /*early return if we are fully floating and we are not doing bottom pressure loads:*/
     6419        if(!computebp){
     6420                if (isoceanonly){
     6421                        #ifdef _ISSM_DEBUG_
     6422                        constant=0;
     6423                        this->AddInput(SealevelBarystaticIceMaskEnum,&constant,P0Enum);
     6424                        this->AddInput(SealevelBarystaticIceWeightsEnum,loadweights,P1DGEnum);
     6425                        this->AddInput(SealevelBarystaticIceAreaEnum,&constant,P0Enum);
     6426                        this->AddInput(SealevelBarystaticHydroMaskEnum,&constant,P0Enum);
     6427                        this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum);
     6428                        this->AddInput(SealevelBarystaticHydroAreaEnum,&constant,P0Enum);
     6429                        this->AddInput(SealevelBarystaticBpMaskEnum,&constant,P0Enum);
     6430                        this->AddInput(SealevelBarystaticBpWeightsEnum,loadweights,P1DGEnum);
     6431                        this->AddInput(SealevelBarystaticBpAreaEnum,&constant,P0Enum);
     6432                        #endif
     6433                        for(int i=0;i<NUMVERTICES;i++){
     6434                                slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid]=0;
     6435                                slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]=0;
     6436                        }
     6437                        slgeom->LoadArea[SLGEOM_ICE][this->lid]=0;
     6438                        slgeom->LoadArea[SLGEOM_WATER][this->lid]=0;
     6439                        return;
     6440                }
     6441        }
     6442
     6443        /*early return if we are not on the ocean and we are not doing ice mass transport of
     6444         * hydrology:*/
     6445        if(!computeice  && !computehydro){
     6446                if(!isocean){
     6447                        #ifdef _ISSM_DEBUG_
     6448                        constant=0;
     6449                        this->AddInput(SealevelBarystaticIceMaskEnum,&constant,P0Enum);
     6450                        this->AddInput(SealevelBarystaticIceWeightsEnum,loadweights,P1DGEnum);
     6451                        this->AddInput(SealevelBarystaticIceAreaEnum,&constant,P0Enum);
     6452                        this->AddInput(SealevelBarystaticHydroMaskEnum,&constant,P0Enum);
     6453                        this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum);
     6454                        this->AddInput(SealevelBarystaticHydroAreaEnum,&constant,P0Enum);
     6455                        this->AddInput(SealevelBarystaticBpMaskEnum,&constant,P0Enum);
     6456                        this->AddInput(SealevelBarystaticBpWeightsEnum,loadweights,P1DGEnum);
     6457                        this->AddInput(SealevelBarystaticBpAreaEnum,&constant,P0Enum);
     6458                        #endif
     6459                        for(int i=0;i<NUMVERTICES;i++){
     6460                                slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid]=0;
     6461                                slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]=0;
     6462                        }
     6463                        slgeom->LoadArea[SLGEOM_ICE][this->lid]=0;
     6464                        slgeom->LoadArea[SLGEOM_WATER][this->lid]=0;
     6465                        return;
     6466                }
     6467        }
     6468       
     6469        /*Deal with ice loads if we are on grounded ice:*/
     6470        if(isice && !isoceanonly && computeice){
     6471                if(isiceonly){
     6472                        slgeom->LoadArea[SLGEOM_ICE][this->lid]=area;
     6473                        for(int i=0;i<NUMVERTICES;i++) slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid]=1.0/3.0;
     6474
     6475                        #ifdef _ISSM_DEBUG_ /*{{{*/
     6476                        /*Inform mask: */
     6477                        constant=1.0;
     6478                        for(int i=0;i<NUMVERTICES;i++) loadweights[i]=1.0/3.0;
     6479                        this->AddInput(SealevelBarystaticIceMaskEnum,&constant,P0Enum);
     6480                        this->AddInput(SealevelBarystaticIceWeightsEnum,loadweights,P1DGEnum);
     6481                        this->AddInput(SealevelBarystaticIceAreaEnum,&area,P0Enum);
     6482                        #endif /*}}}*/
     6483                }
     6484                else{
     6485                        slgeom->issubelement[SLGEOM_ICE][this->lid]=true;
     6486                        slgeom->nsubel[SLGEOM_ICE]++;
     6487                }
     6488        }
     6489
     6490        /*Deal with water loads if we are on ground:*/
     6491        if(!isoceanonly && computehydro){
     6492
     6493                if(!isocean){
     6494                        slgeom->LoadArea[SLGEOM_WATER][this->lid]=area;
     6495                        for(int i=0;i<NUMVERTICES;i++) slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]=1.0/3.0;
     6496
     6497                        #ifdef _ISSM_DEBUG_ /*{{{*/
     6498                        /*Inform mask: */
     6499                        constant=1.0;
     6500                        for(int i=0;i<NUMVERTICES;i++) loadweights[i]=1.0/3.0;
     6501                        this->AddInput(SealevelBarystaticHydroMaskEnum,&constant,P0Enum);
     6502                        this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum);
     6503                        this->AddInput(SealevelBarystaticHydroAreaEnum,&area,P0Enum);
     6504                        #endif /*}}}*/
     6505                }
     6506                else{
     6507                        slgeom->issubelement[SLGEOM_WATER][this->lid]=true;
     6508                        slgeom->nsubel[SLGEOM_WATER]++;
     6509                }
     6510        }
     6511       
     6512}
     6513/*}}}*/
     6514void       Tria::SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){ /*{{{*/
     6515
     6516        /* Classic buildup of load weights, centroids and areas *for elements which are fully inside a mask.
     6517         * At the same time, we'll tag the elements that are fractionally only inside a mask*/
     6518       
     6519        IssmDouble loadweights[3]={0};
     6520        IssmDouble area,loadarea;
     6521        IssmDouble loadareaocean;
     6522        IssmDouble loadweightsocean[3]; //to keep memory of these loads, no need to recompute for bottom pressure.
     6523        IssmDouble xyz_list[NUMVERTICES][3];
     6524        IssmDouble latbar=0;
     6525        IssmDouble longbar=0;
     6526        IssmDouble constant;
     6527        IssmDouble nanconstant=NAN;
     6528       
     6529        /*get vertex and area information:*/
     6530        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6531        area=areae[this->sid];
     6532
     6533        #ifdef _ISSM_DEBUG_
     6534        this->AddInput(SealevelBarystaticIceLatbarEnum,&nanconstant,P0Enum);
     6535        this->AddInput(SealevelBarystaticIceLongbarEnum,&nanconstant,P0Enum);
     6536        this->AddInput(SealevelBarystaticHydroLatbarEnum,&nanconstant,P0Enum);
     6537        this->AddInput(SealevelBarystaticHydroLongbarEnum,&nanconstant,P0Enum);
     6538        this->AddInput(SealevelBarystaticOceanLatbarEnum,&nanconstant,P0Enum);
     6539        this->AddInput(SealevelBarystaticOceanLongbarEnum,&nanconstant,P0Enum);
     6540        #endif
     6541       
     6542        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
     6543                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
     6544               
     6545                this->GetNodalWeightsAndAreaAndCentroidsFromLeveset(&loadweightsocean[0],&loadareaocean,&latbar, &longbar, slgeom->late[this->lid], slgeom->longe[this->lid], area, MaskOceanLevelsetEnum);
     6546                slgeom->LoadArea[SLGEOM_OCEAN][this->lid]=loadareaocean;
     6547                slgeom->vareae_subel[SLGEOM_OCEAN]->SetValue(intj,loadareaocean,INS_VAL);
     6548                slgeom->vlatbarycentre[SLGEOM_OCEAN]->SetValue(intj,latbar,INS_VAL);
     6549                slgeom->vlongbarycentre[SLGEOM_OCEAN]->SetValue(intj,longbar,INS_VAL);
     6550
     6551                for(int i=0;i<NUMVERTICES;i++) slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]=loadweightsocean[i];
     6552                #ifdef _ISSM_DEBUG_ /*{{{*/
     6553                /*Inform mask: */
     6554                constant=loadareaocean/area;
     6555                this->AddInput(SealevelBarystaticOceanMaskEnum,&constant,P0Enum);
     6556                this->AddInput(SealevelBarystaticOceanWeightsEnum,loadweightsocean,P1DGEnum);
     6557                this->AddInput(SealevelBarystaticOceanAreaEnum,&loadareaocean,P0Enum);
     6558
     6559                this->AddInput(SealevelBarystaticOceanLatbarEnum,&latbar,P0Enum);
     6560                this->AddInput(SealevelBarystaticOceanLongbarEnum,&longbar,P0Enum);
     6561                #endif /*}}}*/
     6562        }
     6563        if(slgeom->issubelement[SLGEOM_ICE][this->lid]){
     6564                int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid];
     6565
     6566                this->GetNodalWeightsAndAreaAndCentroidsFromLeveset(&loadweights[0],&loadarea,&latbar, &longbar, slgeom->late[this->lid], slgeom->longe[this->lid], area, -MaskOceanLevelsetEnum,MaskIceLevelsetEnum);
     6567
     6568                slgeom->LoadArea[SLGEOM_ICE][this->lid]=loadarea;
     6569                slgeom->vareae_subel[SLGEOM_ICE]->SetValue(intj,loadarea,INS_VAL);
     6570                slgeom->vlatbarycentre[SLGEOM_ICE]->SetValue(intj,latbar,INS_VAL);
     6571                slgeom->vlongbarycentre[SLGEOM_ICE]->SetValue(intj,longbar,INS_VAL);
     6572
     6573                for(int i=0;i<NUMVERTICES;i++)slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid]=loadweights[i];
     6574
     6575                #ifdef _ISSM_DEBUG_
     6576                /*Inform mask: */
     6577                constant=loadarea/area;
     6578                this->AddInput(SealevelBarystaticIceMaskEnum,&constant,P0Enum);
     6579                this->AddInput(SealevelBarystaticIceWeightsEnum,loadweights,P1DGEnum);
     6580                this->AddInput(SealevelBarystaticIceAreaEnum,&loadarea,P0Enum);
     6581
     6582                this->AddInput(SealevelBarystaticIceLatbarEnum,&latbar,P0Enum);
     6583                this->AddInput(SealevelBarystaticIceLongbarEnum,&longbar,P0Enum);
     6584
     6585                #endif
     6586        }
     6587        if(slgeom->issubelement[SLGEOM_WATER][this->lid]){
     6588                int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid];
     6589
     6590                this->GetNodalWeightsAndAreaAndCentroidsFromLeveset(&loadweights[0],&loadarea,&latbar, &longbar, slgeom->late[this->lid], slgeom->longe[this->lid], area, -MaskOceanLevelsetEnum);
     6591
     6592                slgeom->LoadArea[SLGEOM_WATER][this->lid]=loadarea;
     6593                slgeom->vareae_subel[SLGEOM_WATER]->SetValue(intj,loadarea,INS_VAL);
     6594                slgeom->vlatbarycentre[SLGEOM_WATER]->SetValue(intj,latbar,INS_VAL);
     6595                slgeom->vlongbarycentre[SLGEOM_WATER]->SetValue(intj,longbar,INS_VAL);
     6596
     6597                for(int i=0;i<NUMVERTICES;i++)slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]=loadweights[i];
     6598
     6599                #ifdef _ISSM_DEBUG_
     6600                /*Inform mask: */
     6601                constant=loadarea/area;
     6602                this->AddInput(SealevelBarystaticHydroMaskEnum,&constant,P0Enum);
     6603                this->AddInput(SealevelBarystaticHydroWeightsEnum,loadweights,P1DGEnum);
     6604                this->AddInput(SealevelBarystaticHydroAreaEnum,&loadarea,P0Enum);
     6605               
     6606                this->AddInput(SealevelBarystaticHydroLatbarEnum,&latbar,P0Enum);
     6607                this->AddInput(SealevelBarystaticHydroLongbarEnum,&longbar,P0Enum);
     6608
     6609                #endif
     6610        }
     6611       
     6612}
     6613/*}}}*/
     6614void       Tria::SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom){ /*{{{*/
     6615
     6616        /*Declarations:{{{*/
     6617        int nel;
     6618        IssmDouble planetarea,planetradius;
     6619        IssmDouble constant,ratioe;
     6620        IssmDouble rho_earth;
     6621        IssmDouble lati,longi;
     6622        IssmDouble latitude[NUMVERTICES];
     6623        IssmDouble longitude[NUMVERTICES];
     6624        IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
     6625        IssmDouble xyz_list[NUMVERTICES][3];
     6626
     6627        #ifdef _HAVE_RESTRICT_
     6628        IssmDouble* __restrict__ G_elastic_precomputed=NULL;
     6629        IssmDouble* __restrict__ G_rigid_precomputed=NULL;
     6630        IssmDouble* __restrict__ U_elastic_precomputed=NULL;
     6631        IssmDouble* __restrict__ H_elastic_precomputed=NULL;
     6632        IssmDouble** __restrict__ Gsubel=NULL;
     6633        IssmDouble** __restrict__ GUsubel=NULL;
     6634        IssmDouble** __restrict__ GNsubel=NULL;
     6635        IssmDouble** __restrict__ GEsubel=NULL;
     6636
     6637        #else
     6638        IssmDouble* G_elastic_precomputed=NULL;
     6639        IssmDouble* G_rigid_precomputed=NULL;
     6640        IssmDouble* U_elastic_precomputed=NULL;
     6641        IssmDouble* H_elastic_precomputed=NULL;
     6642        IssmDouble** Gsubel=NULL;
     6643        IssmDouble** GUsubel=NULL;
     6644        IssmDouble** GNsubel=NULL;
     6645        IssmDouble** GEsubel=NULL;
     6646        #endif
     6647       
     6648        /*elastic green function:*/
     6649        int index;
     6650        int         M;
     6651
     6652        /*Computational flags:*/
     6653        bool computerigid = false;
     6654        bool computeelastic = false;
     6655        int  horiz;
     6656
     6657        /*}}}*/
     6658        /*Recover parameters:{{{ */
     6659        rho_earth=FindParam(MaterialsEarthDensityEnum);
     6660        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
     6661        this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
     6662        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
     6663        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
     6664        this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
     6665        this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
     6666        /*}}}*/
     6667
     6668        /*early return:*/
     6669        if(!computerigid)return;
     6670
     6671        /*Recover precomputed green function kernels:{{{*/
     6672        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGRigidEnum)); _assert_(parameter);
     6673        parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
     6674
     6675        if(computeelastic){
     6676                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGElasticEnum)); _assert_(parameter);
     6677                parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
     6678
     6679                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHElasticEnum)); _assert_(parameter);
     6680                parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
     6681
     6682                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUElasticEnum)); _assert_(parameter);
     6683                parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
     6684        }
     6685        /*}}}*/
     6686        /*Compute lat long of all vertices in the element:{{{*/
     6687        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6688        for(int i=0;i<NUMVERTICES;i++){
     6689                latitude[i]= asin(xyz_list[i][2]/planetradius);
     6690                longitude[i]= atan2(xyz_list[i][1],xyz_list[i][0]);
     6691        }
     6692        /*}}}*/
     6693        /*Compute green functions:{{{ */
     6694        constant=3/rho_earth/planetarea;
     6695
     6696        Gsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
     6697        if(computeelastic){
     6698                GUsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
     6699                if(horiz){
     6700                        GNsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
     6701                        GEsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
     6702                }
     6703        }
     6704       
     6705        //Allocate:
     6706        for(int l=0;l<SLGEOM_NUMLOADS;l++){
     6707                int nbar=slgeom->nbar[l];
     6708                Gsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
     6709                if(computeelastic){
     6710                        GUsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
     6711                        if(horiz){
     6712                                GNsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
     6713                                GEsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
     6714                        }
     6715                }
     6716
     6717                for(int e=0;e<nbar;e++){
     6718                        ratioe=constant*slgeom->area_subel[l][e];
     6719                        for (int i=0;i<3;i++){
     6720                                IssmDouble alpha;
     6721                                IssmDouble delPhi,delLambda;
     6722                                /*recover info for this element and vertex:*/
     6723                                IssmDouble late= slgeom->latbarycentre[l][e];
     6724                                IssmDouble longe= slgeom->longbarycentre[l][e];
     6725                                late=late/180*M_PI;
     6726                                longe=longe/180*M_PI;
     6727
     6728                                lati=latitude[i];
     6729                                longi=longitude[i];
     6730
     6731                                /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
     6732                                delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
     6733                                alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0)));
     6734                                index=reCast<int,IssmDouble>( alpha/M_PI*reCast<IssmDouble,int>(M-1) );
     6735
     6736                                /*Rigid earth gravitational perturbation: */
     6737                                Gsubel[l][i*nbar+e]+=G_rigid_precomputed[index];
     6738
     6739                                if(computeelastic){
     6740                                        Gsubel[l][i*nbar+e]+=G_elastic_precomputed[index];
     6741                                }
     6742                                Gsubel[l][i*nbar+e]*=ratioe;
     6743
     6744                                /*Elastic components:*/
     6745                                if(computeelastic){
     6746                                        GUsubel[l][i*nbar+e] =  ratioe * U_elastic_precomputed[index];
     6747                                        if(horiz){
     6748                                                /*Compute azimuths, both north and east components: */
     6749                                                x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
     6750                                                if(lati==90){
     6751                                                        x=1e-12; y=1e-12;
     6752                                                }
     6753                                                if(lati==-90){
     6754                                                        x=1e-12; y=1e-12;
     6755                                                }
     6756                                                IssmDouble xbar=planetradius*cos(late)*cos(longe);
     6757                                                IssmDouble ybar=planetradius*cos(late)*sin(longe);
     6758                                                IssmDouble zbar=planetradius*sin(late);
     6759
     6760                                                dx = xbar-x; dy = ybar-y; dz = zbar-z;
     6761                                                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);
     6762                                                E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
     6763
     6764                                                GNsubel[l][i*nbar+e] = ratioe*H_elastic_precomputed[index]*N_azim;
     6765                                                GEsubel[l][i*nbar+e] = ratioe*H_elastic_precomputed[index]*E_azim;
     6766                                        }
     6767                                }
     6768                        }
     6769                }
     6770        }
     6771
     6772        /*Save all these arrayins for each element:*/
     6773        for (int l=0;l<SLGEOM_NUMLOADS;l++){
     6774                this->inputs->SetArrayInput(slgeom->GEnum(l),this->lid,Gsubel[l],slgeom->nbar[l]*3);
     6775                if(computeelastic){
     6776                        this->inputs->SetArrayInput(slgeom->GUEnum(l),this->lid,GUsubel[l],slgeom->nbar[l]*3);
     6777                        if(horiz){
     6778                                this->inputs->SetArrayInput(slgeom->GNEnum(l),this->lid,GNsubel[l],slgeom->nbar[l]*3);
     6779                                this->inputs->SetArrayInput(slgeom->GEEnum(l),this->lid,GEsubel[l],slgeom->nbar[l]*3);
     6780                        }
     6781                }
     6782        }
     6783        /*}}}*/
     6784        /*Free memory:{{{*/
     6785        for (int l=0;l<SLGEOM_NUMLOADS;l++){
     6786                xDelete<IssmDouble>(Gsubel[l]);
     6787                if(computeelastic){
     6788                        xDelete<IssmDouble>(GUsubel[l]);
     6789                        if(horiz){
     6790                                xDelete<IssmDouble>(GNsubel[l]);
     6791                                xDelete<IssmDouble>(GEsubel[l]);
     6792                        }
     6793                }
     6794        }
     6795        xDelete<IssmDouble*>(Gsubel);
     6796        if(computeelastic){
     6797                xDelete<IssmDouble*>(GUsubel);
     6798                if(horiz){
     6799                        xDelete<IssmDouble*>(GNsubel);
     6800                        xDelete<IssmDouble*>(GEsubel);
     6801                }
     6802        }
     6803        /*}}}*/
     6804        return;
     6805
     6806}
     6807/*}}}*/
     6808void       Tria::SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/
     6809
     6810        /*Inputs:*/
     6811        IssmDouble I[NUMVERTICES];
     6812        IssmDouble W[NUMVERTICES];
     6813        IssmDouble BP[NUMVERTICES];
    59226814
    59236815        /*output: */
    … …  
    59256817        IssmDouble bslchydro=0;
    59266818        IssmDouble bslcbp=0;
     6819        IssmDouble BPavg=0;
     6820        IssmDouble Iavg=0;
     6821        IssmDouble Wavg=0;
    59276822
    59286823        /*ice properties: */
    59296824        IssmDouble rho_ice,rho_water,rho_freshwater;
    5930 
    5931         /*constants:*/
    5932         IssmDouble constant=0;
    5933        
    5934         /*recover parameters:*/
    5935         this->parameters->FindParam(&isice,TransientIsmasstransportEnum);
    5936         this->parameters->FindParam(&isocean,TransientIsoceantransportEnum);
    5937         this->parameters->FindParam(&ishydro,TransientIshydrologyEnum);
    5938 
    5939         /*early return if we are not on an ice cap, and we are not requesting
    5940          *hydrology or bottom pressure loads :*/
    5941         if(!isocean && !ishydro){
    5942                 if(!masks->isiceonly[this->lid]){
    5943                 #ifdef _ISSM_DEBUG_
    5944                         constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    5945                 #endif
    5946                 return;
    5947                 }
    5948         }
    5949 
    5950         /*early return if we are fully floating and we are not doing bottom pressure loads:*/
    5951         if(!isocean){
    5952                 if (masks->isfullyfloating[this->lid]){
    5953                         constant=0;
    5954                         #ifdef _ISSM_DEBUG_
    5955                         this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    5956                         #endif
    5957                         return;
    5958                 }
    5959         }
    5960 
    5961         /*early return if we are not on the ocean and we are not doing ice mass transport of
    5962          * hydrology:*/
    5963         if(!isice  && !ishydro){
    5964                 if(!masks->isoceanin[this->lid]){
    5965                         constant=0;
    5966                         #ifdef _ISSM_DEBUG_
    5967                         this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    5968                         #endif
    5969                         return;
    5970                 }
    5971 
    5972         }
    5973         /*if we are an ice shelf, are we fully grounded or not? (used later):*/
    5974         if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on.
    5975        
    59766825
    59776826        /*recover some parameters:*/
    … …  
    59796828        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    59806829        this->parameters->FindParam(&rho_freshwater,MaterialsRhoFreshwaterEnum);
    5981         this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    5982         this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
    5983         this->parameters->FindParam(&npartice,SolidearthNpartIceEnum);
    5984         this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
    5985 
    5986         /*Get area of element: precomputed in the sealevelchange_geometry:*/
    5987         this->Element::GetInputValue(&area,AreaEnum);
    5988 
    5989         /*Deal with ice loads if we are on grounded ice:*/
    5990         if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid] && isice){
    5991 
    5992                 /*Compute fraction of the element that is grounded: {{{*/
    5993                 if(notfullygrounded){
    5994                         IssmDouble xyz_list[NUMVERTICES][3];
    5995                         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5996 
    5997                         phi_ice=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
    5998                         if(glfraction==0)phi_ice=1;
    5999 #ifdef _ISSM_DEBUG_
    6000                         this->AddInput(SealevelBarystaticMaskEnum,&phi_ice,P0Enum);
    6001 #endif
    6002                 }
    6003                 else phi_ice=1.0;
    6004                 /*}}}*/
    6005 
    6006                 /*Inform mask: */
    6007                 constant+=100; //1 for ice.
    6008                 #ifdef _ISSM_DEBUG_
    6009                 this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    6010                 #endif
    6011 
    6012                 /*Retrieve surface load for ice: */
    6013                 Input* deltathickness_input=this->GetInput(DeltaIceThicknessEnum);
    6014                 if (!deltathickness_input)_error_("delta thickness input needed to compute sea level change!");
    6015 
    6016                 /*/Average ice thickness over grounded area of the element only: {{{*/
    6017                 if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
    6018                 else{
    6019                         IssmDouble total_weight=0;
    6020                         bool mainlyfloating = true;
    6021                         int         point1;
    6022                         IssmDouble  fraction1,fraction2;
    6023 
    6024                         /*Recover portion of element that is grounded*/
    6025                         this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    6026                         Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
    6027 
    6028                         /* Start  looping on the number of gaussian points and average over these gaussian points: */
    6029                         total_weight=0;
    6030                         I=0;
    6031                         while(gauss->next()){
    6032                                 IssmDouble Ig=0;
    6033                                 deltathickness_input->GetInputValue(&Ig,gauss);
    6034                                 I+=Ig*gauss->weight;
    6035                                 total_weight+=gauss->weight;
    6036                         }
    6037                         if(total_weight) I=I/total_weight;
    6038                         delete gauss;
    6039                 }
    6040                 /*}}}*/
    6041 
    6042                 /*Compute barystatic contribution in kg:*/
    6043                 bslcice = -rho_ice*phi_ice*area*I;
    6044                 _assert_(!xIsNan<IssmDouble>(bslcice));
    6045 
    6046                 /*Transfer thickness change into kg/m^2:*/
    6047                 I=I*rho_ice*phi_ice;
    6048         }
    6049 
    6050         /*Deal with water loads if we are on ground:*/
    6051         if(!masks->isfullyfloating[this->lid] && ishydro){
    6052 
    6053                 constant+=10; //1 for water.
    6054                 #ifdef _ISSM_DEBUG_
    6055                 this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    6056                 #endif
    6057 
    6058                 /*Retrieve water height at vertices: */
    6059                 Input* deltathickness_input=this->GetInput(DeltaTwsEnum);
    6060                 if (!deltathickness_input)_error_("DeltaTwsEnum input needed to compute sea level change!");
    6061                 deltathickness_input->GetInputAverage(&W);
    6062 
    6063                 /*Compute barystatic component in kg:*/
    6064                 bslchydro = -rho_freshwater*area*phi_water*W;
    6065                 _assert_(!xIsNan<IssmDouble>(bslchydro));
    6066 
    6067                 /*convert from m to kg/m^2:*/
    6068                 W=W*rho_freshwater*phi_water;
    6069         }
    6070 
    6071         /*Deal with ocean bottom pressures:*/
    6072         if(masks->isoceanin[this->lid] && isocean){
    6073 
    6074                 constant+=1; //1 for bottom pressure.
    6075                 #ifdef _ISSM_DEBUG_
    6076                 this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum);
    6077                 #endif
    6078 
    6079                 /*Retrieve bottom pressure change and average over the element: */
    6080                 Input* bottompressure_change_input=this->GetInput(DeltaBottomPressureEnum);
    6081                 if (!bottompressure_change_input)_error_("DeltaBottomPressureEnum pressure input needed to compute sea level change fingerprint!");
    6082                 bottompressure_change_input->GetInputAverage(&BP);
    6083                
    6084                 /*Compute barystatic component in kg:*/
    6085                 bslcbp = -rho_water*area*phi_bp*BP;
    6086 
    6087                 /*convert from m to kg/m^2:*/
    6088                 BP=BP*rho_water*phi_bp;
    6089         }
    6090 
    6091         /*Plug all loads into total load vector:*/
    6092         loads->SetValue(this->sid,I+W+BP,INS_VAL);
     6830
     6831        /*Retrieve inputs:*/
     6832        Element::GetInputListOnVertices(&I[0],DeltaIceThicknessEnum);
     6833        Element::GetInputListOnVertices(&W[0],DeltaTwsEnum);
     6834        Element::GetInputListOnVertices(&BP[0],DeltaBottomPressureEnum);
     6835
     6836        for(int i=0;i<NUMVERTICES;i++){
     6837                Iavg+=I[i]*slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid];
     6838                Wavg+=W[i]*slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid];
     6839                BPavg+=BP[i]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
     6840        }
     6841
     6842        /*convert from m to kg/m^2:*/
     6843        Iavg*=rho_ice;
     6844        Wavg*=rho_freshwater;
     6845        BPavg*=rho_water;
     6846
     6847        #ifdef _ISSM_DEBUG_
     6848        this->AddInput(SealevelBarystaticIceLoadEnum,&Iavg,P0Enum);
     6849        this->AddInput(SealevelBarystaticHydroLoadEnum,&Wavg,P0Enum);
     6850        this->AddInput(SealevelBarystaticBpLoadEnum,&BPavg,P0Enum);
     6851        #endif
     6852
     6853        /*Compute barystatic component in kg:*/
     6854        bslcice =   -slgeom->LoadArea[SLGEOM_ICE][this->lid]*Iavg;
     6855        bslchydro = -slgeom->LoadArea[SLGEOM_WATER][this->lid]*Wavg;
     6856        bslcbp =    -slgeom->LoadArea[SLGEOM_OCEAN][this->lid]*BPavg;
    60936857       
     6858        _assert_(!xIsNan<IssmDouble>(bslcice));
     6859        _assert_(!xIsNan<IssmDouble>(bslchydro));
     6860        _assert_(!xIsNan<IssmDouble>(bslcbp));
     6861       
     6862        /*Plug values into subelement load vector:*/
     6863        if(slgeom->issubelement[SLGEOM_ICE][this->lid]){
     6864                int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid];
     6865                subelementiceloads->SetValue(intj,Iavg,INS_VAL);
     6866                Iavg=0; //avoid double counting centroid loads and subelement loads!
     6867        }
     6868        if(slgeom->issubelement[SLGEOM_WATER][this->lid]){
     6869                int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid];
     6870                subelementhydroloads->SetValue(intj,Wavg,INS_VAL);
     6871                Wavg=0;
     6872        }
     6873        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
     6874                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
     6875                subelementbploads->SetValue(intj,BPavg,INS_VAL);
     6876                BPavg=0;
     6877        }
     6878        /*Plug remaining values into centroi load vector:*/
     6879        loads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL);
     6880
    60946881        /*Keep track of barystatic contributions:*/
    60956882        barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
    … …  
    60976884}
    60986885/*}}}*/
    6099 void       Tria::SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationvector){ /*{{{*/
     6886void       Tria::SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/
    61006887
    61016888        /*sal green function:*/
    61026889        IssmDouble* G=NULL;
     6890        IssmDouble* GsubelIce=NULL;
     6891        IssmDouble* GsubelHydro=NULL;
     6892        IssmDouble* GsubelOcean=NULL;
    61036893        IssmDouble SealevelGRD[NUMVERTICES]={0};
    6104         IssmDouble oceanaverage,oceanarea=0;
     6894        IssmDouble oceanaverage=0;
     6895        IssmDouble oceanarea=0;
    61056896        IssmDouble rho_water;
    61066897       
    … …  
    61086899        bool rotation= false;
    61096900        int  size;
    6110         int  nel;
     6901        int  nel,nbar;
    61116902        IssmDouble Grotm1[3];
    61126903        IssmDouble Grotm2[3];
    … …  
    61286919        if(sal){
    61296920                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
     6921                this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&GsubelIce,&size);
     6922                this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&GsubelHydro,&size);
     6923                this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&GsubelOcean,&size);
    61306924               
    61316925                if(allsealevelloads){
    … …  
    61346928                                        SealevelGRD[i]+=G[i*nel+e]*(allsealevelloads[e]+allloads[e]);
    61356929                                }
     6930                                nbar=slgeom->nbar[SLGEOM_ICE];
     6931                                for (int e=0;e<nbar;e++){
     6932                                        SealevelGRD[i]+=GsubelIce[i*nbar+e]*(allsubelementiceloads[e]);
     6933                                }
     6934                                nbar=slgeom->nbar[SLGEOM_WATER];
     6935                                for (int e=0;e<nbar;e++){
     6936                                        SealevelGRD[i]+=GsubelHydro[i*nbar+e]*(allsubelementhydroloads[e]);
     6937                                }
     6938                                nbar=slgeom->nbar[SLGEOM_OCEAN];
     6939                                for (int e=0;e<nbar;e++){
     6940                                        SealevelGRD[i]+=GsubelOcean[i*nbar+e]*(allsubelementbploads[e]+allsubelementsealevelloads[e]);
     6941                                }
    61366942                        }
    61376943                }
    … …  
    61416947                                        SealevelGRD[i]+=G[i*nel+e]*(allloads[e]);
    61426948                                }
    6143                         }
    6144                 }
    6145         }
    6146 
    6147         /*compute ocean average over element:*/
    6148         LevelsetAverage(&oceanaverage,&oceanarea,&SealevelGRD[0],MaskOceanLevelsetEnum);
    6149        
     6949                                nbar=slgeom->nbar[SLGEOM_ICE];
     6950                                for (int e=0;e<nbar;e++){
     6951                                        SealevelGRD[i]+=GsubelIce[i*nbar+e]*(allsubelementiceloads[e]);
     6952                                }
     6953                                nbar=slgeom->nbar[SLGEOM_WATER];
     6954                                for (int e=0;e<nbar;e++){
     6955                                        SealevelGRD[i]+=GsubelHydro[i*nbar+e]*(allsubelementhydroloads[e]);
     6956                                }
     6957                                nbar=slgeom->nbar[SLGEOM_OCEAN];
     6958                                for (int e=0;e<nbar;e++){
     6959                                        SealevelGRD[i]+=GsubelOcean[i*nbar+e]*(allsubelementbploads[e]);
     6960                                }
     6961                        }
     6962                }
     6963        }
     6964
     6965        /*retrieve ocean average and area:*/
     6966        for(int i=0;i<NUMVERTICES;i++){
     6967                oceanaverage+=SealevelGRD[i]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
     6968        }
     6969        #ifdef _ISSM_DEBUG_
     6970        this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum);
     6971        #endif
     6972        oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
     6973
    61506974        /*add ocean average in the global sealevelloads vector:*/
    6151         sealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
    6152         if(!allsealevelloads)oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
     6975        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
     6976                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
     6977                subelementsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL);
     6978        }
     6979        else sealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
     6980
     6981        if(!allsealevelloads){
     6982                oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
     6983                if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
     6984                        int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
     6985                        subelementoceanareas->SetValue(intj,oceanarea,INS_VAL);
     6986                }
     6987        }
     6988
    61536989       
    61546990        return;
    61556991} /*}}}*/
    6156 void       Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector){ /*{{{*/
     6992void       Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom){ /*{{{*/
    61576993
    61586994        IssmDouble Sealevel[3]={0,0,0};
    … …  
    61616997        IssmDouble SealevelN[3]={0,0,0};
    61626998        IssmDouble SealevelE[3]={0,0,0};
    6163         int nel;
     6999        int nel,nbar;
    61647000        bool sal = false;
    61657001        IssmDouble* G=NULL;
    … …  
    61677003        IssmDouble* GE=NULL;
    61687004        IssmDouble* GN=NULL;
     7005        IssmDouble* GsubelIce=NULL;
     7006        IssmDouble* GsubelHydro=NULL;
     7007        IssmDouble* GsubelOcean=NULL;
     7008        IssmDouble* GUsubelIce=NULL;
     7009        IssmDouble* GUsubelHydro=NULL;
     7010        IssmDouble* GUsubelOcean=NULL;
     7011        IssmDouble* GNsubelIce=NULL;
     7012        IssmDouble* GNsubelHydro=NULL;
     7013        IssmDouble* GNsubelOcean=NULL;
     7014        IssmDouble* GEsubelIce=NULL;
     7015        IssmDouble* GEsubelHydro=NULL;
     7016        IssmDouble* GEsubelOcean=NULL;
     7017
    61697018        int horiz;
    61707019        int size;
    … …  
    61917040        if(sal){
    61927041                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    6193                 if(elastic) this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size);
    6194                 if(horiz && elastic){
    6195                         this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size);
    6196                         this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size);
     7042                this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&GsubelIce,&size);
     7043                this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&GsubelHydro,&size);
     7044                this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&GsubelOcean,&size);
     7045                if(elastic){
     7046                        this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size);
     7047                        this->inputs->GetArrayPtr(SealevelchangeGUsubelIceEnum,this->lid,&GUsubelIce,&size);
     7048                        this->inputs->GetArrayPtr(SealevelchangeGUsubelHydroEnum,this->lid,&GUsubelHydro,&size);
     7049                        this->inputs->GetArrayPtr(SealevelchangeGUsubelOceanEnum,this->lid,&GUsubelOcean,&size);
     7050                        if(horiz){
     7051                                this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size);
     7052                                this->inputs->GetArrayPtr(SealevelchangeGEsubelIceEnum,this->lid,&GEsubelIce,&size);
     7053                                this->inputs->GetArrayPtr(SealevelchangeGEsubelOceanEnum,this->lid,&GEsubelOcean,&size);
     7054                                this->inputs->GetArrayPtr(SealevelchangeGEsubelHydroEnum,this->lid,&GEsubelHydro,&size);
     7055                                this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size);
     7056                                this->inputs->GetArrayPtr(SealevelchangeGNsubelIceEnum,this->lid,&GNsubelIce,&size);
     7057                                this->inputs->GetArrayPtr(SealevelchangeGNsubelHydroEnum,this->lid,&GNsubelHydro,&size);
     7058                                this->inputs->GetArrayPtr(SealevelchangeGNsubelOceanEnum,this->lid,&GNsubelOcean,&size);
     7059                        }
    61977060                }
    61987061
    … …  
    62017064                                SealevelRSL[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]);
    62027065                        }
     7066                        nbar=slgeom->nbar[SLGEOM_ICE];
     7067                        for (int e=0;e<nbar;e++){
     7068                                SealevelRSL[i]+=GsubelIce[i*nbar+e]*(subelementiceloads[e]);
     7069                        }
     7070                        nbar=slgeom->nbar[SLGEOM_WATER];
     7071                        for (int e=0;e<nbar;e++){
     7072                                SealevelRSL[i]+=GsubelHydro[i*nbar+e]*(subelementhydroloads[e]);
     7073                        }
     7074                        nbar=slgeom->nbar[SLGEOM_OCEAN];
     7075                        for (int e=0;e<nbar;e++){
     7076                                SealevelRSL[i]+=GsubelOcean[i*nbar+e]*(subelementbploads[e]+subelementsealevelloads[e]);
     7077                        }
     7078
    62037079                        if(elastic){
    62047080                                for (int e=0;e<nel;e++){
    62057081                                        SealevelU[i]+=GU[i*nel+e]*(sealevelloads[e]+loads[e]);
    62067082                                }
     7083                                nbar=slgeom->nbar[SLGEOM_ICE];
     7084                                for (int e=0;e<nbar;e++){
     7085                                        SealevelU[i]+=GUsubelIce[i*nbar+e]*(subelementiceloads[e]);
     7086                                }
     7087                                nbar=slgeom->nbar[SLGEOM_WATER];
     7088                                for (int e=0;e<nbar;e++){
     7089                                        SealevelU[i]+=GUsubelHydro[i*nbar+e]*(subelementhydroloads[e]);
     7090                                }
     7091                                nbar=slgeom->nbar[SLGEOM_OCEAN];
     7092                                for (int e=0;e<nbar;e++){
     7093                                        SealevelU[i]+=GUsubelOcean[i*nbar+e]*(subelementbploads[e]+subelementsealevelloads[e]);
     7094                                }
    62077095                        }
    62087096                        if(horiz && elastic){
    … …  
    62107098                                        SealevelN[i]+=GN[i*nel+e]*(sealevelloads[e]+loads[e]);
    62117099                                        SealevelE[i]+=GE[i*nel+e]*(sealevelloads[e]+loads[e]);
     7100                                }
     7101                                nbar=slgeom->nbar[SLGEOM_ICE];
     7102                                for (int e=0;e<nbar;e++){
     7103                                        SealevelN[i]+=GNsubelIce[i*nbar+e]*(subelementiceloads[e]);
     7104                                        SealevelE[i]+=GEsubelIce[i*nbar+e]*(subelementiceloads[e]);
     7105                                }
     7106                                nbar=slgeom->nbar[SLGEOM_WATER];
     7107                                for (int e=0;e<nbar;e++){
     7108                                        SealevelN[i]+=GNsubelHydro[i*nbar+e]*(subelementhydroloads[e]);
     7109                                        SealevelE[i]+=GEsubelHydro[i*nbar+e]*(subelementhydroloads[e]);
     7110                                }
     7111                                nbar=slgeom->nbar[SLGEOM_OCEAN];
     7112                                for (int e=0;e<nbar;e++){
     7113                                        SealevelN[i]+=GNsubelOcean[i*nbar+e]*(subelementbploads[e]+subelementsealevelloads[e]);
     7114                                        SealevelE[i]+=GEsubelOcean[i*nbar+e]*(subelementbploads[e]+subelementsealevelloads[e]);
    62127115                                }
    62137116                        }
    … …  
    62287131
    62297132} /*}}}*/
    6230 void       Tria::LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum){ /*{{{*/
    6231 
    6232         IssmDouble phi=1.0;
    6233         IssmDouble area;
    6234         IssmDouble average = 0;
    6235         IssmDouble total_weight=0;
    6236         int        point1;
    6237         IssmDouble fraction1,fraction2;
    6238         IssmDouble levelsetvalues[NUMVERTICES];
    6239         bool       fractiongeometryflag=true;
    6240        
    6241 
    6242 
    6243         /*retrieve value of levelset:*/
    6244         Input *levelset= this->GetInput(levelsetenum);
    6245         this->Element::GetInputListOnVertices(&levelsetvalues[0],levelsetenum);
    6246        
    6247         /*Early return if no vertices are inside the desired area:*/
    6248         if(levelset->GetInputMin()>=0){
    6249                 *paverage=0;
    6250                 *parea=0;
    6251                 return;
    6252         }
    6253 
    6254         /*Get area of element:*/
    6255         this->Element::GetInputValue(&area,AreaEnum);
    6256 
    6257         /*Are we fully in the desired area, in which case average over the entire area?:*/
    6258         if(levelset->GetInputMax()<=0){
    6259                 for(int i=0;i<NUMVERTICES;i++) average+=field_on_localvertices[i]/NUMVERTICES;
    6260                
    6261                 *parea=area;
    6262                 *paverage=average;
    6263                 return;
    6264         }
    6265 
    6266         /*What fraction of the triangle is in the desired area?:*/
    6267         IssmDouble xyz_list[NUMVERTICES][3];
    6268         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6269         phi=this->GetFractionArea(&xyz_list[0][0],levelsetvalues);
    6270         area=phi*area;
    6271 
    6272         /*Average over  the fraction area only, using the right gaussian points: */
    6273         this->GetFractionGeometry(&point1,&fraction1,&fraction2,&fractiongeometryflag,levelsetvalues);
    6274         Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,fractiongeometryflag,2);
    6275 
    6276         total_weight=0;
    6277         average=0;
    6278         while(gauss->next()){
    6279                 IssmDouble field_gauss=0;
    6280                 TriaRef::GetInputValue(&field_gauss, field_on_localvertices, gauss,P1Enum);
    6281                 average+=field_gauss*gauss->weight;
    6282                 total_weight+=gauss->weight;
    6283         }
    6284         if(total_weight) average=average/total_weight;
    6285        
    6286         /*free ressources:*/
    6287         delete gauss;
    6288 
    6289         *paverage=average;
    6290         *parea=area;
    6291         return;
    6292 
    6293 }
    6294 /*}}}*/
    6295 void       Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){/*{{{*/
     7133void       Tria::SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom){/*{{{*/
    62967134               
    62977135        IssmDouble S=0;
    62987136
    62997137        /*Compute area of element:*/
    6300         IssmDouble area,planetarea;
     7138        IssmDouble area,planetarea,re;
     7139        IssmDouble late,longe;
    63017140        this->Element::GetInputValue(&area,AreaEnum);
    63027141
    6303         /*recover earth area: */
     7142        /*recover parameters: */
    63047143        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    6305 
    6306         /*Compute lat,long,radius of elemental centroid: */
    6307         bool spherical=true;
    6308         IssmDouble llr_list[NUMVERTICES][3];
    6309         IssmDouble late,longe,re;
    6310         /* Where is the centroid of this element?:{{{*/
    6311         ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
    6312 
    6313         IssmDouble minlong=400;
    6314         IssmDouble maxlong=-20;
    6315         for (int i=0;i<NUMVERTICES;i++){
    6316                 llr_list[i][0]=(90-llr_list[i][0]);
    6317                 if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
    6318                 if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
    6319                 if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
    6320         }
    6321         if(minlong==0 && maxlong>180){
    6322                 if (llr_list[0][1]==0)llr_list[0][1]=360;
    6323                 if (llr_list[1][1]==0)llr_list[1][1]=360;
    6324                 if (llr_list[2][1]==0)llr_list[2][1]=360;
    6325         }
    6326 
    6327         // correction at the north pole
    6328         if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    6329         if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    6330         if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
    6331 
    6332         //correction at the south pole
    6333         if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    6334         if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    6335         if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
    6336 
    6337         late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
    6338         longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
    6339 
    6340         late=90.-late;
    6341         if(longe>180.)longe=(longe-180.)-180.;
    6342 
    6343         late=late/180.*M_PI;
    6344         longe=longe/180.*M_PI;
    6345         /*}}}*/
    6346         re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0;
     7144        this->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);
     7145        late=slgeom->late[this->lid]/180*M_PI;
     7146        longe=slgeom->longe[this->lid]/180*M_PI;
    63477147
    63487148
    … …  
    63507150        if(loads) S+=loads[this->Sid()];
    63517151        if(sealevelloads) S+=sealevelloads[this->Sid()];
    6352 
     7152       
    63537153        /* Perturbation terms for moment of inertia (moi_list):
    63547154         * computed analytically (see Wu & Peltier, eqs 10 & 32)
    … …  
    63617161        return;
    63627162}/*}}}*/
    6363 void       Tria::SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks){ /*{{{*/
     7163void       Tria::SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){/*{{{*/
     7164               
     7165        IssmDouble  SA=0;
     7166        IssmDouble* loads=NULL;
     7167        IssmDouble* sealevelloads=NULL;
     7168        IssmDouble  late,longe,re;
     7169        int         intj;
     7170        IssmDouble  area;
     7171        IssmDouble  planetarea;
     7172
     7173        /*recover parameters: */
     7174        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
     7175        this->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);
    63647176       
    6365         if(masks->isoceanin[this->lid]){
    6366                 loads->SetValue(this->sid,offset,ADD_VAL);
    6367         }
     7177        /*Initalize:*/
     7178        for(int i=0;i<3;i++)dI_list[i]=0;
     7179
     7180        /*Go through our loads:*/
     7181        for(int i=0;i<SLGEOM_NUMLOADS;i++){
     7182                if(slgeom->issubelement[i][this->lid]){
     7183                        switch(i){
     7184                                case SLGEOM_ICE:
     7185                                        loads=subelementiceloads;
     7186                                        break;
     7187                                case SLGEOM_WATER:
     7188                                        loads=subelementhydroloads;
     7189                                        break;
     7190                                case SLGEOM_OCEAN:
     7191                                        loads=subelementbploads;
     7192                                        sealevelloads=subelementsealevelloads;
     7193                                        break;
     7194                        }
     7195                        intj=slgeom->subelementmapping[i][this->lid];
     7196                        late=slgeom->latbarycentre[i][intj]/180*M_PI;
     7197                        longe=slgeom->longbarycentre[i][intj]/180*M_PI;
     7198                        area=slgeom->area_subel[i][intj];
     7199
     7200                        /*recover total load: */
     7201                        if(loads) SA+=loads[intj]*area;
     7202                        if(sealevelloads) SA+=sealevelloads[intj]*area;
     7203                }
     7204        }
     7205
     7206        /* Perturbation terms for moment of inertia (moi_list):
     7207         * computed analytically (see Wu & Peltier, eqs 10 & 32)
     7208         * also consistent with my GMD formulation!
     7209         * ALL in geographic coordinates
     7210         * */
     7211        dI_list[0] += -4*M_PI*(SA)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
     7212        dI_list[1] += -4*M_PI*(SA)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
     7213        dI_list[2] += +4*M_PI*(SA)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
     7214
     7215        return;
     7216}/*}}}*/
     7217void       Tria::SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
     7218       
     7219        if(slgeom->isoceanin[this->lid]){
     7220                if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
     7221                        int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
     7222                        subelementloads->SetValue(intj,offset,ADD_VAL);
     7223                }
     7224                else loads->SetValue(this->sid,offset,ADD_VAL);
     7225        }
     7226
    63687227
    63697228} /*}}}*/

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

    r26167 r26222  
    7171                void        CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum);
    7272                void        ElementCoordinates(Vector<IssmDouble>* vxe,Vector<IssmDouble>* vye,Vector<IssmDouble>* vze, Vector<IssmDouble>* vareae, bool spherical=false);
     73                void        ElementCoordinates(Vector<IssmDouble>* vlonge,Vector<IssmDouble>* vlate,Vector<IssmDouble>* vareae);
    7374                int         EdgeOnBaseIndex();
    7475                void        EdgeOnBaseIndices(int* pindex1,int* pindex);
    … …  
    8485                void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
    8586                IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
    86                 void        GetFractionGeometry(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* pmainlynegative, IssmDouble* gl);
    87                 IssmDouble  GetFractionArea(IssmDouble* xyz_list, IssmDouble* gl);
     87                void        GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl);
     88                void        GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelsetenum);
     89                void        GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum);
     90                void        GetBarycenterFromLevelset(IssmDouble* platbar, IssmDouble* plongbar,IssmDouble phi,IssmDouble fraction1,IssmDouble fraction2,IssmDouble late, IssmDouble longe, int point1,int istrapeze1, IssmDouble planetradius);
     91                IssmDouble  GetTriangleAreaSpherical(IssmDouble xyz_list[3][3]);
    8892                IssmDouble  GetIcefrontArea();
    8993                void          GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
    … …  
    164168                #endif
    165169                #ifdef _HAVE_SEALEVELCHANGE_
    166                 void       LevelsetAverage(IssmDouble* paverage, IssmDouble* parea, IssmDouble* field_on_localvertices, int levelsetenum);
    167                 void       SealevelchangeMomentOfInertia(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads);
    168                 void       SealevelchangeGeometry(IssmDouble* lat, IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
    169170                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
    170                 void       SetSealevelMasks(SealevelMasks* masks);
    171                 void       SealevelchangeGeometry(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
    172                 void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks);
    173                 void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* oceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* rotationaxismotionvector);
    174                 void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector);
    175                 void       SealevelchangeShift(Vector<IssmDouble>* loads, IssmDouble offset, SealevelMasks* masks);
     171                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
     172                void       SealevelchangeGeometryFractionKernel(SealevelGeometry* slgeom);
     173                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae);
     174                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
     175                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementiceloads, Vector<IssmDouble>* subelementhydroloads, Vector<IssmDouble>* subelementbploads, BarystaticContributions* barycontrib, SealevelGeometry* slgeom);
     176                void       SealevelchangeConvolution(Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* allsealevelloads, IssmDouble* allloads,IssmDouble* allsubelementiceloads, IssmDouble* allsubelementhydroloads, IssmDouble* allsubelementbploads, IssmDouble* allsubelementsealevelloads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
     177                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
     178                void       SealevelchangeShift(Vector<IssmDouble>* loads, Vector<IssmDouble>* subelementloads, IssmDouble offset, SealevelGeometry* slgeom);
     179                void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads, SealevelGeometry* slgeom);
     180                void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom);
    176181                #endif
    177182                /*}}}*/

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

    r26110 r26222  
    1818#include "./Massfluxatgate.h"
    1919#include "./Misfit.h"
    20 #include "./SealevelMasks.h"
     20#include "./SealevelGeometry.h"
    2121#include "./BarystaticContributions.h"
    2222#include "./Nodalvalue.h"

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

    r26165 r26222  
    99class FemModel;
    1010class Parameters;
    11 class SealevelMasks;
     11class SealevelGeometry;
    1212template <class doubletype> class Matrix;
    1313template <class doubletype> class Vector;
    … …  
    4545void steadystate_core(FemModel* femmodel);
    4646void transient_core(FemModel* femmodel);
     47void transient_precore(FemModel* femmodel);
    4748void dakota_core(FemModel* femmodel);
    4849void ad_core(FemModel* femmodel);
    … …  
    6061#ifdef _HAVE_SEALEVELCHANGE_
    6162void sealevelchange_core(FemModel* femmodel);
    62 void sealevelchange_geometry(FemModel* femmodel);
     63void sealevelchange_initialgeometry(FemModel* femmodel);
     64SealevelGeometry* sealevelchange_geometry(FemModel* femmodel);
    6365#endif
    64 void grd_core(FemModel* femmodel);
     66void grd_core(FemModel* femmodel,SealevelGeometry* slgeom);
    6567void solidearthexternal_core(FemModel* femmodel);
    6668void dynstr_core(FemModel* femmodel);
    … …  
    8284void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype);
    8385void WrapperCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype,bool nodakotacore=false);
     86void WrapperPreCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype);
    8487void AdjointCorePointerFromSolutionEnum(void (**padjointcore)(FemModel*),int solutiontype);
    8588

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

    r26165 r26222  
    22 * \brief: core of the sea-level change solution
    33 */
     4
     5#ifdef HAVE_CONFIG_H
     6        #include <config.h>
     7#else
     8#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
     9#endif
     10
    411
    512#include "./cores.h"
    … …  
    1724void TransferSealevel(FemModel* femmodel,int forcingenum);
    1825void slcconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
    19 IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea);
    20 SealevelMasks* sealevel_masks(FemModel* femmodel);
    21 void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealeveloads);
    22 void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, IssmDouble offset, SealevelMasks* masks);
     26IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble oceanarea);
     27void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom);
     28void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, IssmDouble offset, SealevelGeometry* slgeom);
    2329void ivins_deformation_core(FemModel* femmodel);
     30IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel);
    2431/*}}}*/
    2532
    … …  
    2734void sealevelchange_core(FemModel* femmodel){ /*{{{*/
    2835
     36        SealevelGeometry* slgeom=NULL;
     37
    2938        /*Start profiler*/
    3039        femmodel->profiler->Start(SLRCORE);
    … …  
    4352
    4453        /*run geometry core: */
    45         sealevelchange_geometry(femmodel);
     54        slgeom=sealevelchange_geometry(femmodel);
    4655       
    4756        /*any external forcings?:*/
    … …  
    5261
    5362        /*Run geodetic:*/
    54         grd_core(femmodel);
     63        grd_core(femmodel,slgeom);
    5564
    5665        /*Run steric core for sure:*/
    … …  
    206215       
    207216}; /*}}}*/
    208 void grd_core(FemModel* femmodel) { /*{{{*/
     217void grd_core(FemModel* femmodel, SealevelGeometry* slgeom) { /*{{{*/
    209218
    210219        /*variables:{{{*/
    211220        int nel;
    212221        BarystaticContributions* barycontrib=NULL;
    213         SealevelMasks* masks=NULL;
    214222        GenericParam<BarystaticContributions*>* barycontribparam=NULL;
    215223        IssmDouble rotationaxismotionvector[3]={0};
    … …  
    217225        Vector<IssmDouble>*    loads=NULL;
    218226        IssmDouble*            allloads=NULL;
     227        Vector<IssmDouble>*    subelementiceloads=NULL;
     228        IssmDouble*            allsubelementiceloads=NULL;
     229        Vector<IssmDouble>*    subelementhydroloads=NULL;
     230        IssmDouble*            allsubelementhydroloads=NULL;
     231        Vector<IssmDouble>*    subelementbploads=NULL;
     232        IssmDouble*            allsubelementbploads=NULL;
    219233        Vector<IssmDouble>*    sealevelloads=NULL;
     234        IssmDouble*            allsealevelloads=NULL;
     235        Vector<IssmDouble>*    subelementsealevelloads=NULL;
     236        IssmDouble*            allsubelementsealevelloads=NULL;
    220237        Vector<IssmDouble>*    oldsealevelloads=NULL;
    221         IssmDouble*            allsealevelloads=NULL;
    222238        Vector<IssmDouble>*    oceanareas=NULL;
    223239        IssmDouble             oceanarea;
     240        Vector<IssmDouble>*    subelementoceanareas=NULL;
    224241        IssmDouble             oceanaverage;
    225242        bool                   scaleoceanarea=false;
    … …  
    265282                if(modelid!=earthid)return;
    266283        }
    267 
    268284        /*branch directly to Ivins deformation core if requested:*/
    269285        if(grdmodel==IvinsEnum){
    … …  
    282298        /*}}}*/
    283299
    284         /*initialize matrices and vectors:*/
     300        /*initialize loads and sea level loads:*/
    285301        femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
     302
    286303        loads=new Vector<IssmDouble>(nel);
    287         oceanareas=new Vector<IssmDouble>(nel);
     304        subelementiceloads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_ICE]);
     305        subelementhydroloads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_WATER]);
     306        subelementbploads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
     307
    288308        sealevelloads=new Vector<IssmDouble>(nel);
    289309        sealevelloads->Set(0);sealevelloads->Assemble();
    290 
    291         /*call masks core: */
    292         masks=sealevel_masks(femmodel);
     310        subelementsealevelloads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
     311        subelementoceanareas=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
     312        oceanareas=new Vector<IssmDouble>(nel);
     313
    293314        if(VerboseSolution()) _printf0_("         starting  GRD convolutions\n");
    294        
     315
    295316        /*buildup loads: */
    296317        for(Object* & object : femmodel->elements->objects){
    297318                Element* element = xDynamicCast<Element*>(object);
    298                 element->SealevelchangeBarystaticLoads(loads, barycontrib,masks);
     319                element->SealevelchangeBarystaticLoads(loads, subelementiceloads, subelementhydroloads, subelementbploads, barycontrib,slgeom);
    299320        }
    300321        loads->Assemble();
     322        subelementiceloads->Assemble();
     323        subelementhydroloads->Assemble();
     324        subelementbploads->Assemble();
     325       
    301326
    302327        //broadcast loads
    303328        allloads=loads->ToMPISerial();
     329        allsubelementiceloads=subelementiceloads->ToMPISerial();
     330        allsubelementhydroloads=subelementhydroloads->ToMPISerial();
     331        allsubelementbploads=subelementbploads->ToMPISerial();
    304332
    305333        //compute rotation axis motion:
    306         RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,NULL);
     334        RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsubelementiceloads,allsubelementhydroloads,allsubelementbploads, NULL, NULL,slgeom);
    307335
    308336        /*skip computation of sea level if requested, which means sea level loads should be zeroed */
    309337        if(!computesealevel){
    310338                allsealevelloads=xNewZeroInit<IssmDouble>(nel);
     339                allsubelementsealevelloads=xNewZeroInit<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
    311340                goto deformation;
    312341        }
    … …  
    314343        if(VerboseSolution()) _printf0_("         converging GRD convolutions\n");
    315344        for(;;){
    316                        
     345
    317346                oldsealevelloads=sealevelloads->Duplicate(); sealevelloads->Copy(oldsealevelloads);
    318347
    … …  
    320349                for(Object* & object : femmodel->elements->objects){
    321350                        Element* element = xDynamicCast<Element*>(object);
    322                         element->SealevelchangeConvolution(sealevelloads, oceanareas, allsealevelloads, allloads,rotationaxismotionvector);
    323                 }
     351                        element->SealevelchangeConvolution(sealevelloads, subelementsealevelloads, oceanareas, subelementoceanareas, allsealevelloads, allloads,allsubelementiceloads, allsubelementhydroloads, allsubelementbploads, allsubelementsealevelloads, rotationaxismotionvector,slgeom);
     352                }
     353
    324354                sealevelloads->Assemble();
     355                subelementsealevelloads->Assemble();
    325356
    326357                /*compute ocean areas:*/
    327358                if(!allsealevelloads){ //first time in the loop
    328                         oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
     359                        oceanareas->Assemble();
     360                        subelementoceanareas->Assemble();
     361                        oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
    329362                        if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
    330363                }
     364       
     365                //Conserve ocean mass:
     366                oceanaverage=SealevelloadsOceanAverage(sealevelloads,subelementsealevelloads, oceanareas,subelementoceanareas, oceanarea);
    331367               
    332                 //Conserve ocean mass:
    333                 oceanaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea);
    334                 ConserveOceanMass(femmodel,sealevelloads,barycontrib->Total()/oceanarea - oceanaverage,masks);
     368                ConserveOceanMass(femmodel,sealevelloads,subelementsealevelloads,barycontrib->Total()/oceanarea - oceanaverage,slgeom);
    335369
    336370                //broadcast sea level loads
    337371                allsealevelloads=sealevelloads->ToMPISerial();
     372                allsubelementsealevelloads=subelementsealevelloads->ToMPISerial();
    338373
    339374                //compute rotation axis motion:
    340                 RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads);
     375                RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsubelementiceloads,allsubelementhydroloads,allsubelementbploads, allsealevelloads,allsubelementsealevelloads,slgeom);
    341376
    342377                //convergence?
    … …  
    356391        for(Object* & object : femmodel->elements->objects){
    357392                Element* element = xDynamicCast<Element*>(object);
    358                 element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector);
    359         }
    360        
     393                element->SealevelchangeDeformationConvolution(allsealevelloads, allsubelementsealevelloads,allloads, allsubelementiceloads, allsubelementhydroloads, allsubelementbploads,rotationaxismotionvector,slgeom);
     394        }
     395
    361396        if(VerboseSolution()) _printf0_("         updating GRD fields\n");
    362397
    363398        /*Update bedrock motion and geoid:*/
    364399        if(computesealevel){
    365                 femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water);
    366 
    367                 /*cumulate barystatic contributions and save to results: */
     400                femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/rho_water/oceanarea- oceanaverage/rho_water); //given that we converged, no need to recompute ocean average
     401
     402                //cumulate barystatic contributions and save to results:
    368403                barycontrib->Cumulate(femmodel->parameters);
    369404                barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
    … …  
    518553
    519554//Geometry:
    520 SealevelMasks* sealevel_masks(FemModel* femmodel) {  /*{{{*/
    521 
    522         int grdmodel=0;
    523 
    524         /*early return?:*/
    525         femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
    526         if(grdmodel==IvinsEnum) return NULL;
    527 
    528         if(VerboseSolution()) _printf0_("         computing sea level masks\n");
    529        
    530         /*initialize SealevelMasks structure: */
    531         SealevelMasks* masks=new SealevelMasks(femmodel->elements->Size());
    532 
    533         /*go through elements and fill the masks: */
    534         for(Object* & object : femmodel->elements->objects){
    535                 Element*   element=xDynamicCast<Element*>(object);
    536                 element->SetSealevelMasks(masks);
    537         }
    538 
    539         return masks;
    540 }/*}}}*/
    541 void sealevelchange_geometry(FemModel* femmodel) {  /*{{{*/
    542 
    543         /*Geometry core where we compute indices into tables pre computed in the SealevelRiseAnalysis: */
     555void sealevelchange_initialgeometry(FemModel* femmodel) {  /*{{{*/
     556
     557        /*Geometry core where we compute geometrical kernels and weights:*/
    544558
    545559        /*parameters: */
    546         bool spherical=true;
    547         IssmDouble *latitude  = NULL;
    548         IssmDouble *longitude = NULL;
    549         IssmDouble *radius    = NULL;
    550         IssmDouble *xx    = NULL;
    551         IssmDouble *yy    = NULL;
    552         IssmDouble *zz    = NULL;
    553560        IssmDouble *xxe    = NULL;
    554561        IssmDouble *yye    = NULL;
    555562        IssmDouble *zze    = NULL;
    556563        IssmDouble* areae  = NULL;
    557 
    558         int  horiz;
    559         bool geometrydone = false;
    560         int  optim;
     564        int  nel;
    561565        int  grdmodel=0;
    562566
    563                
    564567        /*retrieve parameters:*/
    565         femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    566         femmodel->parameters->FindParam(&geometrydone,SealevelchangeGeometryDoneEnum);
    567568        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
     569        nel=femmodel->elements->NumberOfElements();
    568570       
    569571        /*early return?:*/
    570572        if(grdmodel==IvinsEnum) return;
    571573
    572         if(geometrydone){
    573                 if(VerboseSolution()) _printf0_("         geometrical offsets have already been computed, skipping \n");
    574                 return; //don't need to run this again.
    575         }
    576 
    577574        /*Verbose: */
    578         if(VerboseSolution()) _printf0_("         computing geometrical offsets into precomputed Green tables \n");
    579 
    580         /*first, recover lat,long and radius vectors from vertices: */
    581         VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
    582         if(horiz) VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
    583        
    584         /*first, recover x,y,z and areas from elements: */
     575        if(VerboseSolution()) _printf0_("         computing initial sea level geometrical kernels and weights.\n");
     576
     577        /*recover x,y,z and areas from elements: */
    585578        ElementCoordinatesx(&xxe,&yye,&zze,&areae,femmodel->elements);
    586579
    … …  
    588581        for(Object* & object : femmodel->elements->objects){
    589582                Element*   element=xDynamicCast<Element*>(object);
    590                 element->SealevelchangeGeometry(latitude,longitude,radius,xx,yy,zz,xxe,yye,zze,areae);
    591         }
    592 
    593         /*Free ressources:*/
    594         if(horiz){
    595                 xDelete<IssmDouble>(xx);
    596                 xDelete<IssmDouble>(yy);
    597                 xDelete<IssmDouble>(zz);
    598         }
    599         xDelete<IssmDouble>(xxe);
    600         xDelete<IssmDouble>(yye);
    601         xDelete<IssmDouble>(zze);
    602         xDelete<IssmDouble>(areae);
    603         xDelete<IssmDouble>(latitude);
    604         xDelete<IssmDouble>(longitude);
    605         xDelete<IssmDouble>(radius);
    606 
    607         /*Record the fact that we ran this module already: */
    608         femmodel->parameters->SetParam(true,SealevelchangeGeometryDoneEnum);
    609 
     583                element->SealevelchangeGeometryInitial(xxe,yye,zze,areae);
     584        }
     585
     586        femmodel->parameters->AddObject(new DoubleVecParam(XxeEnum,xxe,nel));
     587        femmodel->parameters->AddObject(new DoubleVecParam(YyeEnum,yye,nel));
     588        femmodel->parameters->AddObject(new DoubleVecParam(ZzeEnum,zze,nel));
     589        femmodel->parameters->AddObject(new DoubleVecParam(AreaeEnum,areae,nel));
     590
     591
     592        #ifdef _ISSM_DEBUG_
     593        femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,XxeEnum,xxe,nel,1,1,1));
     594        femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,YyeEnum,yye,nel,1,1,1));
     595        femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,ZzeEnum,zze,nel,1,1,1));
     596        femmodel->results->AddResult(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,AreaeEnum,areae,nel,1,1,1));
     597        #endif
     598
     599        return;
     600
     601
     602}/*}}}*/
     603SealevelGeometry* sealevelchange_geometry(FemModel* femmodel) {  /*{{{*/
     604
     605        /*Geometry core where we compute updates to the Green function kernels and weights, dependent
     606         * on the evolution of levelsets: */
     607
     608        /*parameters: */
     609        IssmDouble *xxe    = NULL;
     610        IssmDouble *yye    = NULL;
     611        IssmDouble *zze    = NULL;
     612        IssmDouble* areae  = NULL;
     613
     614        int nel;
     615        int  grdmodel=0;
     616        SealevelGeometry* slgeom=NULL;
     617
     618        /*early return?:*/
     619        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
     620        if(grdmodel==IvinsEnum) return NULL;
     621
     622        /*retrieve parameters:*/
     623        femmodel->parameters->FindParam(&xxe,&nel,XxeEnum);
     624        femmodel->parameters->FindParam(&yye,&nel,YyeEnum);
     625        femmodel->parameters->FindParam(&zze,&nel,ZzeEnum);
     626        femmodel->parameters->FindParam(&areae,&nel,AreaeEnum);
     627       
     628        /*initialize SealevelMasks structure: */
     629        slgeom=new SealevelGeometry(femmodel->elements->Size());
     630       
     631        /*Verbose: */
     632        if(VerboseSolution()) _printf0_("         computing sea level geometrical kernel and weight updates.\n");
     633       
     634       
     635        /*Run sealevel geometry routine for elements with full loading:*/
     636        for(Object* & object : femmodel->elements->objects){
     637                Element*   element=xDynamicCast<Element*>(object);
     638                element->SealevelchangeGeometryCentroidLoads(slgeom,xxe,yye,zze,areae);
     639        }
     640       
     641        /*Initialize fractional loading mapping: */
     642        slgeom->InitializeMappingsAndBarycentres();
     643
     644       
     645        /*Run sealevel geometry routine for elements with fractional loading:*/
     646        for(Object* & object : femmodel->elements->objects){
     647                Element*   element=xDynamicCast<Element*>(object);
     648                element->SealevelchangeGeometrySubElementLoads(slgeom,areae);
     649        }
     650
     651        /*Assemble barycentres of fraction loading elements:*/
     652        slgeom->Assemble();
     653
     654        /*Create fractional green function kernels: */
     655        for(Object* & object : femmodel->elements->objects){
     656                Element*   element=xDynamicCast<Element*>(object);
     657                element->SealevelchangeGeometryFractionKernel(slgeom);
     658        }
     659
     660
     661        femmodel->parameters->AddObject(new DoubleVecParam(XxeEnum,xxe,nel));
     662        femmodel->parameters->AddObject(new DoubleVecParam(YyeEnum,yye,nel));
     663        femmodel->parameters->AddObject(new DoubleVecParam(ZzeEnum,zze,nel));
     664        femmodel->parameters->AddObject(new DoubleVecParam(AreaeEnum,areae,nel));
     665
     666        return slgeom;
    610667
    611668}/*}}}*/
    … …  
    903960
    904961} /*}}}*/
    905 IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea){ /*{{{*/
     962IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* subelementsealevelloads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble oceanarea){ /*{{{*/
    906963
    907964        IssmDouble sealevelloadsaverage;       
     965        IssmDouble subelementsealevelloadsaverage;     
    908966
    909967        Vector<IssmDouble>* sealevelloadsvolume=sealevelloads->Duplicate();
     968        Vector<IssmDouble>* subelementsealevelloadsvolume=subelementsealevelloads->Duplicate();
     969
    910970        sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas);
     971        subelementsealevelloadsvolume->PointwiseMult(subelementsealevelloads,subelementoceanareas);
     972       
    911973        sealevelloadsvolume->Sum(&sealevelloadsaverage);
     974        subelementsealevelloadsvolume->Sum(&subelementsealevelloadsaverage);
    912975        delete sealevelloadsvolume;
    913        
    914         return sealevelloadsaverage/oceanarea;
     976        delete subelementsealevelloadsvolume;
     977       
     978        //return (sealevelloadsaverage+subelementsealevelloadsaverage)/oceanarea;
     979        return (sealevelloadsaverage)/oceanarea;
    915980} /*}}}*/
    916 void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealevelloads){ /*{{{*/
     981void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* subelementiceloads, IssmDouble* subelementhydroloads, IssmDouble* subelementbploads, IssmDouble* sealevelloads, IssmDouble* subelementsealevelloads, SealevelGeometry* slgeom){ /*{{{*/
    917982
    918983        IssmDouble  moi_list[3]={0,0,0};
     984        IssmDouble  moi_list_sub[3]={0,0,0};
    919985        IssmDouble  moi_list_cpu[3]={0,0,0};
    920986        IssmDouble*     tide_love_h  = NULL;
    … …  
    9401006        for(Object* & object : femmodel->elements->objects){
    9411007                Element* element = xDynamicCast<Element*>(object);
    942                 element->SealevelchangeMomentOfInertia(&moi_list[0],loads,sealevelloads);
    943                 moi_list_cpu[0] += moi_list[0];
    944                 moi_list_cpu[1] += moi_list[1];
    945                 moi_list_cpu[2] += moi_list[2];
    946         }
     1008
     1009                element->SealevelchangeMomentOfInertiaCentroid(&moi_list[0],loads,sealevelloads,slgeom);
     1010
     1011                element->SealevelchangeMomentOfInertiaSubElement(&moi_list_sub[0],subelementiceloads, subelementhydroloads, subelementbploads, subelementsealevelloads,slgeom);
     1012
     1013                moi_list_cpu[0] += moi_list[0]+moi_list_sub[0];
     1014                moi_list_cpu[1] += moi_list[1]+moi_list_sub[1];
     1015                moi_list_cpu[2] += moi_list[2]+moi_list_sub[2];
     1016        }
     1017
    9471018
    9481019        ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    … …  
    9651036        m[2]=m3;
    9661037} /*}}}*/
    967 void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, IssmDouble offset, SealevelMasks* masks){ /*{{{*/
     1038void ConserveOceanMass(FemModel* femmodel,Vector<IssmDouble>* sealevelloads, Vector<IssmDouble>* subelementsealevelloads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
    9681039
    9691040        /*Shift sealevel loads by ocean average, only on ocean! :*/
    9701041        for(Object* & object : femmodel->elements->objects){
    9711042                Element* element = xDynamicCast<Element*>(object);
    972                 element->SealevelchangeShift(sealevelloads,offset,masks);
     1043                element->SealevelchangeShift(sealevelloads,subelementsealevelloads,offset,slgeom);
    9731044        }
    9741045        sealevelloads->Assemble();
     1046        subelementsealevelloads->Assemble();
    9751047
    9761048} /*}}}*/
     1049IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel){ /*{{{*/
     1050
     1051        int* indices=xNew<int>(nel);
     1052        for(int i=0;i<nel;i++)indices[i]=i;
     1053       
     1054        Vector<IssmDouble>* vloadcopy=new Vector<IssmDouble>(nel);
     1055        IssmDouble* loadcopy=xNew<IssmDouble>(nel);
     1056       
     1057        vloadcopy->SetValues(nel,indices,load,INS_VAL);
     1058        vloadcopy->Assemble();
     1059
     1060
     1061        if(subload){
     1062                for (int i=0;i<femmodel->elements->Size();i++){
     1063                        if (slgeom->issubelement[loadtype][i]){
     1064                                int se= slgeom->subelementmapping[loadtype][i];
     1065                                IssmDouble subloadi=subload[se];
     1066                                Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
     1067                                vloadcopy->SetValue(element->Sid(),subloadi,ADD_VAL);
     1068                        }
     1069                }
     1070        }
     1071        vloadcopy->Assemble();
     1072        loadcopy=vloadcopy->ToMPISerial();
     1073
     1074        return loadcopy;
     1075
     1076} /*}}}*/
     1077

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

    r26177 r26222  
    2424        /*parameters: */
    2525        IssmDouble finaltime,dt,yts;
    26         bool       isoceancoupling,iscontrol,isautodiff;
     26        bool       iscontrol,isautodiff;
    2727        int        timestepping;
    2828        int        output_frequency,checkpoint_frequency;
    29         int        amr_frequency,amr_restart;
     29        int        amr_frequency;
    3030        char     **requested_outputs = NULL;
    3131
    … …  
    4545        femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum);
    4646        femmodel->parameters->FindParam(&timestepping,TimesteppingTypeEnum);
    47         femmodel->parameters->FindParam(&isoceancoupling,TransientIsoceancouplingEnum);
    4847        femmodel->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
    4948        femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum);
    5049        femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
    5150
    52         #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
    53         if(amr_frequency){
    54                 femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum);
    55                 if(amr_restart) femmodel->ReMesh();
    56         }
    57         #endif
    58 
    59         #if defined(_HAVE_OCEAN_ )
    60         if(isoceancoupling) OceanExchangeDatax(femmodel,true);
    61         #endif
     51        /*call modules that are not dependent on time stepping:*/
     52        transient_precore(femmodel);
    6253
    6354        while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime.
    … …  
    253244        }
    254245}/*}}}*/
     246void transient_precore(FemModel* femmodel){/*{{{*/
     247
     248        bool       isoceancoupling,isslc;
     249        int        amr_frequency,amr_restart;
     250
     251        femmodel->parameters->FindParam(&isoceancoupling,TransientIsoceancouplingEnum);
     252        femmodel->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
     253        femmodel->parameters->FindParam(&isslc,TransientIsslcEnum);
     254
     255        #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
     256        if(amr_frequency){
     257                femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum);
     258                if(amr_restart) femmodel->ReMesh();
     259        }
     260        #endif
     261
     262        #if defined(_HAVE_OCEAN_ )
     263        if(isoceancoupling) OceanExchangeDatax(femmodel,true);
     264        #endif
     265
     266        if(isslc) sealevelchange_initialgeometry(femmodel);
     267}/*}}}*/
    255268
    256269#ifdef _HAVE_CODIPACK_

  • issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.cpp

    r26149 r26222  
    99#include "../../toolkits/toolkits.h"
    1010
    11 void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze, IssmDouble** pareae, Elements* elements,bool spherical) {
     11void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze, IssmDouble** pareae, Elements* elements,bool spherical) { /*{{{*/
    1212
    1313        /*figure out how many vertices we have: */
    … …  
    5353        else xDelete<IssmDouble>(areae);
    5454
    55 }
     55} /*}}}*/
     56void ElementCoordinatesx( IssmDouble** plonge, IssmDouble** plate, IssmDouble** pareae, Elements* elements) { /*{{{*/
     57
     58        /*figure out how many vertices we have: */
     59        int numberofelements=elements->NumberOfElements();
     60
     61        Vector<IssmDouble>* vlonge=new Vector<IssmDouble>(numberofelements);
     62        Vector<IssmDouble>* vlate=new Vector<IssmDouble>(numberofelements);
     63        Vector<IssmDouble>* vareae=new Vector<IssmDouble>(numberofelements);
     64
     65        /*march through our elements: */
     66        for(Object* & object : elements->objects){
     67                Element* element=(Element*)object;
     68                element->ElementCoordinates(vlonge,vlate,vareae);
     69        }
     70
     71        /*Assemble*/
     72        vlonge->Assemble();
     73        vlate->Assemble();
     74        vareae->Assemble();
     75
     76        /*serialize: */
     77        IssmDouble* longe=vlonge->ToMPISerial();
     78        IssmDouble* late=vlate->ToMPISerial();
     79        IssmDouble* areae=vareae->ToMPISerial();
     80
     81        /*Free ressources: */
     82        delete vlonge;
     83        delete vlate;
     84        delete vareae;
     85
     86        /*output: */
     87        if(plonge) *plonge=longe;
     88        else xDelete<IssmDouble>(longe);
     89        if(plate) *plate=late;
     90        else xDelete<IssmDouble>(late);
     91        if(pareae) *pareae=areae;
     92        else xDelete<IssmDouble>(areae);
     93
     94} /*}}}*/

  • issm/trunk-jpl/src/c/modules/ElementCoordinatesx/ElementCoordinatesx.h

    r26149 r26222  
    99/* local prototypes: */
    1010void ElementCoordinatesx( IssmDouble** pxe, IssmDouble** pye, IssmDouble** pze,IssmDouble** pareae, Elements* elements,bool spherical=false);
     11void ElementCoordinatesx( IssmDouble** plonge, IssmDouble** plate, IssmDouble** pareae, Elements* elements);
    1112
    1213#endif  /* _ELEMENT_COORDINATESX_H */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r26208 r26222  
    331331                        parameters->AddObject(iomodel->CopyConstantObject("md.materials.thermal_exchange_velocity",MaterialsThermalExchangeVelocityEnum));
    332332                        parameters->AddObject(iomodel->CopyConstantObject("md.constants.g",ConstantsGEnum));
     333                        parameters->AddObject(iomodel->CopyConstantObject("md.constants.gravitational_constant",ConstantsNewtonGravityEnum));
    333334                        parameters->AddObject(iomodel->CopyConstantObject("md.materials.rheology_law",MaterialsRheologyLawEnum));
    334335                        parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum));
    … …  
    373374                        }
    374375                        parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum));
     376                        parameters->AddObject(iomodel->CopyConstantObject("md.constants.gravitational_constant",ConstantsNewtonGravityEnum));
    375377                        /*Free rssources:*/
    376378                        xDelete<int>(nature);

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

    r26210 r26222  
    106106syn keyword cConstant ConfigurationTypeEnum
    107107syn keyword cConstant ConstantsGEnum
     108syn keyword cConstant ConstantsNewtonGravityEnum
    108109syn keyword cConstant ConstantsReferencetemperatureEnum
    109110syn keyword cConstant ConstantsYtsEnum
    … …  
    365366syn keyword cConstant SolidearthSettingsComputesealevelchangeEnum
    366367syn keyword cConstant SolidearthSettingsGRDEnum
    367 syn keyword cConstant SolidearthSettingsGlfractionEnum
    368368syn keyword cConstant SolidearthSettingsRunFrequencyEnum
    369369syn keyword cConstant SealevelchangeHElasticEnum
    … …  
    506506syn keyword cConstant TransientRequestedOutputsEnum
    507507syn keyword cConstant VelocityEnum
     508syn keyword cConstant XxeEnum
     509syn keyword cConstant YyeEnum
     510syn keyword cConstant ZzeEnum
     511syn keyword cConstant AreaeEnum
    508512syn keyword cConstant WorldCommEnum
    509513syn keyword cConstant ParametersENDEnum
    … …  
    745749syn keyword cConstant SealevelGRDEnum
    746750syn keyword cConstant SealevelBarystaticMaskEnum
     751syn keyword cConstant SealevelBarystaticIceMaskEnum
     752syn keyword cConstant SealevelBarystaticIceWeightsEnum
     753syn keyword cConstant SealevelBarystaticIceAreaEnum
     754syn keyword cConstant SealevelBarystaticIceLatbarEnum
     755syn keyword cConstant SealevelBarystaticIceLongbarEnum
     756syn keyword cConstant SealevelBarystaticIceLoadEnum
     757syn keyword cConstant SealevelBarystaticHydroMaskEnum
     758syn keyword cConstant SealevelBarystaticHydroWeightsEnum
     759syn keyword cConstant SealevelBarystaticHydroAreaEnum
     760syn keyword cConstant SealevelBarystaticHydroLatbarEnum
     761syn keyword cConstant SealevelBarystaticHydroLongbarEnum
     762syn keyword cConstant SealevelBarystaticHydroLoadEnum
     763syn keyword cConstant SealevelBarystaticBpMaskEnum
     764syn keyword cConstant SealevelBarystaticBpWeightsEnum
     765syn keyword cConstant SealevelBarystaticBpAreaEnum
     766syn keyword cConstant SealevelBarystaticBpLoadEnum
    747767syn keyword cConstant SealevelBarystaticOceanMaskEnum
     768syn keyword cConstant SealevelBarystaticOceanWeightsEnum
     769syn keyword cConstant SealevelBarystaticOceanAreaEnum
     770syn keyword cConstant SealevelBarystaticOceanLatbarEnum
     771syn keyword cConstant SealevelBarystaticOceanLongbarEnum
     772syn keyword cConstant SealevelBarystaticOceanLoadEnum
    748773syn keyword cConstant SealevelNEsaEnum
    749774syn keyword cConstant SealevelNEsaRateEnum
    … …  
    769794syn keyword cConstant SealevelchangeGEEnum
    770795syn keyword cConstant SealevelchangeGNEnum
     796syn keyword cConstant SealevelchangeGsubelOceanEnum
     797syn keyword cConstant SealevelchangeGUsubelOceanEnum
     798syn keyword cConstant SealevelchangeGEsubelOceanEnum
     799syn keyword cConstant SealevelchangeGNsubelOceanEnum
     800syn keyword cConstant SealevelchangeGsubelIceEnum
     801syn keyword cConstant SealevelchangeGUsubelIceEnum
     802syn keyword cConstant SealevelchangeGEsubelIceEnum
     803syn keyword cConstant SealevelchangeGNsubelIceEnum
     804syn keyword cConstant SealevelchangeGsubelHydroEnum
     805syn keyword cConstant SealevelchangeGUsubelHydroEnum
     806syn keyword cConstant SealevelchangeGEsubelHydroEnum
     807syn keyword cConstant SealevelchangeGNsubelHydroEnum
    771808syn keyword cConstant SedimentHeadEnum
    772809syn keyword cConstant SedimentHeadOldEnum
    … …  
    15401577syn keyword cType RiftStruct
    15411578syn keyword cType Riftfront
    1542 syn keyword cType SealevelMasks
     1579syn keyword cType SealevelGeometry
    15431580syn keyword cType Seg
    15441581syn keyword cType SegInput

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

    r26210 r26222  
    100100        ConfigurationTypeEnum,
    101101        ConstantsGEnum,
     102        ConstantsNewtonGravityEnum,
    102103        ConstantsReferencetemperatureEnum,
    103104        ConstantsYtsEnum,
    … …  
    359360        SolidearthSettingsComputesealevelchangeEnum,
    360361        SolidearthSettingsGRDEnum,
    361         SolidearthSettingsGlfractionEnum,
    362362        SolidearthSettingsRunFrequencyEnum,
    363363        SealevelchangeHElasticEnum,
    … …  
    500500        TransientRequestedOutputsEnum,
    501501        VelocityEnum,
     502        XxeEnum,
     503        YyeEnum,
     504        ZzeEnum,
     505        AreaeEnum,
    502506        WorldCommEnum,
    503507        /*}}}*/
    … …  
    741745        SealevelGRDEnum,
    742746        SealevelBarystaticMaskEnum,
     747        SealevelBarystaticIceMaskEnum,
     748        SealevelBarystaticIceWeightsEnum,
     749        SealevelBarystaticIceAreaEnum,
     750        SealevelBarystaticIceLatbarEnum,
     751        SealevelBarystaticIceLongbarEnum,
     752        SealevelBarystaticIceLoadEnum,
     753        SealevelBarystaticHydroMaskEnum,
     754        SealevelBarystaticHydroWeightsEnum,
     755        SealevelBarystaticHydroAreaEnum,
     756        SealevelBarystaticHydroLatbarEnum,
     757        SealevelBarystaticHydroLongbarEnum,
     758        SealevelBarystaticHydroLoadEnum,
     759        SealevelBarystaticBpMaskEnum,
     760        SealevelBarystaticBpWeightsEnum,
     761        SealevelBarystaticBpAreaEnum,
     762        SealevelBarystaticBpLoadEnum,
    743763        SealevelBarystaticOceanMaskEnum,
     764        SealevelBarystaticOceanWeightsEnum,
     765        SealevelBarystaticOceanAreaEnum,
     766        SealevelBarystaticOceanLatbarEnum,
     767        SealevelBarystaticOceanLongbarEnum,
     768        SealevelBarystaticOceanLoadEnum,
    744769        SealevelNEsaEnum,
    745770        SealevelNEsaRateEnum,
    … …  
    765790        SealevelchangeGEEnum,
    766791        SealevelchangeGNEnum,
     792        SealevelchangeGsubelOceanEnum,
     793        SealevelchangeGUsubelOceanEnum,
     794        SealevelchangeGEsubelOceanEnum,
     795        SealevelchangeGNsubelOceanEnum,
     796        SealevelchangeGsubelIceEnum,
     797        SealevelchangeGUsubelIceEnum,
     798        SealevelchangeGEsubelIceEnum,
     799        SealevelchangeGNsubelIceEnum,
     800        SealevelchangeGsubelHydroEnum,
     801        SealevelchangeGUsubelHydroEnum,
     802        SealevelchangeGEsubelHydroEnum,
     803        SealevelchangeGNsubelHydroEnum,
    767804        SedimentHeadEnum,
    768805        SedimentHeadOldEnum,

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

    r26210 r26222  
    108108                case ConfigurationTypeEnum : return "ConfigurationType";
    109109                case ConstantsGEnum : return "ConstantsG";
     110                case ConstantsNewtonGravityEnum : return "ConstantsNewtonGravity";
    110111                case ConstantsReferencetemperatureEnum : return "ConstantsReferencetemperature";
    111112                case ConstantsYtsEnum : return "ConstantsYts";
    … …  
    367368                case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange";
    368369                case SolidearthSettingsGRDEnum : return "SolidearthSettingsGRD";
    369                 case SolidearthSettingsGlfractionEnum : return "SolidearthSettingsGlfraction";
    370370                case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency";
    371371                case SealevelchangeHElasticEnum : return "SealevelchangeHElastic";
    … …  
    508508                case TransientRequestedOutputsEnum : return "TransientRequestedOutputs";
    509509                case VelocityEnum : return "Velocity";
     510                case XxeEnum : return "Xxe";
     511                case YyeEnum : return "Yye";
     512                case ZzeEnum : return "Zze";
     513                case AreaeEnum : return "Areae";
    510514                case WorldCommEnum : return "WorldComm";
    511515                case ParametersENDEnum : return "ParametersEND";
    … …  
    747751                case SealevelGRDEnum : return "SealevelGRD";
    748752                case SealevelBarystaticMaskEnum : return "SealevelBarystaticMask";
     753                case SealevelBarystaticIceMaskEnum : return "SealevelBarystaticIceMask";
     754                case SealevelBarystaticIceWeightsEnum : return "SealevelBarystaticIceWeights";
     755                case SealevelBarystaticIceAreaEnum : return "SealevelBarystaticIceArea";
     756                case SealevelBarystaticIceLatbarEnum : return "SealevelBarystaticIceLatbar";
     757                case SealevelBarystaticIceLongbarEnum : return "SealevelBarystaticIceLongbar";
     758                case SealevelBarystaticIceLoadEnum : return "SealevelBarystaticIceLoad";
     759                case SealevelBarystaticHydroMaskEnum : return "SealevelBarystaticHydroMask";
     760                case SealevelBarystaticHydroWeightsEnum : return "SealevelBarystaticHydroWeights";
     761                case SealevelBarystaticHydroAreaEnum : return "SealevelBarystaticHydroArea";
     762                case SealevelBarystaticHydroLatbarEnum : return "SealevelBarystaticHydroLatbar";
     763                case SealevelBarystaticHydroLongbarEnum : return "SealevelBarystaticHydroLongbar";
     764                case SealevelBarystaticHydroLoadEnum : return "SealevelBarystaticHydroLoad";
     765                case SealevelBarystaticBpMaskEnum : return "SealevelBarystaticBpMask";
     766                case SealevelBarystaticBpWeightsEnum : return "SealevelBarystaticBpWeights";
     767                case SealevelBarystaticBpAreaEnum : return "SealevelBarystaticBpArea";
     768                case SealevelBarystaticBpLoadEnum : return "SealevelBarystaticBpLoad";
    749769                case SealevelBarystaticOceanMaskEnum : return "SealevelBarystaticOceanMask";
     770                case SealevelBarystaticOceanWeightsEnum : return "SealevelBarystaticOceanWeights";
     771                case SealevelBarystaticOceanAreaEnum : return "SealevelBarystaticOceanArea";
     772                case SealevelBarystaticOceanLatbarEnum : return "SealevelBarystaticOceanLatbar";
     773                case SealevelBarystaticOceanLongbarEnum : return "SealevelBarystaticOceanLongbar";
     774                case SealevelBarystaticOceanLoadEnum : return "SealevelBarystaticOceanLoad";
    750775                case SealevelNEsaEnum : return "SealevelNEsa";
    751776                case SealevelNEsaRateEnum : return "SealevelNEsaRate";
    … …  
    771796                case SealevelchangeGEEnum : return "SealevelchangeGE";
    772797                case SealevelchangeGNEnum : return "SealevelchangeGN";
     798                case SealevelchangeGsubelOceanEnum : return "SealevelchangeGsubelOcean";
     799                case SealevelchangeGUsubelOceanEnum : return "SealevelchangeGUsubelOcean";
     800                case SealevelchangeGEsubelOceanEnum : return "SealevelchangeGEsubelOcean";
     801                case SealevelchangeGNsubelOceanEnum : return "SealevelchangeGNsubelOcean";
     802                case SealevelchangeGsubelIceEnum : return "SealevelchangeGsubelIce";
     803                case SealevelchangeGUsubelIceEnum : return "SealevelchangeGUsubelIce";
     804                case SealevelchangeGEsubelIceEnum : return "SealevelchangeGEsubelIce";
     805                case SealevelchangeGNsubelIceEnum : return "SealevelchangeGNsubelIce";
     806                case SealevelchangeGsubelHydroEnum : return "SealevelchangeGsubelHydro";
     807                case SealevelchangeGUsubelHydroEnum : return "SealevelchangeGUsubelHydro";
     808                case SealevelchangeGEsubelHydroEnum : return "SealevelchangeGEsubelHydro";
     809                case SealevelchangeGNsubelHydroEnum : return "SealevelchangeGNsubelHydro";
    773810                case SedimentHeadEnum : return "SedimentHead";
    774811                case SedimentHeadOldEnum : return "SedimentHeadOld";

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

    r26210 r26222  
    108108              else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
    109109              else if (strcmp(name,"ConstantsG")==0) return ConstantsGEnum;
     110              else if (strcmp(name,"ConstantsNewtonGravity")==0) return ConstantsNewtonGravityEnum;
    110111              else if (strcmp(name,"ConstantsReferencetemperature")==0) return ConstantsReferencetemperatureEnum;
    111112              else if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum;
    … …  
    136137              else if (strcmp(name,"DamageStabilization")==0) return DamageStabilizationEnum;
    137138              else if (strcmp(name,"DamageStressThreshold")==0) return DamageStressThresholdEnum;
    138               else if (strcmp(name,"DamageStressUBound")==0) return DamageStressUBoundEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
     142              if (strcmp(name,"DamageStressUBound")==0) return DamageStressUBoundEnum;
     143              else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
    143144              else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
    144145              else if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
    … …  
    259260              else if (strcmp(name,"LevelsetStabilization")==0) return LevelsetStabilizationEnum;
    260261              else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum;
    261               else if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"LoveForcingType")==0) return LoveForcingTypeEnum;
     265              if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum;
     266              else if (strcmp(name,"LoveForcingType")==0) return LoveForcingTypeEnum;
    266267              else if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum;
    267268              else if (strcmp(name,"LoveG0")==0) return LoveG0Enum;
    … …  
    373374              else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum;
    374375              else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum;
    375               else if (strcmp(name,"SolidearthSettingsGlfraction")==0) return SolidearthSettingsGlfractionEnum;
    376376              else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
    377377              else if (strcmp(name,"SealevelchangeHElastic")==0) return SealevelchangeHElasticEnum;
    … …  
    520520              else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
    521521              else if (strcmp(name,"Velocity")==0) return VelocityEnum;
     522              else if (strcmp(name,"Xxe")==0) return XxeEnum;
     523              else if (strcmp(name,"Yye")==0) return YyeEnum;
     524              else if (strcmp(name,"Zze")==0) return ZzeEnum;
     525              else if (strcmp(name,"Areae")==0) return AreaeEnum;
    522526              else if (strcmp(name,"WorldComm")==0) return WorldCommEnum;
    523527              else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
    … …  
    625629              else if (strcmp(name,"EffectivePressureTransient")==0) return EffectivePressureTransientEnum;
    626630              else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
    627               else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
    628635              else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
    629636              else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
    630637              else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
     638              else if (strcmp(name,"EplHeadSlopeY")==0) return EplHeadSlopeYEnum;
    635639              else if (strcmp(name,"EplHeadSubstep")==0) return EplHeadSubstepEnum;
    636640              else if (strcmp(name,"EplHeadTransient")==0) return EplHeadTransientEnum;
    … …  
    748752              else if (strcmp(name,"Radar")==0) return RadarEnum;
    749753              else if (strcmp(name,"RadarAttenuationMacGregor")==0) return RadarAttenuationMacGregorEnum;
    750               else if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum;
    751758              else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum;
    752759              else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum;
    753760              else if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
     761              else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
    758762              else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
    759763              else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum;
    … …  
    765769              else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
    766770              else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum;
     771              else if (strcmp(name,"SealevelBarystaticIceMask")==0) return SealevelBarystaticIceMaskEnum;
     772              else if (strcmp(name,"SealevelBarystaticIceWeights")==0) return SealevelBarystaticIceWeightsEnum;
     773              else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum;
     774              else if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum;
     775              else if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum;
     776              else if (strcmp(name,"SealevelBarystaticIceLoad")==0) return SealevelBarystaticIceLoadEnum;
     777              else if (strcmp(name,"SealevelBarystaticHydroMask")==0) return SealevelBarystaticHydroMaskEnum;
     778              else if (strcmp(name,"SealevelBarystaticHydroWeights")==0) return SealevelBarystaticHydroWeightsEnum;
     779              else if (strcmp(name,"SealevelBarystaticHydroArea")==0) return SealevelBarystaticHydroAreaEnum;
     780              else if (strcmp(name,"SealevelBarystaticHydroLatbar")==0) return SealevelBarystaticHydroLatbarEnum;
     781              else if (strcmp(name,"SealevelBarystaticHydroLongbar")==0) return SealevelBarystaticHydroLongbarEnum;
     782              else if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
     783              else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum;
     784              else if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum;
     785              else if (strcmp(name,"SealevelBarystaticBpArea")==0) return SealevelBarystaticBpAreaEnum;
     786              else if (strcmp(name,"SealevelBarystaticBpLoad")==0) return SealevelBarystaticBpLoadEnum;
    767787              else if (strcmp(name,"SealevelBarystaticOceanMask")==0) return SealevelBarystaticOceanMaskEnum;
     788              else if (strcmp(name,"SealevelBarystaticOceanWeights")==0) return SealevelBarystaticOceanWeightsEnum;
     789              else if (strcmp(name,"SealevelBarystaticOceanArea")==0) return SealevelBarystaticOceanAreaEnum;
     790              else if (strcmp(name,"SealevelBarystaticOceanLatbar")==0) return SealevelBarystaticOceanLatbarEnum;
     791              else if (strcmp(name,"SealevelBarystaticOceanLongbar")==0) return SealevelBarystaticOceanLongbarEnum;
     792              else if (strcmp(name,"SealevelBarystaticOceanLoad")==0) return SealevelBarystaticOceanLoadEnum;
    768793              else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum;
    769794              else if (strcmp(name,"SealevelNEsaRate")==0) return SealevelNEsaRateEnum;
    … …  
    789814              else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
    790815              else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum;
     816              else if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum;
     817              else if (strcmp(name,"SealevelchangeGUsubelOcean")==0) return SealevelchangeGUsubelOceanEnum;
     818              else if (strcmp(name,"SealevelchangeGEsubelOcean")==0) return SealevelchangeGEsubelOceanEnum;
     819              else if (strcmp(name,"SealevelchangeGNsubelOcean")==0) return SealevelchangeGNsubelOceanEnum;
     820              else if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum;
     821              else if (strcmp(name,"SealevelchangeGUsubelIce")==0) return SealevelchangeGUsubelIceEnum;
     822              else if (strcmp(name,"SealevelchangeGEsubelIce")==0) return SealevelchangeGEsubelIceEnum;
     823              else if (strcmp(name,"SealevelchangeGNsubelIce")==0) return SealevelchangeGNsubelIceEnum;
     824              else if (strcmp(name,"SealevelchangeGsubelHydro")==0) return SealevelchangeGsubelHydroEnum;
     825              else if (strcmp(name,"SealevelchangeGUsubelHydro")==0) return SealevelchangeGUsubelHydroEnum;
     826              else if (strcmp(name,"SealevelchangeGEsubelHydro")==0) return SealevelchangeGEsubelHydroEnum;
     827              else if (strcmp(name,"SealevelchangeGNsubelHydro")==0) return SealevelchangeGNsubelHydroEnum;
    791828              else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    792829              else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    … …  
    838875              else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
    839876              else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
    840               else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    841881              else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
    842882              else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
    … …  
    875915              else if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum;
    876916              else if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum;
     917              else if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum;
    881918              else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum;
    882919              else if (strcmp(name,"SmbS0t")==0) return SmbS0tEnum;
    … …  
    961998              else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
    962999              else if (strcmp(name,"VxObs")==0) return VxObsEnum;
    963               else if (strcmp(name,"VxShear")==0) return VxShearEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"VxShear")==0) return VxShearEnum;
    9641004              else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum;
    9651005              else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    … …  
    9981038              else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
    9991039              else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
     1040              else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
    10041041              else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
    10051042              else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
    … …  
    10841121              else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum;
    10851122              else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
    1086               else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
    10871127              else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
    10881128              else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum;
    … …  
    11211161              else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    11221162              else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
     1163              else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
    11271164              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    11281165              else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
    … …  
    12071244              else if (strcmp(name,"Gradient1")==0) return Gradient1Enum;
    12081245              else if (strcmp(name,"Gradient2")==0) return Gradient2Enum;
    1209               else if (strcmp(name,"Gradient3")==0) return Gradient3Enum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"Gradient3")==0) return Gradient3Enum;
    12101250              else if (strcmp(name,"Gradient4")==0) return Gradient4Enum;
    12111251              else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum;
    … …  
    12441284              else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
    12451285              else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"IntParam")==0) return IntParamEnum;
     1286              else if (strcmp(name,"IntParam")==0) return IntParamEnum;
    12501287              else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
    12511288              else if (strcmp(name,"Inputs")==0) return InputsEnum;
    … …  
    13301367              else if (strcmp(name,"OceantransportSolution")==0) return OceantransportSolutionEnum;
    13311368              else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
    1332               else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
    13331373              else if (strcmp(name,"Open")==0) return OpenEnum;
    13341374              else if (strcmp(name,"Option")==0) return OptionEnum;
    … …  
    13671407              else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
    13681408              else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
     1409              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    13731410              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    13741411              else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
    … …  
    14531490              else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
    14541491              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
    1455               else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
    14561496              else if (strcmp(name,"Water")==0) return WaterEnum;
    14571497              else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
    … …  
    14761516              else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
    14771517              else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
    1478          else stage=13;
     1518         else stage=14;
    14791519   }
    14801520        /*If we reach this point, the string provided has not been found*/

  • issm/trunk-jpl/src/m/classes/constants.m

    r20896 r26222  
    1010                yts                  = 0.;
    1111                referencetemperature = 0.;
     12                gravitational_constant = 0.;
    1213        end
    1314        methods
    … …  
    3334                        %the reference temperature for enthalpy model (cf Aschwanden)
    3435                        self.referencetemperature=223.15;
     36               
     37                        %gravity:
     38                        self.gravitational_constant = 6.67259e-11;
    3539
    3640                end % }}}
    … …  
    4145                        md = checkfield(md,'fieldname','constants.yts','>',0,'size',[1 1]);
    4246                        md = checkfield(md,'fieldname','constants.referencetemperature','size',[1 1]);
     47                        md = checkfield(md,'fieldname','constants.gravitational_constant','size',[1 1]);
    4348
    4449                end % }}}
    … …  
    5055                        fielddisplay(self,'yts','number of seconds in a year [s/yr]');
    5156                        fielddisplay(self,'referencetemperature','reference temperature used in the enthalpy model [K]');
     57                        fielddisplay(self,'gravitational_constant','Newtonian constant of gravitation [m^3/kg/s^2]');
    5258
    5359                end % }}}
    … …  
    5662                        WriteData(fid,prefix,'object',self,'fieldname','yts','format','Double');
    5763                        WriteData(fid,prefix,'object',self,'fieldname','referencetemperature','format','Double');
     64                        WriteData(fid,prefix,'object',self,'fieldname','gravitational_constant','format','Double');
    5865                end % }}}
    5966                function savemodeljs(self,fid,modelname) % {{{
    … …  
    6370                        writejsdouble(fid,[modelname '.constants.yts'],self.yts);
    6471                        writejsdouble(fid,[modelname '.constants.referencetemperature'],self.referencetemperature);
     72                        writejsdouble(fid,[modelname '.constants.gravitational_constant'],self.gravitational_constant);
    6573
    6674                end % }}}

  • issm/trunk-jpl/src/m/classes/solidearthsettings.m

    r26165 r26222  
    1919                degacc                 = 0; %degree increment for resolution of Green tables
    2020                horiz                  = 0; %compute horizontal deformation
    21                 glfraction             = 1; %barystatic contribution full or fractional (default fractional)
    2221                grdmodel               = 0; %grd model (0 by default, 1 for elastic, 2 for Ivins)
    2322                cross_section_shape    = 0; %cross section only used when grd model is Ivins
    … …  
    5655                self.runfrequency=1;
    5756               
    58                 %fractional contribution:
    59                 self.glfraction=1;
    60        
    6157                %horizontal displacement?  (not by default)
    6258                self.horiz=0;
    … …  
    8076                        md = checkfield(md,'fieldname','solidearth.settings.degacc','size',[1 1],'>=',1e-10);
    8177                        md = checkfield(md,'fieldname','solidearth.settings.horiz','NaN',1,'Inf',1,'values',[0 1]);
    82                         md = checkfield(md,'fieldname','solidearth.settings.glfraction','values',[0 1]);
    8378                        md = checkfield(md,'fieldname','solidearth.settings.grdmodel','values',[1 2]);
    8479                        md = checkfield(md,'fieldname','solidearth.settings.cross_section_shape','numel',[1],'values',[1,2]);
    … …  
    123118                        fielddisplay(self,'rotation','earth rotational potential perturbation');
    124119                        fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
    125                         fielddisplay(self,'glfraction','contribute fractionally (default, 1) to barystatic sea level');
    126120                        fielddisplay(self,'grdmodel','type of deformation model, 1 for elastic, 2 for visco-elastic from Ivins');
    127121                        fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
    … …  
    141135                        WriteData(fid,prefix,'object',self,'fieldname','isgrd','name','md.solidearth.settings.isgrd','format','Integer');
    142136                        WriteData(fid,prefix,'object',self,'fieldname','compute_bp_grd','name','md.solidearth.settings.compute_bp_grd','format','Integer');
    143                         WriteData(fid,prefix,'object',self,'fieldname','glfraction','name','md.solidearth.settings.glfraction','format','Integer');
    144137                        WriteData(fid,prefix,'object',self,'fieldname','grdmodel','name','md.solidearth.settings.grdmodel','format','Integer');
    145138                        WriteData(fid,prefix,'object',self,'fieldname','cross_section_shape','name','md.solidearth.settings.cross_section_shape','format','Integer');
    … …  
    156149                        writejsdouble(fid,[modelname '.solidearth.settings.run_frequency'],self.run_frequency);
    157150                        writejsdouble(fid,[modelname '.solidearth.settings.degacc'],self.degacc);
    158                         writejsdouble(fid,[modelname '.solidearth.settings.glfraction'],self.glfraction);
    159151                        writejsdouble(fid,[modelname '.solidearth.settings.cross_section_shape'],self.cross_section_shape);
    160152                end % }}}

  • issm/trunk-jpl/test/NightlyRun/test2002.m

    r26165 r26222  
    2020md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
    2121%antarctica
    22 late=sum(md.mesh.lat(md.mesh.elements),2)/3;
    23 longe=sum(md.mesh.long(md.mesh.elements),2)/3;
     22xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
     23ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3;
     24ze=md.mesh.z(md.mesh.elements)*[1;1;1]/3;
     25re=sqrt(xe.^2+ye.^2+ze.^2);
     26
     27late=asind(ze./re);
     28longe=atan2d(ye,xe);
    2429pos=find(late < -80);
    2530md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
    2631posant=pos;
    2732%greenland
    28 pos=find(late>70 & late<80 & longe>-60 & longe<-30);
     33pos=find(late>60 & late<90 & longe>-75 & longe<-15);
    2934md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
    3035posgre=pos;
    … …  
    5156md.timestepping.time_step=1;
    5257md.timestepping.final_time=1;
    53 
    5458
    5559md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
    … …  
    9397md.solidearth.settings.elastic=0;
    9498md.solidearth.settings.rotation=0;
     99
    95100md=solve(md,'Transient');
    96101Seustatic=md.results.TransientSolution.Sealevel;
Note: See TracChangeset for help on using the changeset viewer.