Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24937) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24938) @@ -30,4 +30,5 @@ class DatasetInput2; class IoModel; +class SealevelMasks; class Gauss; class ElementVector; @@ -373,13 +374,13 @@ #endif #ifdef _HAVE_SEALEVELRISE_ + virtual void SetSealevelMasks(SealevelMasks* masks)=0; virtual IssmDouble GetArea3D(void)=0; virtual IssmDouble GetAreaSpherical(void)=0; - virtual IssmDouble OceanAverage(IssmDouble* Sg)=0; - virtual IssmDouble OceanArea(void)=0; - virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea)=0; - virtual void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=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 SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0; - virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0; - virtual void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,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,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; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24937) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24938) @@ -212,11 +212,11 @@ #endif #ifdef _HAVE_SEALEVELRISE_ - IssmDouble OceanArea(void){_error_("not implemented yet!");}; - IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");}; + 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 SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,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,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!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24937) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24938) @@ -172,11 +172,11 @@ #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_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,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; - IssmDouble OceanArea(void){_error_("not implemented yet!");}; - IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; + 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!");}; + 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 24937) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24938) @@ -178,11 +178,11 @@ #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_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 SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");}; - IssmDouble OceanArea(void){_error_("not implemented yet!");}; - IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; + 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!");}; + 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 24937) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24938) @@ -3016,5 +3016,4 @@ if(!IsIceInElement())return 0.; - //if(!IsIceOnlyInElement()) return 0; int domaintype; @@ -5447,14 +5446,7 @@ #endif #ifdef _HAVE_SEALEVELRISE_ -IssmDouble Tria::OceanArea(void){ /*{{{*/ - - if(IsOceanInElement()) return GetAreaSpherical(); - else return 0; - -} -/*}}}*/ -IssmDouble Tria::OceanAverage(IssmDouble* Sg){ /*{{{*/ - - if(IsOceanInElement()){ +IssmDouble Tria::OceanAverage(IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/ + + if(masks->isoceanin[this->lid]){ IssmDouble area; @@ -5473,7 +5465,7 @@ } /*}}}*/ -void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){/*{{{*/ +void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks, IssmDouble eartharea){/*{{{*/ /*early return if we are not on an ice cap OR ocean:*/ - if(!IsIceOnlyInElement() && !IsOceanInElement()){ + if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){ dI_list[0] = 0.0; // this is important!!! dI_list[1] = 0.0; // this is important!!! @@ -5528,5 +5520,5 @@ re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0; - if(IsOceanInElement()){ + if(masks->isoceanin[this->lid]){ IssmDouble rho_water, S; @@ -5546,5 +5538,5 @@ dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea; } - else if(IsIceOnlyInElement()){ + else if(masks->isiceonly[this->lid]){ IssmDouble rho_ice, I; @@ -5564,4 +5556,20 @@ return; }/*}}}*/ +void Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/ + + masks->isiceonly[this->lid]=this->IsIceOnlyInElement(); + masks->isoceanin[this->lid]=this->IsOceanInElement(); + + /*are we fully floating:*/ + Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input); + if (gr_input->GetInputMax()<=0)masks->isfullyfloating[this->lid]=true; + else masks->isfullyfloating[this->lid]=false; + + /*are we not fully grounded: */ + if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true; + else masks->notfullygrounded[this->lid]=false; + +} +/*}}}*/ void Tria::SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){ /*{{{*/ @@ -5653,5 +5661,5 @@ } /*}}}*/ -void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,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,IssmDouble eartharea){ /*{{{*/ /*Computational flags:*/ @@ -5661,5 +5669,5 @@ this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum); - if(!IsOceanInElement()){ + if(!masks->isoceanin[this->lid]){ /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested *bottom pressure fingerprints:*/ @@ -5668,10 +5676,10 @@ //if(!IsIceInElement()){ /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/ - this->SealevelriseEustaticIce(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); + this->SealevelriseEustaticIce(Sgi,peustatic,masks, latitude,longitude,radius,oceanarea,eartharea); //} } /*}}}*/ -void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,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,IssmDouble eartharea){ /*{{{*/ /*diverse:*/ @@ -5708,5 +5716,5 @@ /*early return if we are not on an ice cap:*/ - if(!IsIceOnlyInElement()){ + if(!masks->isiceonly[this->lid]){ constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! @@ -5715,6 +5723,5 @@ /*early return if we are fully floating:*/ - Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input); - if (gr_input->GetInputMax()<=0){ + if (masks->isfullyfloating[this->lid]){ constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! @@ -5723,7 +5730,5 @@ /*If we are here, we are on ice that is fully grounded or half-way to floating:*/ - if ((gr_input->GetInputMin())<0){ - notfullygrounded=true; //used later on. - } + if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. /*Inform mask: */ @@ -5831,5 +5836,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,IssmDouble eartharea){ /*{{{*/ /*diverse:*/ @@ -5977,5 +5982,5 @@ } /*}}}*/ -void Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/ /*diverse:*/ @@ -6011,5 +6016,5 @@ /*early return if we are not on the ocean:*/ - if (!IsOceanInElement()){ + if (!masks->isoceanin[this->lid]){ constant=0; this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum); return; @@ -6076,5 +6081,5 @@ } /*}}}*/ -void Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg,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,IssmDouble eartharea,int horiz){ /*{{{*/ /*diverse:*/ @@ -6116,9 +6121,8 @@ /*early return if we are not on the ocean or on an ice cap:*/ - if(!IsIceOnlyInElement() && !IsOceanInElement()) return; + if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]) return; /*early return if we are fully floating: */ - Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input); - if(gr_input->GetInputMax()<=0)return; + if(masks->isfullyfloating[this->lid])return; /*recover computational flags: */ @@ -6185,6 +6189,6 @@ /*we are going to use the result of these masks a lot: */ - bool isiceonlyinelement=IsIceOnlyInElement(); - bool isoceaninelement=IsOceanInElement(); + bool isiceonlyinelement=masks->isiceonly[this->lid]; + bool isoceaninelement=masks->isoceanin[this->lid]; for(int i=0;i* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ +void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ /*serialized vectors:*/ @@ -4717,5 +4717,5 @@ IssmDouble area=element->GetAreaSpherical(); eartharea_cpu += area; - if(element->IsOceanInElement()) oceanarea_cpu += area; + if (masks->isoceanin[i]) oceanarea_cpu += area; } ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); @@ -4729,5 +4729,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); + element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, latitude,longitude,radius,oceanarea,eartharea); eustatic_cpu+=eustatic_cpu_e; } @@ -4753,5 +4753,5 @@ } /*}}}*/ -void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ +void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ /*serialized vectors:*/ @@ -4774,5 +4774,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseNonEustatic(RSLgo,RSLg_old,latitude,longitude,radius,eartharea); + element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius,eartharea); } pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL); @@ -4786,5 +4786,5 @@ } /*}}}*/ -void FemModel::SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ +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){/*{{{*/ /*serialized vectors:*/ @@ -4802,5 +4802,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,eartharea); + element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks, eartharea); moi_list_cpu[0] += moi_list[0]; moi_list_cpu[1] += moi_list[1]; @@ -4862,5 +4862,5 @@ } /*}}}*/ -void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, IssmDouble eartharea, 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, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ /*serialized vectors:*/ @@ -4888,5 +4888,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseGeodetic(Up,North,East,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz); + element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,eartharea,horiz); } @@ -4908,5 +4908,5 @@ } /*}}}*/ -IssmDouble FemModel::SealevelriseOceanAverage(Vector* RSLg,IssmDouble oceanarea) { /*{{{*/ +IssmDouble FemModel::SealevelriseOceanAverage(Vector* RSLg,SealevelMasks* masks, IssmDouble oceanarea) { /*{{{*/ IssmDouble* RSLg_serial=NULL; @@ -4922,5 +4922,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - oceanvalue_cpu += element->OceanAverage(RSLg_serial); + oceanvalue_cpu += element->OceanAverage(RSLg_serial,masks); } Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24937) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24938) @@ -18,4 +18,5 @@ class Loads; class Materials; +class SealevelMasks; class Profiler; class Elements; @@ -162,10 +163,10 @@ #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseEustatic(Vector* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop); + void SealevelriseEustatic(Vector* pSgi, IssmDouble* peartharea, 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, 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,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); - void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); - IssmDouble SealevelriseOceanAverage(Vector* Sg,IssmDouble oceanarea); + 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); + IssmDouble SealevelriseOceanAverage(Vector* Sg,SealevelMasks* masks, IssmDouble oceanarea); #endif void HydrologyEPLupdateDomainx(IssmDouble* pEplcount); Index: /issm/trunk-jpl/src/c/classes/SealevelMasks.h =================================================================== --- /issm/trunk-jpl/src/c/classes/SealevelMasks.h (revision 24938) +++ /issm/trunk-jpl/src/c/classes/SealevelMasks.h (revision 24938) @@ -0,0 +1,34 @@ +/*!\file SealevelMasks.h + * \brief: header file for SealevelMask object + */ + +#ifndef _SEALEVELMASKS_H_ +#define _SEALEVELMASKS_H_ + +/*Headers:*/ + +class SealevelMasks { + + public: + + bool* isiceonly; + bool* isfullyfloating; + bool* notfullygrounded; + bool* isoceanin; + + /*SealevelMasks constructors, destructors :*/ + SealevelMasks(int nel){ + /*allocate fields:*/ + this->isiceonly=xNew(nel); + this->isfullyfloating=xNew(nel); + this->notfullygrounded=xNew(nel); + this->isoceanin=xNew(nel); + }; + ~SealevelMasks(){ + xDelete(this->isiceonly); + xDelete(this->isfullyfloating); + xDelete(this->notfullygrounded); + xDelete(this->isoceanin); + }; +}; +#endif /* _SEALEVELMASKS_H_ */ Index: /issm/trunk-jpl/src/c/classes/classes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/classes.h (revision 24937) +++ /issm/trunk-jpl/src/c/classes/classes.h (revision 24938) @@ -18,4 +18,5 @@ #include "./Massfluxatgate.h" #include "./Misfit.h" +#include "./SealevelMasks.h" #include "./Nodalvalue.h" #include "./Numberedcostfunction.h" Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 24937) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 24938) @@ -9,4 +9,5 @@ class FemModel; class Parameters; +class SealevelMasks; template class Matrix; template class Vector; @@ -56,7 +57,8 @@ void steric_core(FemModel* femmodel); void sealevelrise_core_geometry(FemModel* femmodel); -Vector* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea); -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,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* 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); 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 24937) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24938) @@ -85,4 +85,5 @@ Vector *N_gia_rate= NULL; Vector *U_gia_rate= NULL; + SealevelMasks* masks=NULL; /*parameters:*/ @@ -148,12 +149,15 @@ if(modelid==earthid){ + /*call masks core: */ + masks=sealevelrise_core_masks(femmodel); + /*call eustatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */ - RSLg_eustatic=sealevelrise_core_eustatic(femmodel,&eartharea,&oceanarea); + RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&eartharea,&oceanarea); /*call non-eustatic core (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */ - RSLg=sealevelrise_core_noneustatic(femmodel,RSLg_eustatic,eartharea,oceanarea); + RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,eartharea,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); + sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,eartharea,masks); /*compute viscosus (GIA) geodetic signatures:*/ @@ -223,4 +227,5 @@ delete N_gia_rate; delete U_gia_rate; + //delete masks; /*}}}*/ @@ -307,4 +312,20 @@ } /*}}}*/ +SealevelMasks* sealevelrise_core_masks(FemModel* femmodel) { /*{{{*/ + + if(VerboseSolution()) _printf0_(" computing sea level masks\n"); + + /*initialize SealevelMasks structure: */ + SealevelMasks* masks=new SealevelMasks(femmodel->elements->Size()); + + /*go through elements and fill the masks: */ + for (int i=0;ielements->Size();i++){ + + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->SetSealevelMasks(masks); + } + + return masks; +}/*}}}*/ void sealevelrise_core_geometry(FemModel* femmodel) { /*{{{*/ @@ -333,5 +354,5 @@ }/*}}}*/ -Vector* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/ +Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/ /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ @@ -367,8 +388,8 @@ /*call the eustatic main module: */ - femmodel->SealevelriseEustatic(RSLgi,&eartharea, &oceanarea,&eustatic, latitude, longitude, radius,loop); //this computes + femmodel->SealevelriseEustatic(RSLgi,&eartharea, &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: */ - RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,oceanarea); + RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,masks, oceanarea); /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the @@ -389,5 +410,5 @@ return RSLgi; }/*}}}*/ -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea){ /*{{{*/ +Vector* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea){ /*{{{*/ /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. @@ -458,5 +479,5 @@ /*call the non eustatic module: */ - femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, eartharea,latitude,longitude, radius,verboseconvolution,loop); + femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, eartharea, masks, latitude,longitude, radius,verboseconvolution,loop); /*assemble solution vector: */ @@ -467,5 +488,5 @@ /*call rotational feedback module: */ RSLgo_rot = new Vector(gsize); RSLgo_rot->Assemble(); - femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz,eartharea,latitude,longitude,radius); + femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, eartharea, masks, latitude,longitude,radius); RSLgo_rot->Assemble(); @@ -479,5 +500,5 @@ /*we need to average RSLgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */ - RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo,oceanarea); + RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo,masks, oceanarea); /*RSLg is the sum of the eustatic term, and the ocean terms: */ @@ -516,5 +537,5 @@ return RSLg; } /*}}}*/ -void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg,IssmDouble eartharea){ /*{{{*/ +void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg,IssmDouble eartharea, SealevelMasks* masks){ /*{{{*/ Vector *U_esa = NULL; @@ -556,5 +577,5 @@ /*call the elastic main modlule:*/ - femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,eartharea,latitude,longitude,radius,xx,yy,zz,loop,horiz); + femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,eartharea, masks, latitude,longitude,radius,xx,yy,zz,loop,horiz); /*Assign output pointers:*/