Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Element.h =================================================================== --- /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Element.h (revision 21391) +++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Element.h (revision 21392) @@ -301,4 +301,5 @@ 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(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; Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Penta.h =================================================================== --- /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Penta.h (revision 21391) +++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Penta.h (revision 21392) @@ -186,4 +186,5 @@ 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!");}; 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!");}; Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Seg.h =================================================================== --- /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Seg.h (revision 21391) +++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Seg.h (revision 21392) @@ -170,4 +170,5 @@ #ifdef _HAVE_SEALEVELRISE_ + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_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!");}; Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tetra.h (revision 21391) +++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tetra.h (revision 21392) @@ -176,4 +176,5 @@ #endif #ifdef _HAVE_SEALEVELRISE_ + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_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!");}; Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tria.cpp (revision 21391) +++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tria.cpp (revision 21392) @@ -3611,4 +3611,91 @@ } /*}}}*/ +void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){/*{{{*/ + + /*early return if we are not on an ice cap OR ocean:*/ + if(!(this->inputs->Max(MaskIceLevelsetEnum)<0) && !IsWaterInElement()) return; + + /*Compute area of element:*/ + IssmDouble area; + area=GetAreaSpherical(); + + /*Compute lat,long,radius of elemental centroid: */ + bool spherical=true; + IssmDouble llr_list[NUMVERTICES][3]; + IssmDouble late,longe,re; + /* Where is the centroid of this element?:{{{*/ + ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); + + IssmDouble minlong=400; + IssmDouble maxlong=-20; + for (int i=0;imaxlong)maxlong=llr_list[i][1]; + if(llr_list[i][1]180){ + if (llr_list[0][1]==0)llr_list[0][1]=360; + if (llr_list[1][1]==0)llr_list[1][1]=360; + if (llr_list[2][1]==0)llr_list[2][1]=360; + } + + // correction at the north pole + if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; + if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; + if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; + + //correction at the south pole + if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; + if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; + if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; + + late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0; + longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0; + + late=90-late; + if(longe>180)longe=(longe-180)-180; + + late=late/180*PI; + longe=longe/180*PI; + /*}}}*/ + re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0; + + if(IsWaterInElement()){ + IssmDouble rho_water, S; + + /*recover material parameters: */ + rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + + /*From Sg_old, recover water sea level rise:*/ + S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; + + /* Perturbation terms for moment of inertia (moi_list): + * computed analytically (see Wu & Peltier, eqs 10 & 32) + * also consistent with my GMD formulation! + * ALL in geographic coordinates + * */ + dI_list[0] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea; + dI_list[1] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea; + dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea; + } + else if(this->inputs->Max(MaskIceLevelsetEnum)<0){ + IssmDouble rho_ice, I; + + /*recover material parameters: */ + rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum); + + /*Compute ice thickness change: */ + Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); + if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); + deltathickness_input->GetInputAverage(&I); + + dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea; + dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea; + dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea; + } + + return; +}/*}}}*/ void Tria::SealevelriseEustatic(Vector* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tria.h =================================================================== --- /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tria.h (revision 21391) +++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tria.h (revision 21392) @@ -147,4 +147,5 @@ IssmDouble OceanArea(void); IssmDouble OceanAverage(IssmDouble* Sg); + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea); 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); Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/EnumDefinitions.h (revision 21391) +++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/EnumDefinitions.h (revision 21392) @@ -772,4 +772,11 @@ SealevelriseRigidEnum, SealevelriseElasticEnum, + SealevelriseRotationEnum, + SealevelriseTidalLoveHEnum, + SealevelriseTidalLoveKEnum, + SealevelriseFluidLoveEnum, + SealevelriseEquatorialMoiEnum, + SealevelrisePolarMoiEnum, + SealevelriseAngularVelocityEnum, SealevelriseGElasticEnum, SealevelriseUElasticEnum, Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/EnumToStringx.cpp (revision 21391) +++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/EnumToStringx.cpp (revision 21392) @@ -753,4 +753,11 @@ case SealevelriseRigidEnum : return "SealevelriseRigid"; case SealevelriseElasticEnum : return "SealevelriseElastic"; + case SealevelriseRotationEnum : return "SealevelriseRotation"; + case SealevelriseTidalLoveHEnum : return "SealevelriseTidalLoveH"; + case SealevelriseTidalLoveKEnum : return "SealevelriseTidalLoveK"; + case SealevelriseFluidLoveEnum : return "SealevelriseFluidLove"; + case SealevelriseEquatorialMoiEnum : return "SealevelriseEquatorialMoi"; + case SealevelrisePolarMoiEnum : return "SealevelrisePolarMoi"; + case SealevelriseAngularVelocityEnum : return "SealevelriseAngularVelocity"; case SealevelriseGElasticEnum : return "SealevelriseGElastic"; case SealevelriseUElasticEnum : return "SealevelriseUElastic"; Index: /issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/StringToEnumx.cpp (revision 21391) +++ /issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/StringToEnumx.cpp (revision 21392) @@ -771,4 +771,11 @@ else if (strcmp(name,"SealevelriseRigid")==0) return SealevelriseRigidEnum; else if (strcmp(name,"SealevelriseElastic")==0) return SealevelriseElasticEnum; + else if (strcmp(name,"SealevelriseRotation")==0) return SealevelriseRotationEnum; + else if (strcmp(name,"SealevelriseTidalLoveH")==0) return SealevelriseTidalLoveHEnum; + else if (strcmp(name,"SealevelriseTidalLoveK")==0) return SealevelriseTidalLoveKEnum; + else if (strcmp(name,"SealevelriseFluidLove")==0) return SealevelriseFluidLoveEnum; + else if (strcmp(name,"SealevelriseEquatorialMoi")==0) return SealevelriseEquatorialMoiEnum; + else if (strcmp(name,"SealevelrisePolarMoi")==0) return SealevelrisePolarMoiEnum; + else if (strcmp(name,"SealevelriseAngularVelocity")==0) return SealevelriseAngularVelocityEnum; else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum; else if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum; @@ -868,5 +875,8 @@ else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum; else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum; - else if (strcmp(name,"HydrologySommersAnalysis")==0) return HydrologySommersAnalysisEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"HydrologySommersAnalysis")==0) return HydrologySommersAnalysisEnum; else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum; @@ -875,8 +885,5 @@ else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; + else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; Index: /issm/branches/trunk-larour-NatClimateChange2016/src/m/classes/slr.m =================================================================== --- /issm/branches/trunk-larour-NatClimateChange2016/src/m/classes/slr.m (revision 21391) +++ /issm/branches/trunk-larour-NatClimateChange2016/src/m/classes/slr.m (revision 21392) @@ -16,4 +16,8 @@ tide_love_k = 0; %ideam tide_love_h = 0; %ideam + fluid_love = 0; + equatorial_moi = 0; + polar_moi = 0; + angular_velocity = 0; rigid = 0; elastic = 0; @@ -44,9 +48,19 @@ self.rigid=1; self.elastic=1; - self.rotation=1; + self.rotation=0; %tidal love numbers: self.tide_love_h=0.6149; %degree 2 self.tide_love_k=0.3055; % degree 2 + + %secular fluid love number: + self.fluid_love=0.942; + + %moment of inertia: + self.equatorial_moi=8.0077*10^37; % [kg m^2] + self.polar_moi =8.0345*10^37; % [kg m^2] + + % mean rotational velocity of earth + self.angular_velocity=7.2921*10^-5; % [s^-1] %numerical discretization accuracy @@ -70,4 +84,8 @@ md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','slr.fluid_love','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','slr.equatorial_moi','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','slr.polar_moi','NaN',1,'Inf',1); + md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.reltol','size',[1 1]); md = checkfield(md,'fieldname','slr.abstol','size',[1 1]); @@ -106,7 +124,11 @@ 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,'fluid_love','secular fluid Love number'); + fielddisplay(self,'equatorial_moi','mean equatorial moment of inertia [kg m^2]'); + fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]'); + fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]'); 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,'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'); @@ -125,4 +147,8 @@ WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','fluid_love','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','equatorial_moi','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','polar_moi','format','Double'); + WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean'); WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean'); @@ -153,7 +179,11 @@ writejsdouble(fid,[modelname '.slr.tide_love_k'],self.tide_love_k); writejsdouble(fid,[modelname '.slr.tide_love_h'],self.tide_love_h); + writejsdouble(fid,[modelname '.slr.fluid_love'],self.fluid_love); + writejsdouble(fid,[modelname '.slr.equatorial_moi'],self.equatorial_moi); + writejsdouble(fid,[modelname '.slr.polar_moi'],self.polar_moi); + writejsdouble(fid,[modelname '.slr.angular_velocity'],self.angular_velocity); 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.degacc'],self.degacc); writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs);