Index: /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 21907) +++ /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 21908) @@ -80,6 +80,13 @@ U_elastic_local[i-lower_row]= (love_h[nl-1])/2.0/sin(alpha/2.0); H_elastic_local[i-lower_row]= 0; - IssmDouble Pn,Pn1,Pn2; - IssmDouble Pn_p,Pn_p1,Pn_p2; + //IssmDouble Pn,Pn1,Pn2; + //IssmDouble Pn_p,Pn_p1,Pn_p2; + IssmDouble Pn = 0.; + IssmDouble Pn1 = 0.; + IssmDouble Pn2 = 0.; + IssmDouble Pn_p = 0.; + IssmDouble Pn_p1 = 0.; + IssmDouble Pn_p2 = 0.; + for (int n=0;n* pSgo_rot, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ +void FemModel::SealevelriseRotationalFeedback(Vector* pSgo_rot, Vector* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ /*serialized vectors:*/ @@ -3996,4 +3996,9 @@ pSgo_rot->Assemble(); + /*Assign output pointers:*/ + *pIxz=moi_list[0]; + *pIyz=moi_list[1]; + *pIzz=moi_list[2]; + /*Free ressources:*/ xDelete(Sg_old); Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21907) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 21908) @@ -137,5 +137,5 @@ void SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution); - void SealevelriseRotationalFeedback(Vector* pSgo_rot, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); + void SealevelriseRotationalFeedback(Vector* pSgo_rot, Vector* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); void SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); IssmDouble SealevelriseOceanAverage(Vector* Sg); Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 21907) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 21908) @@ -83,6 +83,4 @@ Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //ocean loading tems (2nd and 5th terms on the RHS of Farrel and Clark) - - /*compute other geodetic signatures, such as absolute sea level chagne, components of 3-D crustal motion: */ /*Initialize:*/ Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp (revision 21907) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core_noneustatic.cpp (revision 21908) @@ -37,4 +37,7 @@ IssmDouble *radius = NULL; IssmDouble eustatic; + IssmDouble Ixz; + IssmDouble Iyz; + IssmDouble Izz; /*Recover some parameters: */ @@ -83,8 +86,14 @@ /*call rotational feedback module: */ Sgo_rot = new Vector(gsize); Sgo_rot->Assemble(); - femmodel->SealevelriseRotationalFeedback(Sgo_rot,Sg_old,latitude,longitude,radius); + femmodel->SealevelriseRotationalFeedback(Sgo_rot,Sg_old,&Ixz,&Iyz,&Izz,latitude,longitude,radius); Sgo_rot->Assemble(); + /*save changes in inertia tensor as results: */ + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,Izz)); + Sgo->AXPY(Sgo_rot,1); + delete Sgo_rot; } @@ -93,5 +102,5 @@ /*Sg is the sum of the eustatic term, and the ocean terms: */ - Sg_eustatic->Copy(Sg); Sg->AXPY(Sgo,1); + Sg_eustatic->Copy(Sg); Sg->AXPY(Sgo,1); delete Sgo; Sg->Shift(-Sgo_oceanaverage); @@ -112,5 +121,4 @@ /*some minor verbosing adjustment:*/ if(count>1)verboseconvolution=false; - } if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n"); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21907) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21908) @@ -795,4 +795,7 @@ SealevelAbsoluteEnum, SealevelEustaticEnum, + SealevelInertiaTensorXZEnum, + SealevelInertiaTensorYZEnum, + SealevelInertiaTensorZZEnum, SealevelObsEnum, SealevelWeightsEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21907) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21908) @@ -775,4 +775,7 @@ case SealevelAbsoluteEnum : return "SealevelAbsolute"; case SealevelEustaticEnum : return "SealevelEustatic"; + case SealevelInertiaTensorXZEnum : return "SealevelInertiaTensorXZ"; + case SealevelInertiaTensorYZEnum : return "SealevelInertiaTensorYZ"; + case SealevelInertiaTensorZZEnum : return "SealevelInertiaTensorZZ"; case SealevelObsEnum : return "SealevelObs"; case SealevelWeightsEnum : return "SealevelWeights"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21907) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21908) @@ -793,4 +793,7 @@ else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum; else if (strcmp(name,"SealevelEustatic")==0) return SealevelEustaticEnum; + else if (strcmp(name,"SealevelInertiaTensorXZ")==0) return SealevelInertiaTensorXZEnum; + else if (strcmp(name,"SealevelInertiaTensorYZ")==0) return SealevelInertiaTensorYZEnum; + else if (strcmp(name,"SealevelInertiaTensorZZ")==0) return SealevelInertiaTensorZZEnum; else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum; else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum; @@ -872,11 +875,11 @@ else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; else if (strcmp(name,"Matestar")==0) return MatestarEnum; - else if (strcmp(name,"Matpar")==0) return MatparEnum; - else if (strcmp(name,"Node")==0) return NodeEnum; - else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; + if (strcmp(name,"Matpar")==0) return MatparEnum; + else if (strcmp(name,"Node")==0) return NodeEnum; + else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; + else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; else if (strcmp(name,"Param")==0) return ParamEnum; else if (strcmp(name,"Moulin")==0) return MoulinEnum; @@ -995,11 +998,11 @@ else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum; else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; - else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum; - else if (strcmp(name,"IceMass")==0) return IceMassEnum; - else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; + if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum; + else if (strcmp(name,"IceMass")==0) return IceMassEnum; + else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; + else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum; else if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum;