Index: ../trunk-jpl/src/m/classes/slr.m =================================================================== --- ../trunk-jpl/src/m/classes/slr.m (revision 24939) +++ ../trunk-jpl/src/m/classes/slr.m (revision 24940) @@ -32,6 +32,7 @@ horiz = 0; Ngia = NaN; Ugia = NaN; + eartharea = 4*pi*planetradius('earth')^2; requested_outputs = {}; transitions = {}; @@ -101,7 +102,8 @@ return; end - md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); + %md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]); + md = checkfield(md,'fieldname','slr.deltathickness','timeseries',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1); md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1); @@ -131,12 +133,12 @@ end %cross check that whereever we have an ice load, the mask is <0 on each vertex: - pos=find(self.deltathickness); - maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:)); - [els,vertices]=find(maskpos>0); - if length(els), - warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!'); - end + %pos=find(self.deltathickness); + %maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:)); + %[els,vertices]=find(maskpos>0); + %if length(els), + % warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!'); + %end %check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not, %a coupler to a planet model is provided. @@ -220,6 +222,7 @@ WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','eartharea','format','Double'); %process requested outputs outputs = self.requested_outputs; Index: ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24939) +++ ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 24940) @@ -320,6 +320,7 @@ parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum)); 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)); /*Deal with dsl multi-model ensembles: {{{*/ iomodel->FetchData(&dslmodel,"md.dsl.model"); Index: ../trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24939) +++ ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24940) @@ -301,6 +301,7 @@ syn keyword cConstant ResultsEnum syn keyword cConstant RootPathEnum syn keyword cConstant SaveResultsEnum +syn keyword cConstant SealevelEarthAreaEnum syn keyword cConstant SealevelEustaticEnum syn keyword cConstant SealevelriseAbstolEnum syn keyword cConstant SealevelriseAngularVelocityEnum @@ -1430,6 +1431,7 @@ syn keyword cType Results syn keyword cType RiftStruct syn keyword cType Riftfront +syn keyword cType SealevelMasks syn keyword cType Seg syn keyword cType SegInput2 syn keyword cType SegRef Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24939) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24940) @@ -295,6 +295,7 @@ ResultsEnum, RootPathEnum, SaveResultsEnum, + SealevelEarthAreaEnum, SealevelEustaticEnum, SealevelriseAbstolEnum, SealevelriseAngularVelocityEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24939) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24940) @@ -303,6 +303,7 @@ case ResultsEnum : return "Results"; case RootPathEnum : return "RootPath"; case SaveResultsEnum : return "SaveResults"; + case SealevelEarthAreaEnum : return "SealevelEarthArea"; case SealevelEustaticEnum : return "SealevelEustatic"; case SealevelriseAbstolEnum : return "SealevelriseAbstol"; case SealevelriseAngularVelocityEnum : return "SealevelriseAngularVelocity"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24939) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24940) @@ -309,6 +309,7 @@ else if (strcmp(name,"Results")==0) return ResultsEnum; 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,"SealevelEustatic")==0) return SealevelEustaticEnum; else if (strcmp(name,"SealevelriseAbstol")==0) return SealevelriseAbstolEnum; else if (strcmp(name,"SealevelriseAngularVelocity")==0) return SealevelriseAngularVelocityEnum; @@ -381,11 +382,11 @@ else if (strcmp(name,"SmbK")==0) return SmbKEnum; 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; else if (strcmp(name,"SmbRunoffalti")==0) return SmbRunoffaltiEnum; @@ -504,11 +505,11 @@ else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; 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; else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; @@ -627,11 +628,11 @@ else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; 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; else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; @@ -750,11 +751,11 @@ else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; else if (strcmp(name,"SmbMassBalanceClimate")==0) return SmbMassBalanceClimateEnum; - else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; + if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; + else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum; else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum; else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum; @@ -873,11 +874,11 @@ else if (strcmp(name,"Outputdefinition10")==0) return Outputdefinition10Enum; else if (strcmp(name,"Outputdefinition11")==0) return Outputdefinition11Enum; else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum; - else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum; + if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum; + else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum; else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum; else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum; @@ -996,11 +997,11 @@ else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum; else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum; else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum; - else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; + if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum; + else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum; @@ -1119,11 +1120,11 @@ else if (strcmp(name,"Hydrologyshakti")==0) return HydrologyshaktiEnum; else if (strcmp(name,"Hydrologyshreve")==0) return HydrologyshreveEnum; else if (strcmp(name,"IceMass")==0) return IceMassEnum; - else if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; + if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum; + else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum; else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum; @@ -1242,11 +1243,11 @@ else if (strcmp(name,"Paterson")==0) return PatersonEnum; else if (strcmp(name,"Pengrid")==0) return PengridEnum; else if (strcmp(name,"Penpair")==0) return PenpairEnum; - else if (strcmp(name,"Penta")==0) return PentaEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"PentaInput")==0) return PentaInputEnum; + if (strcmp(name,"Penta")==0) return PentaEnum; + else if (strcmp(name,"PentaInput")==0) return PentaInputEnum; else if (strcmp(name,"Profiler")==0) return ProfilerEnum; else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum; else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; Index: ../trunk-jpl/src/c/cores/cores.h =================================================================== --- ../trunk-jpl/src/c/cores/cores.h (revision 24939) +++ ../trunk-jpl/src/c/cores/cores.h (revision 24940) @@ -57,9 +57,9 @@ void steric_core(FemModel* femmodel); void sealevelrise_core_geometry(FemModel* femmodel); 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); +Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea); +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble oceanarea); +void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks); void sealevelrise_core_viscous(Vector** pU_gia,Vector** pN_gia,FemModel* femmodel,Vector* RSLg); void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg); IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel); Index: ../trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24939) +++ ../trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24940) @@ -95,7 +95,7 @@ int horiz; int geodetic=0; IssmDouble dt; - IssmDouble eartharea,oceanarea; + IssmDouble oceanarea; /*Should we even be here?:*/ femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return; @@ -152,13 +152,13 @@ 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,masks,&eartharea,&oceanarea); + RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&oceanarea); /*call non-eustatic core (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */ - RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,eartharea,oceanarea); + RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,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,masks); + sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,masks); /*compute viscosus (GIA) geodetic signatures:*/ sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg); @@ -353,7 +353,7 @@ }/*}}}*/ -Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/ +Vector* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/ /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ @@ -367,7 +367,7 @@ IssmDouble *longitude = NULL; IssmDouble *radius = NULL; int loop; - IssmDouble eartharea,oceanarea; + IssmDouble oceanarea; /*outputs:*/ IssmDouble eustatic; @@ -387,7 +387,7 @@ RSLgi = new Vector(gsize); /*call the eustatic main module: */ - femmodel->SealevelriseEustatic(RSLgi,&eartharea, &oceanarea,&eustatic, masks, latitude, longitude, radius,loop); //this computes + femmodel->SealevelriseEustatic(RSLgi,&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,masks, oceanarea); @@ -405,11 +405,10 @@ xDelete(radius); /*Assign output pointers and return: */ - *peartharea=eartharea; *poceanarea=oceanarea; return RSLgi; }/*}}}*/ -Vector* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea){ /*{{{*/ +Vector* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector* RSLg_eustatic,IssmDouble oceanarea){ /*{{{*/ /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. non eustatic core of the SLR solution */ @@ -478,7 +477,7 @@ RSLgo = new Vector(gsize); RSLgo->Assemble(); /*call the non eustatic module: */ - femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, eartharea, masks, latitude,longitude, radius,verboseconvolution,loop); + femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, masks, latitude,longitude, radius,verboseconvolution,loop); /*assemble solution vector: */ RSLgo->Assemble(); @@ -487,7 +486,7 @@ /*call rotational feedback module: */ RSLgo_rot = new Vector(gsize); RSLgo_rot->Assemble(); - femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, eartharea, masks, latitude,longitude,radius); + femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks, latitude,longitude,radius); RSLgo_rot->Assemble(); /*save changes in inertia tensor as results: */ @@ -536,7 +535,7 @@ return RSLg; } /*}}}*/ -void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg,IssmDouble eartharea, SealevelMasks* masks){ /*{{{*/ +void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg, SealevelMasks* masks){ /*{{{*/ Vector *U_esa = NULL; Vector *U_north_esa = NULL; @@ -576,7 +575,7 @@ VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices); /*call the elastic main modlule:*/ - femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,eartharea, masks, latitude,longitude,radius,xx,yy,zz,loop,horiz); + femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg, masks, latitude,longitude,radius,xx,yy,zz,loop,horiz); /*Assign output pointers:*/ *pU_esa=U_esa; Index: ../trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 24939) +++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 24940) @@ -208,16 +208,16 @@ #endif #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; - void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; #endif #ifdef _HAVE_SEALEVELRISE_ 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 SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_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!");}; + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_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,int horiz){_error_("not implemented yet!");}; #endif /*}}}*/ Index: ../trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 24939) +++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 24940) @@ -168,15 +168,15 @@ #endif #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; - void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_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, 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!");}; + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_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,int horiz){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; #endif Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24939) +++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24940) @@ -174,15 +174,15 @@ #endif #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast,Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");}; - void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");}; + void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; #endif #ifdef _HAVE_SEALEVELRISE_ 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 SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_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!");}; + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_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,int horiz){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; #endif Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 24939) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 24940) @@ -370,7 +370,7 @@ #endif #ifdef _HAVE_ESA_ virtual void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast, Vector* pX, Vector* pY,IssmDouble* xx,IssmDouble* yy)=0; - virtual void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0; + virtual void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz)=0; #endif #ifdef _HAVE_SEALEVELRISE_ virtual void SetSealevelMasks(SealevelMasks* masks)=0; @@ -377,11 +377,11 @@ virtual IssmDouble GetArea3D(void)=0; virtual IssmDouble GetAreaSpherical(void)=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 SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0; + virtual void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea)=0; virtual void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=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; + virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=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,int horiz)=0; #endif }; Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 24939) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 24940) @@ -159,18 +159,18 @@ #endif #ifdef _HAVE_ESA_ void EsaGeodetic2D(Vector* pUp,Vector* pNorth,Vector* pEast, Vector* pX,Vector* pY,IssmDouble* xx,IssmDouble* yy); - void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea); + void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz); #endif #ifdef _HAVE_SEALEVELRISE_ IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks); - void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble eartharea); + void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks); void SetSealevelMasks(SealevelMasks* masks); void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius); - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); - void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea); - 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); + void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea); + void SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea); + void SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea); + void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius); + void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz); #endif /*}}}*/ /*Tria specific routines:{{{*/ Index: ../trunk-jpl/src/c/classes/FemModel.h =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.h (revision 24939) +++ ../trunk-jpl/src/c/classes/FemModel.h (revision 24940) @@ -162,11 +162,11 @@ void EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseEustatic(Vector* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop); + void SealevelriseEustatic(Vector* pSgi, 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, 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); + void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, 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, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); + void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, 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: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 24939) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 24940) @@ -4634,9 +4634,6 @@ /*}}}*/ void FemModel::EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/ - IssmDouble eartharea=0; - IssmDouble eartharea_cpu=0; - int ns,nsmax; /*Go through elements, and add contribution from each element to the deflection vector wg:*/ @@ -4645,10 +4642,7 @@ /*First, figure out the surface area of Earth: */ for(int i=0;i(elements->GetObjectByOffset(i)); - eartharea_cpu += element->GetAreaSpherical(); } - 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()); /*Figure out max of ns: */ ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); @@ -4658,7 +4652,7 @@ for(int i=0;i(elements->GetObjectByOffset(i)); - element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea); + element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz); } if(i%100==0){ pUp->Assemble(); @@ -4693,7 +4687,7 @@ } /*}}}*/ -void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ +void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ /*serialized vectors:*/ IssmDouble eustatic = 0.; @@ -4701,8 +4695,6 @@ IssmDouble eustatic_cpu_e = 0.; IssmDouble oceanarea = 0.; IssmDouble oceanarea_cpu = 0.; - IssmDouble eartharea = 0.; - IssmDouble eartharea_cpu = 0.; /*Initialize temporary vector that will be used to sum eustatic components * on all local elements, prior to assembly:*/ @@ -4715,7 +4707,6 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); IssmDouble area=element->GetAreaSpherical(); - eartharea_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() ); @@ -4722,13 +4713,10 @@ ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); _assert_(oceanarea>0.); - 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()); - /*Call the sea level rise core: */ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, latitude,longitude,radius,oceanarea,eartharea); + element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, latitude,longitude,radius,oceanarea); eustatic_cpu+=eustatic_cpu_e; } @@ -4746,13 +4734,12 @@ xDelete(RSLgi); /*Assign output pointers:*/ - *peartharea = eartharea; *poceanarea = oceanarea; *peustatic = eustatic; } /*}}}*/ -void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ +void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ /*serialized vectors:*/ IssmDouble* RSLg_old=NULL; @@ -4773,7 +4760,7 @@ /*Call the sea level rise core: */ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius,eartharea); + element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius); } pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL); pRSLgo->Assemble(); @@ -4785,7 +4772,7 @@ } /*}}}*/ -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){/*{{{*/ +void FemModel::SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ /*serialized vectors:*/ IssmDouble* RSLg_old=NULL; @@ -4801,7 +4788,7 @@ IssmDouble moi_list_cpu[3]={0,0,0}; for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks, eartharea); + element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks ); moi_list_cpu[0] += moi_list[0]; moi_list_cpu[1] += moi_list[1]; moi_list_cpu[2] += moi_list[2]; @@ -4861,7 +4848,7 @@ } /*}}}*/ -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){/*{{{*/ +void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ /*serialized vectors:*/ IssmDouble* RSLg=NULL; @@ -4887,7 +4874,7 @@ /*Call the sea level rise core: */ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,eartharea,horiz); + element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,horiz); } pUp->SetValues(gsize,indices,Up,ADD_VAL); Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24939) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24940) @@ -5288,7 +5288,7 @@ return; } /*}}}*/ -void Tria::EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/ +void Tria::EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){ /*{{{*/ /*diverse:*/ int gsize; @@ -5295,7 +5295,7 @@ bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble I; //ice/water loading IssmDouble late,longe,re; IssmDouble lati,longi,ri; @@ -5327,6 +5327,9 @@ /*recover material parameters: */ rho_ice=FindParam(MaterialsRhoIceEnum); rho_earth=FindParam(MaterialsEarthDensityEnum); + + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); /*how many dofs are we working with here? */ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); @@ -5464,7 +5467,7 @@ } /*}}}*/ -void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks, IssmDouble eartharea){/*{{{*/ +void Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/ /*early return if we are not on an ice cap OR ocean:*/ if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){ dI_list[0] = 0.0; // this is important!!! @@ -5474,9 +5477,12 @@ } /*Compute area of element:*/ - IssmDouble area; + IssmDouble area,eartharea; area=GetAreaSpherical(); + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + /*Compute lat,long,radius of elemental centroid: */ bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; @@ -5660,7 +5666,7 @@ return; } /*}}}*/ -void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, 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){ /*{{{*/ /*Computational flags:*/ int bp_compute_fingerprints= 0; @@ -5671,22 +5677,22 @@ 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:*/ - if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); + if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea); } //if(!IsIceInElement()){ /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/ - this->SealevelriseEustaticIce(Sgi,peustatic,masks, latitude,longitude,radius,oceanarea,eartharea); + this->SealevelriseEustaticIce(Sgi,peustatic,masks, latitude,longitude,radius,oceanarea); //} } /*}}}*/ -void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, 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){ /*{{{*/ /*diverse:*/ int gsize,dummy; bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) IssmDouble rho; IssmDouble late,longe,re; @@ -5744,6 +5750,9 @@ this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum); + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + /*recover precomputed green function kernels:*/ DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M); @@ -5835,13 +5844,13 @@ return; } /*}}}*/ -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){ /*{{{*/ /*diverse:*/ int gsize; bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) IssmDouble rho; IssmDouble late,longe,re; @@ -5875,6 +5884,9 @@ rho_water=FindParam(MaterialsRhoSeawaterEnum); rho_earth=FindParam(MaterialsEarthDensityEnum); + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + /*recover love numbers and computational flags: */ this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); @@ -5981,13 +5993,13 @@ return; } /*}}}*/ -void Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){ /*{{{*/ /*diverse:*/ int gsize,dummy; bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble S; //change in water water level(Farrel and Clarke, Equ. 4) IssmDouble late,longe; IssmDouble lati,longi,ri; @@ -6025,6 +6037,9 @@ this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + /*early return if rigid or elastic not requested:*/ if(!computerigid && !computeelastic) return; @@ -6080,7 +6095,7 @@ return; } /*}}}*/ -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){ /*{{{*/ +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,int horiz){ /*{{{*/ /*diverse:*/ int gsize,dummy; @@ -6087,7 +6102,7 @@ bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble area; + IssmDouble area,eartharea; IssmDouble I, S; //change in relative ice thickness and sea level IssmDouble late,longe,re; IssmDouble lati,longi,ri; @@ -6140,6 +6155,9 @@ /*how many dofs are we working with here? */ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + /*recover earth area: */ + this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + /*retrieve indices:*/ if(computerigid){this->inputs2->GetArrayPtr(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);}