Index: /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 24958) +++ /issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp (revision 24959) @@ -47,4 +47,6 @@ int M,m,lower_row,upper_row; IssmDouble degacc=.01; + IssmDouble planetradius=0; + IssmDouble planetarea=0; int numoutputs; @@ -59,10 +61,15 @@ /*some constant parameters: */ parameters->AddObject(iomodel->CopyConstantObject("md.esa.hemisphere",EsaHemisphereEnum)); - parameters->AddObject(iomodel->CopyConstantObject("md.slr.eartharea",SealevelEarthAreaEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.horiz",SealevelriseHorizEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum)); - + + /*deal with planet radius and area: */ + parameters->AddObject(iomodel->CopyConstantObject("md.slr.planetradius",SealevelPlanetRadiusEnum)); + iomodel->FetchData(&planetradius,"md.slr.planetradius"); + planetarea=4*PI*planetradius*planetradius; + parameters->AddObject(new DoubleParam(SealevelPlanetAreaEnum,planetarea)); + /*love numbers: */ iomodel->FetchData(&love_h,&nl,NULL,"md.esa.love_h"); Index: /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24958) +++ /issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24959) @@ -1,3 +1,4 @@ #include "./SealevelriseAnalysis.h" +#include #include "../toolkits/toolkits.h" #include "../classes/classes.h" @@ -292,4 +293,6 @@ int M,m,lower_row,upper_row; IssmDouble degacc=.01; + IssmDouble planetradius=0; + IssmDouble planetarea=0; int numoutputs; @@ -320,5 +323,10 @@ 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)); + parameters->AddObject(iomodel->CopyConstantObject("md.slr.planetradius",SealevelPlanetRadiusEnum)); + + /*compute planet area and plug into parameters:*/ + iomodel->FetchData(&planetradius,"md.slr.planetradius"); + planetarea=4*PI*planetradius*planetradius; + parameters->AddObject(new DoubleParam(SealevelPlanetAreaEnum,planetarea)); /*Deal with dsl multi-model ensembles: {{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24958) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24959) @@ -5297,5 +5297,5 @@ IssmDouble llr_list[NUMVERTICES][3]; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble area,eartharea; + IssmDouble area,planetarea; IssmDouble I; //ice/water loading IssmDouble late,longe,re; @@ -5331,5 +5331,5 @@ /*recover earth area: */ - this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + this->parameters->FindParam(&planetarea,SealevelPlanetAreaEnum); /*how many dofs are we working with here? */ @@ -5356,10 +5356,12 @@ } - // correction at the north pole + // correction at the north pole: given longitude of the North pole a definition + // closer to the other two vertices. 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 + // correction at the north pole: given longitude of the North pole a definition + // closer to the other two vertices. 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; @@ -5370,5 +5372,5 @@ late=90-late; - if(longe>180)longe=(longe-180)-180; + if(longe>180)longe=longe-360; late=late/180*PI; @@ -5432,7 +5434,7 @@ /*Add all components to the pUp solution vectors:*/ - U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i]; - N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i]; - E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i]; + U_values[i]+=3*rho_ice/rho_earth*area/planetarea*I*U_elastic[i]; + N_values[i]+=3*rho_ice/rho_earth*area/planetarea*I*N_elastic[i]; + E_values[i]+=3*rho_ice/rho_earth*area/planetarea*I*E_elastic[i]; } pUp->SetValues(gsize,indices,U_values,ADD_VAL); @@ -5479,9 +5481,9 @@ /*Compute area of element:*/ - IssmDouble area,eartharea; + IssmDouble area,planetarea; this->GetInput2Value(&area,AreaEnum); /*recover earth area: */ - this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + this->parameters->FindParam(&planetarea,SealevelPlanetAreaEnum); /*Compute lat,long,radius of elemental centroid: */ @@ -5541,7 +5543,7 @@ * 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; + dI_list[0] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; + dI_list[1] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; + dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; } else if(masks->isiceonly[this->lid]){ @@ -5556,7 +5558,7 @@ 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; + dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; + dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; + dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; } @@ -5585,5 +5587,5 @@ IssmDouble llr_list[NUMVERTICES][3]; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble area,eartharea; + IssmDouble area,planetarea,planetradius; IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) IssmDouble rho_earth; @@ -5617,5 +5619,6 @@ this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); - this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + this->parameters->FindParam(&planetarea,SealevelPlanetAreaEnum); + this->parameters->FindParam(&planetradius,SealevelPlanetRadiusEnum); this->parameters->FindParam(&horiz,SealevelriseHorizEnum); @@ -5645,50 +5648,20 @@ area=GetAreaSpherical(); + + /* Where is the centroid of this element:*/ + + /*retrieve coordinates: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + /*figure out gravity center of our element (Cartesian): */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; - /* Where is the centroid of this element?:{{{*/ - - /*retrieve coordinates: */ - ::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; - /*}}}*/ - - constant=3/rho_earth*area/eartharea; + /*compute gravity center in lat,long: */ + late= asin(z_element/planetradius); + longe = atan2(y_element,x_element); + + constant=3/rho_earth*area/planetarea; for(int i=0;iPI)delLambda=2*PI-delLambda; alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); indices[i]=alpha/PI*reCast(M-1); @@ -5904,5 +5877,5 @@ bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area,eartharea; + IssmDouble area,planetarea; IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) IssmDouble rho; @@ -5939,5 +5912,5 @@ /*recover earth area: */ - this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + this->parameters->FindParam(&planetarea,SealevelPlanetAreaEnum); /*recover love numbers and computational flags: */ @@ -6038,5 +6011,5 @@ /*Add all components to the pSgi or pSgo solution vectors:*/ - Sgi[i]+=3*rho_water/rho_earth*area/eartharea*I*(G_rigid+G_elastic); + Sgi[i]+=3*rho_water/rho_earth*area/planetarea*I*(G_rigid+G_elastic); } } Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24958) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24959) @@ -302,5 +302,6 @@ syn keyword cConstant RootPathEnum syn keyword cConstant SaveResultsEnum -syn keyword cConstant SealevelEarthAreaEnum +syn keyword cConstant SealevelPlanetRadiusEnum +syn keyword cConstant SealevelPlanetAreaEnum syn keyword cConstant SealevelEustaticEnum syn keyword cConstant SealevelriseAbstolEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24958) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24959) @@ -296,5 +296,6 @@ RootPathEnum, SaveResultsEnum, - SealevelEarthAreaEnum, + SealevelPlanetRadiusEnum, + SealevelPlanetAreaEnum, SealevelEustaticEnum, SealevelriseAbstolEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24958) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24959) @@ -304,5 +304,6 @@ case RootPathEnum : return "RootPath"; case SaveResultsEnum : return "SaveResults"; - case SealevelEarthAreaEnum : return "SealevelEarthArea"; + case SealevelPlanetRadiusEnum : return "SealevelPlanetRadius"; + case SealevelPlanetAreaEnum : return "SealevelPlanetArea"; 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 24958) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24959) @@ -310,5 +310,6 @@ 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,"SealevelPlanetRadius")==0) return SealevelPlanetRadiusEnum; + else if (strcmp(name,"SealevelPlanetArea")==0) return SealevelPlanetAreaEnum; 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,"SmbGsp")==0) return SmbGspEnum; else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum; - else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; + if (strcmp(name,"SmbHref")==0) return SmbHrefEnum; + else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; else if (strcmp(name,"SmbMassBalanceClimate")==0) return SmbMassBalanceClimateEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Waterheight")==0) return WaterheightEnum; else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum; - else if (strcmp(name,"Outputdefinition1")==0) return Outputdefinition1Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition10")==0) return Outputdefinition10Enum; + if (strcmp(name,"Outputdefinition1")==0) return Outputdefinition1Enum; + else if (strcmp(name,"Outputdefinition10")==0) return Outputdefinition10Enum; else if (strcmp(name,"Outputdefinition11")==0) return Outputdefinition11Enum; else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum; @@ -997,9 +998,9 @@ else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum; else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum; - else if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum; + if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum; + else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum; else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum; else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum; else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum; - else if (strcmp(name,"Hydrologypism")==0) return HydrologypismEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Hydrologyshakti")==0) return HydrologyshaktiEnum; + if (strcmp(name,"Hydrologypism")==0) return HydrologypismEnum; + else if (strcmp(name,"Hydrologyshakti")==0) return HydrologyshaktiEnum; else if (strcmp(name,"Hydrologyshreve")==0) return HydrologyshreveEnum; else if (strcmp(name,"IceMass")==0) return IceMassEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum; else if (strcmp(name,"P2xP1")==0) return P2xP1Enum; - else if (strcmp(name,"P2xP4")==0) return P2xP4Enum; else stage=11; } if(stage==11){ - if (strcmp(name,"Paterson")==0) return PatersonEnum; + if (strcmp(name,"P2xP4")==0) return P2xP4Enum; + else if (strcmp(name,"Paterson")==0) return PatersonEnum; else if (strcmp(name,"Pengrid")==0) return PengridEnum; else if (strcmp(name,"Penpair")==0) return PenpairEnum; @@ -1366,9 +1367,9 @@ else if (strcmp(name,"StressbalanceViscosityOvershoot")==0) return StressbalanceViscosityOvershootEnum; else if (strcmp(name,"SubelementMigration4")==0) return SubelementMigration4Enum; - else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; + if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum; + else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum; else stage=13;