Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 21294) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 21295) @@ -70,4 +70,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum)); iomodel->FetchData(&elastic,"md.slr.elastic"); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 21294) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 21295) @@ -306,6 +306,6 @@ virtual IssmDouble OceanArea(void)=0; virtual void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; - virtual void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; - virtual void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea)=0; + virtual void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0; + virtual void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 21294) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 21295) @@ -191,6 +191,6 @@ IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_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 21294) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 21295) @@ -174,6 +174,6 @@ #ifdef _HAVE_SEALEVELRISE_ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21294) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 21295) @@ -181,6 +181,6 @@ #ifdef _HAVE_SEALEVELRISE_ void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; IssmDouble OceanArea(void){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21294) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21295) @@ -3895,4 +3895,5 @@ bool computerigid = true; bool computeelastic= true; + bool scaleoceanarea= true; /*early return if we are not on an ice cap:*/ @@ -3910,4 +3911,5 @@ this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); + this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum); /*recover elastic green function:*/ @@ -3970,4 +3972,5 @@ /*Compute eustatic compoent:*/ _assert_(oceanarea>0.); + if(scaleoceanarea) oceanarea=3.6e+14; // use true ocean area, m^2 eustatic += rho_ice*area*I/(oceanarea*rho_water); @@ -4015,5 +4018,5 @@ } /*}}}*/ -void Tria::SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/ /*diverse:*/ @@ -4158,5 +4161,5 @@ } /*}}}*/ -void Tria::SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/ /*diverse:*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 21294) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 21295) @@ -152,6 +152,6 @@ IssmDouble OceanAverage(IssmDouble* Sg); void SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble oceanarea,IssmDouble eartharea); + void SealevelriseNonEustatic(Vector* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea); + void SealevelriseGeodetic(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21294) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 21295) @@ -2449,6 +2449,4 @@ IssmDouble* Sg_old=NULL; - IssmDouble oceanarea=0; - IssmDouble oceanarea_cpu=0; IssmDouble eartharea=0; IssmDouble eartharea_cpu=0; @@ -2465,10 +2463,6 @@ for(int i=0;i(elements->GetObjectByOffset(i)); - oceanarea_cpu += element->OceanArea(); eartharea_cpu += element->GetAreaSpherical(); } - ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); @@ -2485,5 +2479,5 @@ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseNonEustatic(pSgo, Sg_old,latitude,longitude,radius,oceanarea,eartharea); + element->SealevelriseNonEustatic(pSgo, Sg_old,latitude,longitude,radius,eartharea); } if(i%100==0)pSgo->Assemble(); @@ -2500,6 +2494,4 @@ IssmDouble* Sg=NULL; - IssmDouble oceanarea=0; - IssmDouble oceanarea_cpu=0; IssmDouble eartharea=0; IssmDouble eartharea_cpu=0; @@ -2516,10 +2508,6 @@ for(int i=0;i(elements->GetObjectByOffset(i)); - oceanarea_cpu += element->OceanArea(); eartharea_cpu += element->GetAreaSpherical(); } - ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); @@ -2533,5 +2521,5 @@ if(i(elements->GetObjectByOffset(i)); - element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,oceanarea,eartharea); + element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea); } if(i%100==0){ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21294) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21295) @@ -773,4 +773,5 @@ SealevelriseRigidEnum, SealevelriseElasticEnum, + SealevelriseOceanAreaScalingEnum, SealevelriseGElasticEnum, SealevelriseUElasticEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21294) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21295) @@ -754,4 +754,5 @@ case SealevelriseRigidEnum : return "SealevelriseRigid"; case SealevelriseElasticEnum : return "SealevelriseElastic"; + case SealevelriseOceanAreaScalingEnum : return "SealevelriseOceanAreaScaling"; case SealevelriseGElasticEnum : return "SealevelriseGElastic"; case SealevelriseUElasticEnum : return "SealevelriseUElastic"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21294) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21295) @@ -772,4 +772,5 @@ else if (strcmp(name,"SealevelriseRigid")==0) return SealevelriseRigidEnum; else if (strcmp(name,"SealevelriseElastic")==0) return SealevelriseElasticEnum; + else if (strcmp(name,"SealevelriseOceanAreaScaling")==0) return SealevelriseOceanAreaScalingEnum; else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum; else if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum; else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum; - else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum; + if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; + else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum; else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum; else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; Index: /issm/trunk-jpl/src/m/classes/slr.js =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.js (revision 21294) +++ /issm/trunk-jpl/src/m/classes/slr.js (revision 21295) @@ -19,4 +19,5 @@ this.elastic=1; this.rotation=1; + this.ocean_area_scaling=1; //tidal love numbers: @@ -81,4 +82,5 @@ fielddisplay(this,'elastic','elastic earth graviational potential perturbation'); fielddisplay(this,'rotation','rotational earth potential perturbation'); + fielddisplay(this,'ocean_area_scaling','correction for model representation of ocean area'); fielddisplay(this,'degacc',"accuracy (default .01 deg) for numerical discretization of the Green's functions"); fielddisplay(this,'transitions','indices into parts of the mesh that will be icecaps'); @@ -100,4 +102,5 @@ WriteData(fid,prefix,'object',this,'fieldname','elastic','format','Boolean'); WriteData(fid,prefix,'object',this,'fieldname','rotation','format','Boolean'); + WriteData(fid,prefix,'object',this,'fieldname','ocean_area_scaling','format','Boolean'); WriteData(fid,prefix,'object',this,'fieldname','degacc','format','Double'); WriteData(fid,prefix,'object',this,'fieldname','transitions','format','MatArray'); @@ -128,4 +131,5 @@ this.elastic=NullFix(this.elastic,NaN); this.rotation=NullFix(this.rotation,NaN); + this.ocean_area_scaling=NullFix(this.ocean_area_scaling,NaN); this.degacc=NullFix(this.degacc,NaN); }//}}} @@ -145,4 +149,5 @@ this.elastic = 0; this.rotation = 0; + this.ocean_area_scaling = 0; this.degacc = 0; this.requested_outputs = []; Index: /issm/trunk-jpl/src/m/classes/slr.m =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.m (revision 21294) +++ /issm/trunk-jpl/src/m/classes/slr.m (revision 21295) @@ -19,4 +19,5 @@ elastic = 0; rotation = 0; + ocean_area_scaling = 0; degacc = 0; requested_outputs = {}; @@ -45,4 +46,5 @@ self.elastic=1; self.rotation=1; + self.ocean_area_scaling=1; %tidal love numbers: @@ -106,7 +108,8 @@ fielddisplay(self,'tide_love_k','tidal load Love number (deg 2)'); fielddisplay(self,'tide_love_h','tidal load Love number (deg 2)'); - fielddisplay(self,'rotation','earth rotational potential perturbation'); fielddisplay(self,'rigid','rigid earth graviational potential perturbation'); fielddisplay(self,'elastic','elastic earth graviational potential perturbation'); + fielddisplay(self,'rotation','earth rotational potential perturbation'); + fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area'); fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'); fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); @@ -128,4 +131,5 @@ WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean'); + WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray'); @@ -154,6 +158,7 @@ writejsdouble(fid,[modelname '.slr.tide_love_h'],self.tide_love_h); writejsdouble(fid,[modelname '.slr.rigid'],self.rigid); + writejsdouble(fid,[modelname '.slr.elastic'],self.elastic); writejsdouble(fid,[modelname '.slr.rotation'],self.rotation); - writejsdouble(fid,[modelname '.slr.elastic'],self.elastic); + writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling); writejsdouble(fid,[modelname '.slr.degacc'],self.degacc); writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs); Index: /issm/trunk-jpl/src/m/classes/slr.py =================================================================== --- /issm/trunk-jpl/src/m/classes/slr.py (revision 21294) +++ /issm/trunk-jpl/src/m/classes/slr.py (revision 21295) @@ -28,4 +28,5 @@ self.elastic = 0 self.rotation = 0 + self.ocean_area_scaling = 0; self.degacc = 0 self.requested_outputs = [] @@ -49,4 +50,5 @@ string="%s\n%s"%(string,fielddisplay(self,'elastic','elastic earth graviational potential perturbation')) string="%s\n%s"%(string,fielddisplay(self,'rotation','earth rotational potential perturbation')) + string="%s\n%s"%(string,fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area')) string="%s\n%s"%(string,fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions')) string="%s\n%s"%(string,fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps')) @@ -68,4 +70,5 @@ self.elastic=1 self.rotation=1 + self.ocean_area_scaling=1 #tidal love numbers: @@ -127,4 +130,5 @@ WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean') WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean') + WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean') WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double') WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray')