Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27284) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27285) @@ -256,4 +256,7 @@ iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.malag_coefs"); parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAmalagcoefsEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.monthly_effects"); + parameters->AddObject(new DoubleMatParam(ThermalForcingMonthlyEffectsEnum,transparam,M,N)); xDelete(transparam); /*Do not break here, generic FrontalForcingsRignot parameters still to be retrieved*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27284) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27285) @@ -2775,4 +2775,52 @@ this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); +}/*}}}*/ +void Element::MonthlyEffectBasin(IssmDouble* monthlyeff, int enum_type){/*{{{*/ + + /*Variable declaration*/ + const int numvertices = this->GetNumberOfVertices(); + int basinid,mindex,basinenum_type,varenum_type; + IssmDouble time,fracyear,yts; + IssmDouble monthsteps[12] = {0.,1./12,2./12,3./12,4./12,5./12,6./12,7./12,8./12,9./12,10./12,11./12}; + IssmDouble* monthlyeff_b = xNew(12); + IssmDouble* varlist = xNew(numvertices); + + /*Get field-specific enums*/ + switch(enum_type){ + case(FrontalForcingsRignotarmaEnum): + basinenum_type = FrontalForcingsBasinIdEnum; + varenum_type = ThermalForcingEnum; + break; + } + + /*Evaluate the month index*/ + this->parameters->FindParam(&yts,ConstantsYtsEnum); + this->parameters->FindParam(&time,TimeEnum); + fracyear = time/yts-floor(time/yts); + for(int i=1;i<12;i++){ + if(fracyear>=monthsteps[i-1]) mindex = i-1; + } + if(fracyear>=monthsteps[11]) mindex = 11; + + /*Get basin-specific parameters of the element*/ + this->GetInputValue(&basinid,basinenum_type); + for(int ii=0;ii<12;ii++){ + monthlyeff_b[ii] = monthlyeff[basinid*12+ii]; + } + + /*Retrieve values non-adjusted for monthly effects*/ + this->GetInputListOnVertices(varlist,varenum_type); + + /*Adjust values using monthly effects*/ + for(int v=0;vAddInput(varenum_type,varlist,P1Enum); + + /*Clean-up*/ + xDelete(monthlyeff_b); + xDelete(varlist); }/*}}}*/ void Element::BeckmannGoosseFloatingiceMeltingRate(){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27284) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27285) @@ -164,4 +164,5 @@ void MigrateGroundingLine(IssmDouble* sheet_ungrounding); void MismipFloatingiceMeltingRate(); + void MonthlyEffectBasin(IssmDouble* monthlyeff, int enum_type); void BeckmannGoosseFloatingiceMeltingRate(); void MungsmtpParameterization(void); Index: /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 27284) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 27285) @@ -59,4 +59,5 @@ IssmDouble* arlagcoefs = NULL; IssmDouble* malagcoefs = NULL; + IssmDouble* monthlyeff = NULL; femmodel->parameters->FindParam(&tinit_arma,FrontalForcingsARMAInitialTimeEnum); @@ -65,4 +66,5 @@ femmodel->parameters->FindParam(&arlagcoefs,&M,&Narlagcoefs,FrontalForcingsARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(Narlagcoefs==arorder); femmodel->parameters->FindParam(&malagcoefs,&M,&Nmalagcoefs,FrontalForcingsARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(Nmalagcoefs==maorder); + femmodel->parameters->FindParam(&monthlyeff,&M,&N,ThermalForcingMonthlyEffectsEnum); _assert_(M==numbasins); _assert_(N==12); femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); @@ -83,6 +85,9 @@ for(Object* &object:femmodel->elements->objects){ Element* element = xDynamicCast(object); + /*Compute ARMA*/ element->ArmaProcess(isstepforarma,arorder,maorder,telapsed_arma,tstep_arma,termconstant,trend,arlagcoefs,malagcoefs,istfstochastic,FrontalForcingsRignotarmaEnum); - } + /*Compute monthly effects*/ + element->MonthlyEffectBasin(monthlyeff,FrontalForcingsRignotarmaEnum); + } /*Cleanup*/ @@ -91,3 +96,4 @@ xDelete(arlagcoefs); xDelete(malagcoefs); + xDelete(monthlyeff); }/*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27284) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27285) @@ -564,4 +564,5 @@ syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum syn keyword cConstant StressbalanceShelfDampeningEnum +syn keyword cConstant ThermalForcingMonthlyEffectsEnum syn keyword cConstant ThermalIsdrainicecolumnEnum syn keyword cConstant ThermalIsdynamicbasalspcEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27284) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27285) @@ -558,4 +558,5 @@ StressbalanceRiftPenaltyThresholdEnum, StressbalanceShelfDampeningEnum, + ThermalForcingMonthlyEffectsEnum, ThermalIsdrainicecolumnEnum, ThermalIsdynamicbasalspcEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27284) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27285) @@ -566,4 +566,5 @@ case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold"; case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening"; + case ThermalForcingMonthlyEffectsEnum : return "ThermalForcingMonthlyEffects"; case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn"; case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27284) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27285) @@ -557,4 +557,5 @@ syn keyword juliaConstC StressbalanceRiftPenaltyThresholdEnum syn keyword juliaConstC StressbalanceShelfDampeningEnum +syn keyword juliaConstC ThermalForcingMonthlyEffectsEnum syn keyword juliaConstC ThermalIsdrainicecolumnEnum syn keyword juliaConstC ThermalIsdynamicbasalspcEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27284) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27285) @@ -578,4 +578,5 @@ else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum; else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum; + else if (strcmp(name,"ThermalForcingMonthlyEffects")==0) return ThermalForcingMonthlyEffectsEnum; else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum; else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"Xxe")==0) return XxeEnum; else if (strcmp(name,"Yye")==0) return YyeEnum; - else if (strcmp(name,"Zze")==0) return ZzeEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"Areae")==0) return AreaeEnum; + if (strcmp(name,"Zze")==0) return ZzeEnum; + else if (strcmp(name,"Areae")==0) return AreaeEnum; else if (strcmp(name,"WorldComm")==0) return WorldCommEnum; else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"DrivingStressX")==0) return DrivingStressXEnum; else if (strcmp(name,"DrivingStressY")==0) return DrivingStressYEnum; - else if (strcmp(name,"Dummy")==0) return DummyEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum; + if (strcmp(name,"Dummy")==0) return DummyEnum; + else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum; else if (strcmp(name,"EffectivePressureSubstep")==0) return EffectivePressureSubstepEnum; else if (strcmp(name,"EffectivePressureTransient")==0) return EffectivePressureTransientEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Misfit")==0) return MisfitEnum; else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum; - else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; + if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum; + else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum; else if (strcmp(name,"NewDamage")==0) return NewDamageEnum; else if (strcmp(name,"Node")==0) return NodeEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum; else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum; - else if (strcmp(name,"SmbD")==0) return SmbDEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum; + if (strcmp(name,"SmbD")==0) return SmbDEnum; + else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum; else if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum; else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum; @@ -1120,9 +1121,9 @@ 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=10; } if(stage==10){ - if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; + if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; + else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum; else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"Outputdefinition68")==0) return Outputdefinition68Enum; else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum; - else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum; else stage=11; } if(stage==11){ - if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; + if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum; + else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum; else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum; else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum; @@ -1366,9 +1367,9 @@ else if (strcmp(name,"Element")==0) return ElementEnum; else if (strcmp(name,"ElementHook")==0) return ElementHookEnum; - else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; + if (strcmp(name,"ElementSId")==0) return ElementSIdEnum; + else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum; else if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum; @@ -1489,9 +1490,9 @@ else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum; else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum; - else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; else stage=13; } if(stage==13){ - if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; + if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum; + else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum; @@ -1612,9 +1613,9 @@ else if (strcmp(name,"SmoothAnalysis")==0) return SmoothAnalysisEnum; else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum; - else if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum; else stage=14; } if(stage==14){ - if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum; + if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum; + else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum; else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum; else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum; Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m (revision 27284) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m (revision 27285) @@ -15,4 +15,5 @@ arlag_coefs = NaN; malag_coefs = NaN; + monthly_effects = NaN; basin_id = NaN; subglacial_discharge = NaN; @@ -64,4 +65,10 @@ md = checkfield(md,'fieldname','frontalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]); md = checkfield(md,'fieldname','frontalforcings.malag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ma_order]); + if(any(~isnan(md.frontalforcings.monthly_effects))) + md = checkfield(md,'fieldname','frontalforcings.monthly_effects','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,12]); + if(md.timestepping.time_step>=1) + error('md.frontalforcings.monthly_effects are provided but md.timestepping.time_step>=1'); + end + end end % }}} @@ -79,7 +86,14 @@ fielddisplay(self,'arlag_coefs','basin-specific vectors of AR lag coefficients [unitless]'); fielddisplay(self,'malag_coefs','basin-specific vectors of MA lag coefficients [unitless]'); + fielddisplay(self,'monthly_effects','basin-specific monthly values of TF added at corresponding month (default: all 0) [°C]'); end % }}} function marshall(self,prefix,md,fid) % {{{ yts=md.constants.yts; + if(any(isnan(md.frontalforcings.monthly_effects))) %monthly effects not provided, set to 0 + meffects = zeros(md.frontalforcings.num_basins,12); + else + meffects = md.frontalforcings.monthly_effects; + end + WriteData(fid,prefix,'name','md.frontalforcings.parameterization','data',3,'format','Integer'); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_basins','format','Integer'); @@ -94,4 +108,5 @@ WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','arlag_coefs','format','DoubleMat','name','md.frontalforcings.arlag_coefs','yts',yts); WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','malag_coefs','format','DoubleMat','name','md.frontalforcings.malag_coefs','yts',yts); + WriteData(fid,prefix,'data',meffects,'name','md.frontalforcings.monthly_effects','format','DoubleMat'); end % }}} end Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py (revision 27284) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py (revision 27285) @@ -22,4 +22,6 @@ self.arma_timestep = 0 self.arlag_coefs = np.nan + self.malag_coefs = np.nan + self.monthly_effects = np.nan self.basin_id = np.nan self.subglacial_discharge = np.nan @@ -43,4 +45,5 @@ s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of AR lag coefficients [unitless]')) s += '{}\n'.format(fielddisplay(self, 'malag_coefs', 'basin-specific vectors of MA lag coefficients [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'monthly_effects', 'basin-specific monthly values of TF added at corresponding month (default: all 0) [°C]')) return s #}}} @@ -74,4 +77,8 @@ md = checkfield(md, 'fieldname', 'frontalforcings.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ar_order]) md = checkfield(md, 'fieldname', 'frontalforcings.malag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ma_order]) + if(np.any(np.isnan(md.frontalforcings.monthly_effects)==False)): + md = checkfield(md, 'fieldname', 'frontalforcings.monthly_effects', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, 12]) + if(md.timestepping.time_step>=1): + raise RuntimeError('md.frontalforcings.monthly_effects are provided but md.timestepping.time_step>=1') return md # }}} @@ -84,4 +91,8 @@ def marshall(self, prefix, md, fid): # {{{ yts = md.constants.yts + if(np.any(np.isnan(md.frontalforcings.monthly_effects))): #monthly effects not provided, set to 0 + meffects = np.zeros((md.frontalforcings.num_basins,12)) + else: + meffects = 1*md.frontalforcings.monthly_effects WriteData(fid, prefix, 'name', 'md.frontalforcings.parameterization', 'data', 3, 'format', 'Integer') WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'num_basins', 'format', 'Integer') @@ -96,3 +107,4 @@ WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.arlag_coefs', 'yts', yts) WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'malag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.malag_coefs', 'yts', yts) + WriteData(fid, prefix, 'data', meffects, 'name', 'md.frontalforcings.monthly_effects', 'format', 'DoubleMat') # }}}