Index: ../trunk-jpl/src/m/geometry/planetradius.m =================================================================== --- ../trunk-jpl/src/m/geometry/planetradius.m (revision 25750) +++ ../trunk-jpl/src/m/geometry/planetradius.m (revision 25751) @@ -9,6 +9,8 @@ if strcmpi(planet,'earth'), radius=6.371012*10^6; +elseif strcmpi(planet,'europa'), + radius=1.5008*10^6; else error(['planet type ' planet ' not supported yet!']); end Index: ../trunk-jpl/src/m/classes/solidearth.m =================================================================== --- ../trunk-jpl/src/m/classes/solidearth.m (revision 25750) +++ ../trunk-jpl/src/m/classes/solidearth.m (revision 25751) @@ -13,17 +13,21 @@ planetradius = planetradius('earth'); requested_outputs = {}; transitions = {}; + partitionice = []; + partitionhydro = []; end methods function self = solidearth(varargin) % {{{ switch nargin case 0 - self=setdefaultparameters(self); + self=setdefaultparameters(self,earth); + case 1 + self=setdefaultparameters(self,varargin{:}); otherwise - error('constructor not supported'); + error('solidearth constructor error message: zero or one argument only!'); end end % }}} - function self = setdefaultparameters(self) % {{{ + function self = setdefaultparameters(self,planet) % {{{ %output default: self.requested_outputs={'default'}; @@ -31,8 +35,12 @@ %transitions should be a cell array of vectors: self.transitions={}; + %no partitions requested for barystatic contribution: + self.partitionice=[]; + self.partitionhydro=[]; + %earth radius - self.planetradius= planetradius('earth'); + self.planetradius= planetradius(planet); end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -61,6 +69,8 @@ fielddisplay(self,'planetradius','planet radius [m]'); fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'); fielddisplay(self,'requested_outputs','additional outputs requested'); + fielddisplay(self,'partitionice','ice partition vector for barystatic contribution'); + fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution'); self.settings.disp(); self.surfaceload.disp(); self.lovenumbers.disp(); @@ -72,7 +82,24 @@ WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','planetradius','format','Double'); WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray'); + + if ~isempty(self.partitionice), + npartice=max(self.partitionice)+2; + else + npartice=0; + end + if ~isempty(self.partitionhydro), + nparthydro=max(self.partitionhydro)+2; + else + nparthydro=0; + end + + WriteData(fid,prefix,'object',self,'fieldname','partitionice','mattype',1,'format','DoubleMat'); + WriteData(fid,prefix,'data',npartice,'format','Integer','name','md.solidearth.npartice'); + WriteData(fid,prefix,'object',self,'fieldname','partitionhydro','mattype',1,'format','DoubleMat'); + WriteData(fid,prefix,'data',nparthydro,'format','Integer','name','md.solidearth.nparthydro'); + self.settings.marshall(prefix,md,fid); self.surfaceload.marshall(prefix,md,fid); self.lovenumbers.marshall(prefix,md,fid); @@ -97,6 +124,7 @@ self.rotational.savemodeljs(fid,modelname); writejscellstring(fid,[modelname '.solidearth.requested_outputs'],self.requested_outputs); writejscellarray(fid,[modelname '.solidearth.transitions'],self.transitions); + writejscellarray(fid,[modelname '.solidearth.partition'],self.partition); end % }}} function self = extrude(self,md) % {{{ self.sealevel=project3d(md,'vector',self.sealevel,'type','node'); Index: ../trunk-jpl/src/m/classes/model.m =================================================================== --- ../trunk-jpl/src/m/classes/model.m (revision 25750) +++ ../trunk-jpl/src/m/classes/model.m (revision 25751) @@ -191,15 +191,15 @@ end methods function md = model(varargin) % {{{ + switch nargin case 0 - md=setdefaultparameters(md); - case 1 - error('model constructor not supported yet'); - + md=setdefaultparameters(md,'earth'); otherwise - error('model constructor error message: 0 of 1 argument only in input.'); + options=pairoptions(varargin{:}); + planet=getfieldvalue(options,'planet','earth'); + md=setdefaultparameters(md,planet); end end %}}} @@ -1286,7 +1286,7 @@ md.solidearth.rotational.angularvelocity=structmd.slr.angular_velocity; end end% }}} - function md = setdefaultparameters(md) % {{{ + function md = setdefaultparameters(md,planet) % {{{ %initialize subclasses md.mesh = mesh2d(); @@ -1298,7 +1298,7 @@ md.basalforcings = basalforcings(); md.friction = friction(); md.rifts = rifts(); - md.solidearth = solidearth(); + md.solidearth = solidearth(planet); md.dsl = dsl(); md.timestepping = timestepping(); md.groundingline = groundingline(); Index: ../trunk-jpl/src/m/classes/solidearthsettings.m =================================================================== --- ../trunk-jpl/src/m/classes/solidearthsettings.m (revision 25750) +++ ../trunk-jpl/src/m/classes/solidearthsettings.m (revision 25751) @@ -63,6 +63,11 @@ md = checkfield(md,'fieldname','solidearth.settings.degacc','size',[1 1],'>=',1e-10); md = checkfield(md,'fieldname','solidearth.settings.horiz','NaN',1,'Inf',1,'values',[0 1]); + %checks on computational flags: + if self.elastic==1 & self.rigid==0, + error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set'); + end + %a coupler to a planet model is provided. if self.computesealevelchange, if md.transient.iscoupler, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25750) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25751) @@ -115,6 +115,8 @@ case CumBslrEnum : return "CumBslr"; case CumBslrIceEnum : return "CumBslrIce"; case CumBslrHydroEnum : return "CumBslrHydro"; + case CumBslrIcePartitionEnum : return "CumBslrIcePartition"; + case CumBslrHydroPartitionEnum : return "CumBslrHydroPartition"; case CumGmtslrEnum : return "CumGmtslr"; case CumGmslrEnum : return "CumGmslr"; case DamageC1Enum : return "DamageC1"; @@ -329,6 +331,10 @@ case RootPathEnum : return "RootPath"; case ModelnameEnum : return "Modelname"; case SaveResultsEnum : return "SaveResults"; + case SolidearthPartitionIceEnum : return "SolidearthPartitionIce"; + case SolidearthPartitionHydroEnum : return "SolidearthPartitionHydro"; + case SolidearthNpartIceEnum : return "SolidearthNpartIce"; + case SolidearthNpartHydroEnum : return "SolidearthNpartHydro"; case SolidearthPlanetRadiusEnum : return "SolidearthPlanetRadius"; case SolidearthPlanetAreaEnum : return "SolidearthPlanetArea"; case SolidearthSettingsAbstolEnum : return "SolidearthSettingsAbstol"; @@ -346,6 +352,7 @@ case SealevelriseGRigidEnum : return "SealevelriseGRigid"; case SealevelriseGElasticEnum : return "SealevelriseGElastic"; case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange"; + case SolidearthSettingsGlfractionEnum : return "SolidearthSettingsGlfraction"; case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency"; case SealevelriseHElasticEnum : return "SealevelriseHElastic"; case SolidearthSettingsHorizEnum : return "SolidearthSettingsHoriz"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25750) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25751) @@ -115,6 +115,8 @@ else if (strcmp(name,"CumBslr")==0) return CumBslrEnum; else if (strcmp(name,"CumBslrIce")==0) return CumBslrIceEnum; else if (strcmp(name,"CumBslrHydro")==0) return CumBslrHydroEnum; + else if (strcmp(name,"CumBslrIcePartition")==0) return CumBslrIcePartitionEnum; + else if (strcmp(name,"CumBslrHydroPartition")==0) return CumBslrHydroPartitionEnum; else if (strcmp(name,"CumGmtslr")==0) return CumGmtslrEnum; else if (strcmp(name,"CumGmslr")==0) return CumGmslrEnum; else if (strcmp(name,"DamageC1")==0) return DamageC1Enum; @@ -134,12 +136,12 @@ else if (strcmp(name,"DamageStressUBound")==0) return DamageStressUBoundEnum; else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum; else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum; - else if (strcmp(name,"DomainType")==0) return DomainTypeEnum; - else if (strcmp(name,"DslModel")==0) return DslModelEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"DslModelid")==0) return DslModelidEnum; + if (strcmp(name,"DomainType")==0) return DomainTypeEnum; + else if (strcmp(name,"DslModel")==0) return DslModelEnum; + else if (strcmp(name,"DslModelid")==0) return DslModelidEnum; else if (strcmp(name,"DslNummodels")==0) return DslNummodelsEnum; else if (strcmp(name,"DslComputeFingerprints")==0) return DslComputeFingerprintsEnum; else if (strcmp(name,"EarthId")==0) return EarthIdEnum; @@ -257,12 +259,12 @@ else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum; else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum; else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; - else if (strcmp(name,"LoveR0")==0) return LoveR0Enum; - else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum; + if (strcmp(name,"LoveR0")==0) return LoveR0Enum; + else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; + else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum; else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum; else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum; else if (strcmp(name,"MasstransportHydrostaticAdjustment")==0) return MasstransportHydrostaticAdjustmentEnum; @@ -335,6 +337,10 @@ else if (strcmp(name,"RootPath")==0) return RootPathEnum; else if (strcmp(name,"Modelname")==0) return ModelnameEnum; else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum; + else if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum; + else if (strcmp(name,"SolidearthPartitionHydro")==0) return SolidearthPartitionHydroEnum; + else if (strcmp(name,"SolidearthNpartIce")==0) return SolidearthNpartIceEnum; + else if (strcmp(name,"SolidearthNpartHydro")==0) return SolidearthNpartHydroEnum; else if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum; else if (strcmp(name,"SolidearthPlanetArea")==0) return SolidearthPlanetAreaEnum; else if (strcmp(name,"SolidearthSettingsAbstol")==0) return SolidearthSettingsAbstolEnum; @@ -352,6 +358,7 @@ else if (strcmp(name,"SealevelriseGRigid")==0) return SealevelriseGRigidEnum; else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum; else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum; + else if (strcmp(name,"SolidearthSettingsGlfraction")==0) return SolidearthSettingsGlfractionEnum; else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum; else if (strcmp(name,"SealevelriseHElastic")==0) return SealevelriseHElasticEnum; else if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum; @@ -375,7 +382,10 @@ else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum; else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum; else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; - else if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum; else if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum; else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum; else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum; @@ -382,10 +392,7 @@ else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum; else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum; else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; + else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum; else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum; else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum; @@ -498,7 +505,10 @@ else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; else if (strcmp(name,"InputsSTART")==0) return InputsSTARTEnum; else if (strcmp(name,"Adjoint")==0) return AdjointEnum; - else if (strcmp(name,"Adjointp")==0) return AdjointpEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"Adjointp")==0) return AdjointpEnum; else if (strcmp(name,"Adjointx")==0) return AdjointxEnum; else if (strcmp(name,"Adjointy")==0) return AdjointyEnum; else if (strcmp(name,"Adjointz")==0) return AdjointzEnum; @@ -505,10 +515,7 @@ else if (strcmp(name,"Air")==0) return AirEnum; else if (strcmp(name,"Approximation")==0) return ApproximationEnum; else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum; + else if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum; else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum; else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum; else if (strcmp(name,"BalancethicknessThickeningRate")==0) return BalancethicknessThickeningRateEnum; @@ -621,7 +628,10 @@ else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; - else if (strcmp(name,"Frictionf")==0) return FrictionfEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"Frictionf")==0) return FrictionfEnum; else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum; else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum; else if (strcmp(name,"FrontalForcingsThermalForcing")==0) return FrontalForcingsThermalForcingEnum; @@ -628,10 +638,7 @@ else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum; else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum; else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"NGia")==0) return NGiaEnum; + else if (strcmp(name,"NGia")==0) return NGiaEnum; else if (strcmp(name,"NGiaRate")==0) return NGiaRateEnum; else if (strcmp(name,"UGia")==0) return UGiaEnum; else if (strcmp(name,"UGiaRate")==0) return UGiaRateEnum; @@ -744,7 +751,10 @@ else if (strcmp(name,"SealevelriseGE")==0) return SealevelriseGEEnum; else if (strcmp(name,"SealevelriseGN")==0) return SealevelriseGNEnum; else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; - else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; else if (strcmp(name,"SedimentHeadSubstep")==0) return SedimentHeadSubstepEnum; else if (strcmp(name,"SedimentHeadTransient")==0) return SedimentHeadTransientEnum; else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum; @@ -751,10 +761,7 @@ else if (strcmp(name,"SedimentHeadStacked")==0) return SedimentHeadStackedEnum; else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum; else if (strcmp(name,"SigmaVM")==0) return SigmaVMEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"SmbA")==0) return SmbAEnum; + else if (strcmp(name,"SmbA")==0) return SmbAEnum; else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum; else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum; else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; @@ -867,7 +874,10 @@ else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum; else if (strcmp(name,"Area")==0) return AreaEnum; - else if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"SealevelArea")==0) return SealevelAreaEnum; else if (strcmp(name,"SurfaceArea")==0) return SurfaceAreaEnum; else if (strcmp(name,"SurfaceAverageVelMisfit")==0) return SurfaceAverageVelMisfitEnum; else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum; @@ -874,10 +884,7 @@ else if (strcmp(name,"Surface")==0) return SurfaceEnum; else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum; else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; + else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum; else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum; @@ -990,7 +997,10 @@ else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum; else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum; else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum; - else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum; else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum; else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum; else if (strcmp(name,"Outputdefinition80")==0) return Outputdefinition80Enum; @@ -997,10 +1007,7 @@ else if (strcmp(name,"Outputdefinition81")==0) return Outputdefinition81Enum; else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum; else if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; + else if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum; else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum; else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum; @@ -1113,7 +1120,10 @@ else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; else if (strcmp(name,"FemModel")==0) return FemModelEnum; else if (strcmp(name,"FileParam")==0) return FileParamEnum; - else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum; else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum; else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum; else if (strcmp(name,"FloatingMeltRate")==0) return FloatingMeltRateEnum; @@ -1120,10 +1130,7 @@ else if (strcmp(name,"Free")==0) return FreeEnum; else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; + else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum; else if (strcmp(name,"Fset")==0) return FsetEnum; else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum; @@ -1236,7 +1243,10 @@ else if (strcmp(name,"MeshX")==0) return MeshXEnum; else if (strcmp(name,"MeshY")==0) return MeshYEnum; else if (strcmp(name,"MinVel")==0) return MinVelEnum; - else if (strcmp(name,"MinVx")==0) return MinVxEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"MinVx")==0) return MinVxEnum; else if (strcmp(name,"MinVy")==0) return MinVyEnum; else if (strcmp(name,"MinVz")==0) return MinVzEnum; else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum; @@ -1243,10 +1253,7 @@ else if (strcmp(name,"Moulin")==0) return MoulinEnum; else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; else if (strcmp(name,"Mpi")==0) return MpiEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; + else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; else if (strcmp(name,"Mumps")==0) return MumpsEnum; else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum; else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum; @@ -1359,7 +1366,10 @@ else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum; else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum; else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum; - else if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum; else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum; else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum; else if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum; @@ -1366,10 +1376,7 @@ else if (strcmp(name,"TotalGroundedBmbScaled")==0) return TotalGroundedBmbScaledEnum; else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; + else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; Index: ../trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25750) +++ ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25751) @@ -113,6 +113,8 @@ syn keyword cConstant CumBslrEnum syn keyword cConstant CumBslrIceEnum syn keyword cConstant CumBslrHydroEnum +syn keyword cConstant CumBslrIcePartitionEnum +syn keyword cConstant CumBslrHydroPartitionEnum syn keyword cConstant CumGmtslrEnum syn keyword cConstant CumGmslrEnum syn keyword cConstant DamageC1Enum @@ -327,6 +329,10 @@ syn keyword cConstant RootPathEnum syn keyword cConstant ModelnameEnum syn keyword cConstant SaveResultsEnum +syn keyword cConstant SolidearthPartitionIceEnum +syn keyword cConstant SolidearthPartitionHydroEnum +syn keyword cConstant SolidearthNpartIceEnum +syn keyword cConstant SolidearthNpartHydroEnum syn keyword cConstant SolidearthPlanetRadiusEnum syn keyword cConstant SolidearthPlanetAreaEnum syn keyword cConstant SolidearthSettingsAbstolEnum @@ -344,6 +350,7 @@ syn keyword cConstant SealevelriseGRigidEnum syn keyword cConstant SealevelriseGElasticEnum syn keyword cConstant SolidearthSettingsComputesealevelchangeEnum +syn keyword cConstant SolidearthSettingsGlfractionEnum syn keyword cConstant SolidearthSettingsRunFrequencyEnum syn keyword cConstant SealevelriseHElasticEnum syn keyword cConstant SolidearthSettingsHorizEnum Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25750) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25751) @@ -107,6 +107,8 @@ CumBslrEnum, CumBslrIceEnum, CumBslrHydroEnum, + CumBslrIcePartitionEnum, + CumBslrHydroPartitionEnum, CumGmtslrEnum, CumGmslrEnum, DamageC1Enum, @@ -321,6 +323,10 @@ RootPathEnum, ModelnameEnum, SaveResultsEnum, + SolidearthPartitionIceEnum, + SolidearthPartitionHydroEnum, + SolidearthNpartIceEnum, + SolidearthNpartHydroEnum, SolidearthPlanetRadiusEnum, SolidearthPlanetAreaEnum, SolidearthSettingsAbstolEnum, @@ -338,6 +344,7 @@ SealevelriseGRigidEnum, SealevelriseGElasticEnum, SolidearthSettingsComputesealevelchangeEnum, + SolidearthSettingsGlfractionEnum, SolidearthSettingsRunFrequencyEnum, SealevelriseHElasticEnum, SolidearthSettingsHorizEnum, Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25750) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25751) @@ -5462,160 +5462,121 @@ area=GetAreaSpherical(); this->inputs->SetDoubleInput(AreaEnum,this->lid,area); this->AddInput(SealevelAreaEnum,&area,P0Enum); + if(!computerigid)return; /*recover precomputed green function kernels:*/ - if(computerigid){ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M); - /*allocate indices:*/ - indices=xNew(gsize); - G=xNewZeroInit(gsize); + /*allocate indices:*/ + indices=xNew(gsize); + G=xNewZeroInit(gsize); - if(computeelastic){ - parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); + if(computeelastic){ + parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); - parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); + parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); - parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); + parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); - /*allocate indices:*/ - GU=xNewZeroInit(gsize); - if(horiz){ - GN=xNewZeroInit(gsize); - GE=xNewZeroInit(gsize); - } + /*allocate indices:*/ + GU=xNewZeroInit(gsize); + if(horiz){ + GN=xNewZeroInit(gsize); + GE=xNewZeroInit(gsize); + } + } - /* Where is the centroid of this element:*/ + /* Where is the centroid of this element:{{{*/ - /*retrieve coordinates: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + /*retrieve coordinates: */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - /*figure out gravity center of our element (Cartesian): */ - 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; + /*figure out gravity center of our element (Cartesian): */ + 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; - /*compute gravity center in lat,long: */ - late= asin(z_element/planetradius); - longe = atan2(y_element,x_element); + /*compute gravity center in lat,long: */ + late= asin(z_element/planetradius); + longe = atan2(y_element,x_element); + /*}}}*/ - constant=3/rho_earth*area/planetarea; + 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); + index=reCast(indices[i]); - /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */ - lati=latitude[i]/180*PI; longi=longitude[i]/180*PI; - delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>PI)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); - index=reCast(indices[i]); + /*Rigid earth gravitational perturbation: */ + G[i]+=G_rigid_precomputed[index]; + if(elastic) G[i]+=G_elastic_precomputed[index]; + G[i]=G[i]*constant; - /*Rigid earth gravitational perturbation: */ - G[i]+=G_rigid_precomputed[index]; - G[i]+=G_elastic_precomputed[index]; - G[i]=G[i]*constant; - - /*Elastic components:*/ - GU[i] = constant * U_elastic_precomputed[index]; - if(horiz){ - /*Compute azimuths, both north and east components: */ - x = xx[i]; y = yy[i]; z = zz[i]; - if(latitude[i]==90){ - x=1e-12; y=1e-12; - } - if(latitude[i]==-90){ - x=1e-12; y=1e-12; - } - dx = x_element-x; dy = y_element-y; dz = z_element-z; - N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - - GN[i] = constant*H_elastic_precomputed[index]*N_azim; - GE[i] = constant*H_elastic_precomputed[index]*E_azim; + /*Elastic components:*/ + if(computeelastic){ + GU[i] = constant * U_elastic_precomputed[index]; + if(horiz){ + /*Compute azimuths, both north and east components: */ + x = xx[i]; y = yy[i]; z = zz[i]; + if(latitude[i]==90){ + x=1e-12; y=1e-12; } - } + if(latitude[i]==-90){ + x=1e-12; y=1e-12; + } + dx = x_element-x; dy = y_element-y; dz = z_element-z; + N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - /*Add in inputs:*/ - this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); - this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); - this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); - if(horiz){ - this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize); - this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); + GN[i] = constant*H_elastic_precomputed[index]*N_azim; + GE[i] = constant*H_elastic_precomputed[index]*E_azim; } - this->inputs->SetDoubleInput(AreaEnum,this->lid,area); - this->AddInput(SealevelAreaEnum,&area,P0Enum); + } + } - /*Free allocations:*/ - #ifdef _HAVE_RESTRICT_ - delete GU; - if(horiz){ - delete GN; - delete GE; - } - #else - xDelete(GU); - if(horiz){ - xDelete(GN); - xDelete(GE); - } - #endif - } else { /*computeelastic==false*/ - /* Where is the centroid of this element:*/ + /*Add in inputs:*/ + this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); + this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); + if(computeelastic){ + this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); + if(horiz){ + this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize); + this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); + } + } + this->inputs->SetDoubleInput(AreaEnum,this->lid,area); + this->AddInput(SealevelAreaEnum,&area,P0Enum); - /*retrieve coordinates: */ - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - - /*figure out gravity center of our element (Cartesian): */ - 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; - - /*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); - index=reCast(indices[i]); - - /*Rigid earth gravitational perturbation: */ - G[i]+=G_rigid_precomputed[index]; - G[i]=G[i]*constant; - } - - /*Add in inputs:*/ - this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); - this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); - this->inputs->SetDoubleInput(AreaEnum,this->lid,area); - this->AddInput(SealevelAreaEnum,&area,P0Enum); + /*Free allocations:*/ + #ifdef _HAVE_RESTRICT_ + delete G; + if(computeelastic){ + delete GU; + if(horiz){ + delete GN; + delete GE; } - - /*Free allocations:*/ - #ifdef _HAVE_RESTRICT_ - delete indices; - delete G; - #else - xDelete(indices); - xDelete(G); - #endif - } else { /*computerigid==false*/ - // do nothing } + #else + xDelete(G); + if(elastic){ + xDelete(GU); + if(horiz){ + xDelete(GN); + xDelete(GE); + } + } + #endif return; } @@ -5630,6 +5591,7 @@ bool notfullygrounded=false; bool scaleoceanarea= false; bool computerigid= false; + int glfraction=1; /*output: */ IssmDouble bslrice=0; @@ -5689,6 +5651,10 @@ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias + if(glfraction==0)phi=1; + #ifdef _ISSM_DEBUG_ + this->AddInput2(SealevelEustaticMaskEnum,&phi,P0Enum); + #endif } else phi=1.0; @@ -5736,6 +5702,11 @@ for(int i=0;iSid()]+1; + barystatic_contribution->SetValue(id,bslrice,ADD_VAL); + } /*return :*/ return bslrice; } @@ -5761,11 +5732,12 @@ /*output:*/ IssmDouble bslrhydro = 0; - /*early return if we are on an ice cap:*/ - if(masks->isiceonly[this->lid]){ + /*early return if we are on an ice cap. Nop, we may have hydro + * and ice on the same masscon:*/ + /*if(masks->isiceonly[this->lid]){ bslrhydro=0; return bslrhydro; - } + }*/ /*early return if we are fully floating:*/ if (masks->isfullyfloating[this->lid]){ @@ -5806,6 +5778,11 @@ for(int i=0;iSid()]+1; + barystatic_contribution->SetValue(id,bslrhydro,ADD_VAL); + } /*output:*/ return bslrhydro; } Index: ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 25750) +++ ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp (revision 25751) @@ -176,6 +176,7 @@ IssmDouble planetarea=0; bool rigid=false; bool elastic=false; + bool rotation=false; int numoutputs; char** requestedoutputs = NULL; @@ -202,6 +203,7 @@ parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.ocean_area_scaling",SolidearthSettingsOceanAreaScalingEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.glfraction",SolidearthSettingsGlfractionEnum)); parameters->AddObject(new DoubleParam(CumBslrEnum,0.0)); parameters->AddObject(new DoubleParam(CumBslrIceEnum,0.0)); parameters->AddObject(new DoubleParam(CumBslrHydroEnum,0.0)); @@ -212,6 +214,29 @@ planetarea=4*PI*planetradius*planetradius; parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea)); + /*Deal with partition of the barystatic contribution:*/ + iomodel->FetchData(&npartice,"md.solidearth.npartice"); + parameters->AddObject(new IntParam(SolidearthNpartIceEnum,npartice)); + if(npartice){ + iomodel->FetchData(&partitionice,&nel,NULL,"md.solidearth.partitionice"); + parameters->AddObject(new DoubleMatParam(SolidearthPartitionIceEnum,partitionice,nel,1)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.npartice",SolidearthNpartIceEnum)); + bslrice_partition=xNewZeroInit(npartice); + parameters->AddObject(new DoubleMatParam(CumBslrIcePartitionEnum,bslrice_partition,npartice,1)); + xDelete(partitionice); + } + iomodel->FetchData(&nparthydro,"md.solidearth.nparthydro"); + parameters->AddObject(new IntParam(SolidearthNpartHydroEnum,nparthydro)); + if(nparthydro){ + iomodel->FetchData(&partitionhydro,&nel,NULL,"md.solidearth.partitionhydro"); + parameters->AddObject(new DoubleMatParam(SolidearthPartitionHydroEnum,partitionhydro,nel,1)); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.nparthydro",SolidearthNpartHydroEnum)); + bslrhydro_partition=xNewZeroInit(nparthydro); + parameters->AddObject(new DoubleMatParam(CumBslrHydroPartitionEnum,bslrhydro_partition,nparthydro,1)); + xDelete(partitionhydro); + } + + /*Deal with dsl multi-model ensembles: {{{*/ iomodel->FetchData(&dslmodel,"md.dsl.model"); parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum)); @@ -233,136 +258,170 @@ /*Deal with elasticity {{{*/ iomodel->FetchData(&rigid,"md.solidearth.settings.rigid"); iomodel->FetchData(&elastic,"md.solidearth.settings.elastic"); - if (rigid) { - if (elastic) { - /*love numbers: */ - iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); - iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); - iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l"); - iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th"); - iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); - iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); + iomodel->FetchData(&rotation,"md.solidearth.settings.rotation"); - parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); - parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1)); - parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); + if(elastic | rigid){ + /*compute green functions for a range of angles*/ + iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); + M=reCast(180./degacc+1.); + } - /*compute elastic green function for a range of angles*/ - iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); - M=reCast(180./degacc+1.); + /*love numbers: */ + if(elastic){ + iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); + iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); + iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l"); + iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th"); + iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); + iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_rigid=xNew(M,"t"); - G_elastic=xNew(M,"t"); - U_elastic=xNew(M,"t"); - H_elastic=xNew(M,"t"); - #else - G_rigid=xNew(M); - G_elastic=xNew(M); - U_elastic=xNew(M); - H_elastic=xNew(M); - #endif + parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); + parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1)); + parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); - /*compute combined legendre + love number (elastic green function:*/ - m=DetermineLocalSize(M,IssmComm::GetComm()); - GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_elastic_local=xNew(m,"t"); - G_rigid_local=xNew(m,"t"); - U_elastic_local=xNew(m,"t"); - H_elastic_local=xNew(m,"t"); - #else - G_elastic_local=xNew(m); - G_rigid_local=xNew(m); - U_elastic_local=xNew(m); - H_elastic_local=xNew(m); - #endif + // AD performance is sensitive to calls to ensurecontiguous. + // // Providing "t" will cause ensurecontiguous to be called. + #ifdef _HAVE_AD_ + G_elastic=xNew(M,"t"); + U_elastic=xNew(M,"t"); + H_elastic=xNew(M,"t"); + #else + G_elastic=xNew(M); + U_elastic=xNew(M); + H_elastic=xNew(M); + #endif + } + if(rigid){ + #ifdef _HAVE_AD_ + G_rigid=xNew(M,"t"); + #else + G_rigid=xNew(M); + #endif + } + + if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); - for(int i=lower_row;i(i)*degacc * PI / 180.0; + if(rigid | elastic){ - G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); - G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; - U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; - H_elastic_local[i-lower_row]= 0; - IssmDouble Pn = 0.; - IssmDouble Pn1 = 0.; - IssmDouble Pn2 = 0.; - IssmDouble Pn_p = 0.; - IssmDouble Pn_p1 = 0.; - IssmDouble Pn_p2 = 0.; + /*compute combined legendre + love number (elastic green function:*/ + m=DetermineLocalSize(M,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); + } + if(elastic){ + #ifdef _HAVE_AD_ + G_elastic_local=xNew(m,"t"); + U_elastic_local=xNew(m,"t"); + H_elastic_local=xNew(m,"t"); + #else + G_elastic_local=xNew(m); + U_elastic_local=xNew(m); + H_elastic_local=xNew(m); + #endif + } + if(rigid){ + #ifdef _HAVE_AD_ + G_rigid_local=xNew(m,"t"); + #else + G_rigid_local=xNew(m); + #endif + } - for (int n=0;n(i)*degacc * PI / 180.0; + G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); + } + } + if(elastic){ + for(int i=lower_row;i(i)*degacc * PI / 180.0; - deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]); - deltalove_U = (love_h[n]-love_h[nl-1]); + G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; + U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; + H_elastic_local[i-lower_row]= 0; + IssmDouble Pn = 0.; + IssmDouble Pn1 = 0.; + IssmDouble Pn2 = 0.; + IssmDouble Pn_p = 0.; + IssmDouble Pn_p1 = 0.; + IssmDouble Pn_p2 = 0.; - /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */ - if(n==0){ - Pn=1; - Pn_p=0; - } - else if(n==1){ - Pn = cos(alpha); - Pn_p = 1; - } - else{ - Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n; - Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n; - } - Pn2=Pn1; Pn1=Pn; - Pn_p2=Pn_p1; Pn_p1=Pn_p; + for (int n=0;n(IssmComm::GetSize()); - int* displs=xNew(IssmComm::GetSize()); + /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/ + int* recvcounts=xNew(IssmComm::GetSize()); + int* displs=xNew(IssmComm::GetSize()); - //recvcounts: - ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); + //recvcounts: + ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); - /*displs: */ - ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); + /*displs: */ + ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); - /*All gather:*/ - ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + /*All gather:*/ + ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + if(elastic){ ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - - /*free resources: */ - xDelete(recvcounts); - xDelete(displs); + } + + /*free resources: */ + xDelete(recvcounts); + xDelete(displs); - /*}}}*/ + /*Avoid singularity at 0: */ + G_rigid[0]=G_rigid[1]; + G_elastic[0]=G_elastic[1]; + U_elastic[0]=U_elastic[1]; + H_elastic[0]=H_elastic[1]; - /*Avoid singularity at 0: */ - G_rigid[0]=G_rigid[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); - G_elastic[0]=G_elastic[1]; + /*Save our precomputed tables into parameters*/ + parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); + if(elastic){ parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M)); - U_elastic[0]=U_elastic[1]; parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M)); - H_elastic[0]=H_elastic[1]; parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M)); + } - /*free resources: */ + /*free resources: */ + xDelete(G_rigid); + xDelete(G_rigid_local); + if(elastic){ xDelete(love_h); xDelete(love_k); xDelete(love_l); @@ -369,8 +428,6 @@ xDelete(love_th); xDelete(love_tk); xDelete(love_tl); - xDelete(G_rigid); - xDelete(G_rigid_local); xDelete(G_elastic); xDelete(G_elastic_local); xDelete(U_elastic); @@ -377,69 +434,10 @@ xDelete(U_elastic_local); xDelete(H_elastic); xDelete(H_elastic_local); - } else { /*elastic==false*/ - /*compute elastic green function for a range of angles*/ - iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); - M=reCast(180./degacc+1.); - - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_rigid=xNew(M,"t"); - #else - G_rigid=xNew(M); - #endif - - /*compute combined legendre + love number (elastic green function:*/ - m=DetermineLocalSize(M,IssmComm::GetComm()); - GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_rigid_local=xNew(m,"t"); - #else - G_rigid_local=xNew(m); - #endif - - for(int i=lower_row;i(i)*degacc*PI/180.0; - G_rigid_local[i-lower_row]=.5/sin(alpha/2.0); - } - - /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/ - int* recvcounts=xNew(IssmComm::GetSize()); - int* displs=xNew(IssmComm::GetSize()); - - //recvcounts: - ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*displs: */ - ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*All gather:*/ - ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - - /*free resources: */ - xDelete(recvcounts); - xDelete(displs); - - /*}}}*/ - - /*Avoid singularity at 0: */ - G_rigid[0]=G_rigid[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); - - /*free resources: */ - xDelete(G_rigid); - xDelete(G_rigid_local); } - } else { /*rigid==false*/ - if (elastic) { - _error_("Must set md.solidearth.settings.rigid=1 when setting md.solidearth.settings.elastic=1"); - } } - + /*}}}*/ + /*Indicate we have not yet run the Geometry Core module: */ parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false));