Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24939) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24940) @@ -321,4 +321,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.slr.eartharea",SealevelEarthAreaEnum)); /*Deal with dsl multi-model ensembles: {{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24939) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24940) @@ -371,5 +371,5 @@ #ifdef _HAVE_ESA_ virtual void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast, Vector* pX, Vector* pY,IssmDouble* xx,IssmDouble* yy)=0; - virtual void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0; + virtual void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz)=0; #endif #ifdef _HAVE_SEALEVELRISE_ @@ -378,9 +378,9 @@ virtual IssmDouble GetAreaSpherical(void)=0; virtual IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks)=0; - virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea)=0; - virtual void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; + virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0; + virtual void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea)=0; virtual void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0; - virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0; - virtual void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0; + virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0; + virtual void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24939) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24940) @@ -209,14 +209,14 @@ #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; - void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; #endif #ifdef _HAVE_SEALEVELRISE_ IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; + void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24939) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24940) @@ -169,13 +169,13 @@ #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; - void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; + void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24939) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24940) @@ -175,13 +175,13 @@ #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; - void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; #endif #ifdef _HAVE_SEALEVELRISE_ void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; + void SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24939) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24940) @@ -5289,5 +5289,5 @@ } /*}}}*/ -void Tria::EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/ +void Tria::EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){ /*{{{*/ /*diverse:*/ @@ -5296,5 +5296,5 @@ IssmDouble llr_list[NUMVERTICES][3]; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble I; //ice/water loading IssmDouble late,longe,re; @@ -5328,4 +5328,7 @@ rho_ice=FindParam(MaterialsRhoIceEnum); rho_earth=FindParam(MaterialsEarthDensityEnum); + + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); /*how many dofs are we working with here? */ @@ -5465,5 +5468,5 @@ } /*}}}*/ -void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks, IssmDouble eartharea){/*{{{*/ +void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/ /*early return if we are not on an ice cap OR ocean:*/ if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){ @@ -5475,6 +5478,9 @@ /*Compute area of element:*/ - IssmDouble area; + IssmDouble area,eartharea; area=GetAreaSpherical(); + + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); /*Compute lat,long,radius of elemental centroid: */ @@ -5661,5 +5667,5 @@ } /*}}}*/ -void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){ /*{{{*/ /*Computational flags:*/ @@ -5672,14 +5678,14 @@ /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested *bottom pressure fingerprints:*/ - if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); + if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea); } //if(!IsIceInElement()){ /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/ - this->SealevelriseEustaticIce(Sgi,peustatic,masks, latitude,longitude,radius,oceanarea,eartharea); + this->SealevelriseEustaticIce(Sgi,peustatic,masks, latitude,longitude,radius,oceanarea); //} } /*}}}*/ -void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){ /*{{{*/ /*diverse:*/ @@ -5687,5 +5693,5 @@ bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) IssmDouble rho; @@ -5744,4 +5750,7 @@ this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum); + + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); /*recover precomputed green function kernels:*/ @@ -5836,5 +5845,5 @@ } /*}}}*/ -void Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){ /*{{{*/ /*diverse:*/ @@ -5842,5 +5851,5 @@ bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) IssmDouble rho; @@ -5875,4 +5884,7 @@ rho_water=FindParam(MaterialsRhoSeawaterEnum); rho_earth=FindParam(MaterialsEarthDensityEnum); + + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); /*recover love numbers and computational flags: */ @@ -5982,5 +5994,5 @@ } /*}}}*/ -void Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){ /*{{{*/ /*diverse:*/ @@ -5988,5 +6000,5 @@ bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble S; //change in water water level(Farrel and Clarke, Equ. 4) IssmDouble late,longe; @@ -6026,4 +6038,7 @@ this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + /*early return if rigid or elastic not requested:*/ if(!computerigid && !computeelastic) return; @@ -6081,5 +6096,5 @@ } /*}}}*/ -void Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/ +void Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){ /*{{{*/ /*diverse:*/ @@ -6088,5 +6103,5 @@ IssmDouble llr_list[NUMVERTICES][3]; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble I, S; //change in relative ice thickness and sea level IssmDouble late,longe,re; @@ -6140,4 +6155,7 @@ /*how many dofs are we working with here? */ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); /*retrieve indices:*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24939) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24940) @@ -160,16 +160,16 @@ #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast, Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy); - void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea); + void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz); #endif #ifdef _HAVE_SEALEVELRISE_ IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks); - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble eartharea); + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks); void SetSealevelMasks(SealevelMasks* masks); void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius); - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea); - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz); + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea); + void SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea); + void SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea); + void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius); + void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24939) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24940) @@ -4635,7 +4635,4 @@ void FemModel::EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ - IssmDouble eartharea=0; - IssmDouble eartharea_cpu=0; - int ns,nsmax; @@ -4646,8 +4643,5 @@ for(int i=0;i(elements->GetObjectByOffset(i)); - eartharea_cpu += element->GetAreaSpherical(); - } - ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + } /*Figure out max of ns: */ @@ -4659,5 +4653,5 @@ if(i(elements->GetObjectByOffset(i)); - element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea); + element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz); } if(i%100==0){ @@ -4694,5 +4688,5 @@ } /*}}}*/ -void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ +void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ /*serialized vectors:*/ @@ -4702,6 +4696,4 @@ IssmDouble oceanarea = 0.; IssmDouble oceanarea_cpu = 0.; - IssmDouble eartharea = 0.; - IssmDouble eartharea_cpu = 0.; /*Initialize temporary vector that will be used to sum eustatic components @@ -4716,5 +4708,4 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); IssmDouble area=element->GetAreaSpherical(); - eartharea_cpu += area; if (masks->isoceanin[i]) oceanarea_cpu += area; } @@ -4723,11 +4714,8 @@ _assert_(oceanarea>0.); - ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - /*Call the sea level rise core: */ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, latitude,longitude,radius,oceanarea,eartharea); + element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, latitude,longitude,radius,oceanarea); eustatic_cpu+=eustatic_cpu_e; } @@ -4747,5 +4735,4 @@ /*Assign output pointers:*/ - *peartharea = eartharea; *poceanarea = oceanarea; *peustatic = eustatic; @@ -4753,5 +4740,5 @@ } /*}}}*/ -void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ +void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ /*serialized vectors:*/ @@ -4774,5 +4761,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius,eartharea); + element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius); } pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL); @@ -4786,5 +4773,5 @@ } /*}}}*/ -void FemModel::SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ +void FemModel::SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ /*serialized vectors:*/ @@ -4802,5 +4789,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks, eartharea); + element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks ); moi_list_cpu[0] += moi_list[0]; moi_list_cpu[1] += moi_list[1]; @@ -4862,5 +4849,5 @@ } /*}}}*/ -void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ +void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ /*serialized vectors:*/ @@ -4888,5 +4875,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,eartharea,horiz); + element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,horiz); } Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24939) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24940) @@ -163,9 +163,9 @@ #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseEustatic(Vector* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop); + void SealevelriseEustatic(Vector* pSgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop); void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); - void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop); - void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea,SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); - void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble eartharea,SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); + void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop); + void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); + void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); IssmDouble SealevelriseOceanAverage(Vector* Sg,SealevelMasks* masks, IssmDouble oceanarea); #endif Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 24939) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 24940) @@ -58,7 +58,7 @@ void sealevelrise_core_geometry(FemModel* femmodel); SealevelMasks* sealevelrise_core_masks(FemModel* femmodel); -Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* peartharea,IssmDouble* poceanarea); -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea); -void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* RSLg,IssmDouble eartharea, SealevelMasks* masks); +Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea); +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble oceanarea); +void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks); void sealevelrise_core_viscous(Vector** pU_gia,Vector** pN_gia,FemModel* femmodel,Vector* RSLg); void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg); Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24939) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24940) @@ -96,5 +96,5 @@ int geodetic=0; IssmDouble dt; - IssmDouble eartharea,oceanarea; + IssmDouble oceanarea; /*Should we even be here?:*/ @@ -153,11 +153,11 @@ /*call eustatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */ - RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&eartharea,&oceanarea); + RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&oceanarea); /*call non-eustatic core (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */ - RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,eartharea,oceanarea); + RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,oceanarea); /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */ - sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,eartharea,masks); + sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,masks); /*compute viscosus (GIA) geodetic signatures:*/ @@ -354,5 +354,5 @@ }/*}}}*/ -Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/ +Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/ /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ @@ -368,5 +368,5 @@ IssmDouble *radius = NULL; int loop; - IssmDouble eartharea,oceanarea; + IssmDouble oceanarea; /*outputs:*/ @@ -388,5 +388,5 @@ /*call the eustatic main module: */ - femmodel->SealevelriseEustatic(RSLgi,&eartharea, &oceanarea,&eustatic, masks, latitude, longitude, radius,loop); //this computes + femmodel->SealevelriseEustatic(RSLgi,&oceanarea,&eustatic, masks, latitude, longitude, radius,loop); //this computes /*we need to average RSLgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */ @@ -406,9 +406,8 @@ /*Assign output pointers and return: */ - *peartharea=eartharea; *poceanarea=oceanarea; return RSLgi; }/*}}}*/ -Vector* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea){ /*{{{*/ +Vector* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble oceanarea){ /*{{{*/ /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. @@ -479,5 +478,5 @@ /*call the non eustatic module: */ - femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, eartharea, masks, latitude,longitude, radius,verboseconvolution,loop); + femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, masks, latitude,longitude, radius,verboseconvolution,loop); /*assemble solution vector: */ @@ -488,5 +487,5 @@ /*call rotational feedback module: */ RSLgo_rot = new Vector(gsize); RSLgo_rot->Assemble(); - femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, eartharea, masks, latitude,longitude,radius); + femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks, latitude,longitude,radius); RSLgo_rot->Assemble(); @@ -537,5 +536,5 @@ return RSLg; } /*}}}*/ -void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg,IssmDouble eartharea, SealevelMasks* masks){ /*{{{*/ +void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks){ /*{{{*/ Vector *U_esa = NULL; @@ -577,5 +576,5 @@ /*call the elastic main modlule:*/ - femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,eartharea, masks, latitude,longitude,radius,xx,yy,zz,loop,horiz); + femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg, masks, latitude,longitude,radius,xx,yy,zz,loop,horiz); /*Assign output pointers:*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24939) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24940) @@ -302,4 +302,5 @@ syn keyword cConstant RootPathEnum syn keyword cConstant SaveResultsEnum +syn keyword cConstant SealevelEarthAreaEnum syn keyword cConstant SealevelEustaticEnum syn keyword cConstant SealevelriseAbstolEnum @@ -1431,4 +1432,5 @@ syn keyword cType RiftStruct syn keyword cType Riftfront +syn keyword cType SealevelMasks syn keyword cType Seg syn keyword cType SegInput2 Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24939) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24940) @@ -296,4 +296,5 @@ RootPathEnum, SaveResultsEnum, + SealevelEarthAreaEnum, SealevelEustaticEnum, SealevelriseAbstolEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24939) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24940) @@ -304,4 +304,5 @@ case RootPathEnum : return "RootPath"; case SaveResultsEnum : return "SaveResults"; + case SealevelEarthAreaEnum : return "SealevelEarthArea"; case SealevelEustaticEnum : return "SealevelEustatic"; case SealevelriseAbstolEnum : return "SealevelriseAbstol"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24939) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24940) @@ -310,4 +310,5 @@ else if (strcmp(name,"RootPath")==0) return RootPathEnum; else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum; + else if (strcmp(name,"SealevelEarthArea")==0) return SealevelEarthAreaEnum; else if (strcmp(name,"SealevelEustatic")==0) return SealevelEustaticEnum; else if (strcmp(name,"SealevelriseAbstol")==0) return SealevelriseAbstolEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum; else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum; - else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum; + if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum; + else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum; else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum; else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum; - else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; + if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum; + else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum; else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum; else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"Ice")==0) return IceEnum; else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; - else if (strcmp(name,"Input")==0) return InputEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; + if (strcmp(name,"Input")==0) return InputEnum; + else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum; else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum; else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; else if (strcmp(name,"SmbMassBalanceClimate")==0) return SmbMassBalanceClimateEnum; - else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; + if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; + else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum; else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Outputdefinition11")==0) return Outputdefinition11Enum; else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum; - else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum; + if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum; + else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum; else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum; else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; @@ -997,9 +998,9 @@ else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum; else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum; - else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; + if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum; + else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"Hydrologyshreve")==0) return HydrologyshreveEnum; else if (strcmp(name,"IceMass")==0) return IceMassEnum; - else if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; + if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum; + else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum; else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"Pengrid")==0) return PengridEnum; else if (strcmp(name,"Penpair")==0) return PenpairEnum; - else if (strcmp(name,"Penta")==0) return PentaEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"PentaInput")==0) return PentaInputEnum; + if (strcmp(name,"Penta")==0) return PentaEnum; + else if (strcmp(name,"PentaInput")==0) return PentaInputEnum; else if (strcmp(name,"Profiler")==0) return ProfilerEnum; else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum; Index: /issm/trunk-jpl/src/m/classes/slr.m =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.m (revision 24939) +++ /issm/trunk-jpl/src/m/classes/slr.m (revision 24940) @@ -33,4 +33,5 @@ Ngia = NaN; Ugia = NaN; + eartharea = 4*pi*planetradius('earth')^2; requested_outputs = {}; @@ -102,5 +103,6 @@ end - md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); + %md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'fieldname','slr.deltathickness','timeseries',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1); @@ -132,10 +134,10 @@ %cross check that whereever we have an ice load, the mask is <0 on each vertex: - pos=find(self.deltathickness); - maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:)); - [els,vertices]=find(maskpos>0); - if length(els), - warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!'); - end + %pos=find(self.deltathickness); + %maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:)); + %[els,vertices]=find(maskpos>0); + %if length(els), + % warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!'); + %end %check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not, @@ -221,4 +223,5 @@ WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','eartharea','format','Double'); %process requested outputs