Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27333) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27334) @@ -248,4 +248,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_params",FrontalForcingsNumberofParamsEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_breaks",FrontalForcingsNumberofBreaksEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.monthlyvals_numbreaks",FrontalForcingsNumberofMonthBreaksEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_order",FrontalForcingsARMAarOrderEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ma_order",FrontalForcingsARMAmaOrderEnum)); @@ -263,6 +264,12 @@ 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)); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.monthlyvals_datebreaks"); + parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAmonthdatebreaksEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.monthlyvals_intercepts"); + parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAmonthinterceptsEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.monthlyvals_trends"); + parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAmonthtrendsEnum,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 27333) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27334) @@ -2714,13 +2714,16 @@ }/*}}}*/ -void Element::MonthlyEffectBasin(IssmDouble* monthlyeff, int enum_type){/*{{{*/ +void Element::MonthlyPiecewiseLinearEffectBasin(int nummonthbreaks,IssmDouble* monthlyintercepts,IssmDouble* monthlytrends,IssmDouble* monthlydatebreaks,int enum_type){/*{{{*/ /*Variable declaration*/ const int numvertices = this->GetNumberOfVertices(); - int basinid,mindex,basinenum_type,varenum_type; + int basinid,mindex,basinenum_type,varenum_type,indperiod; + int numperiods = nummonthbreaks+1; 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); + IssmDouble* datebreaks_b = xNew(nummonthbreaks); + IssmDouble* intercepts_b = xNew(nummonthbreaks*12); + IssmDouble* trends_b = xNew(nummonthbreaks*12); + IssmDouble* varlist = xNew(numvertices); /*Get field-specific enums*/ @@ -2743,7 +2746,26 @@ /*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]; - } + if(nummonthbreaks>0){ + for(int i=0;i0){ + /*Find index of time compared to the breakpoints*/ + indperiod = 0; + for(int i=0;idatebreaks_b[i]) indperiod = i+1; + } + /*Compute intercept+trend terms with parameters of indperiod*/ + if(indperiod==0) telapsed_break = time; + else telapsed_break = time-datebreaks_b[indperiod]; + piecewiselin = intercepts_b[indperiod]+trends_b[indperiod]*telapsed_break; + } + else piecewiselin = intercepts_b[indperiod]+trends_b[indperiod]*time; /*Retrieve values non-adjusted for monthly effects*/ @@ -2751,7 +2773,5 @@ /*Adjust values using monthly effects*/ - for(int v=0;v(monthlyeff_b); + xDelete(datebreaks_b); + xDelete(intercepts_b); + xDelete(trends_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 27333) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 27334) @@ -164,5 +164,5 @@ void MigrateGroundingLine(IssmDouble* sheet_ungrounding); void MismipFloatingiceMeltingRate(); - void MonthlyEffectBasin(IssmDouble* monthlyeff, int enum_type); + void MonthlyPiecewiseLinearEffectBasin(int nummonthbreaks,IssmDouble* monthlyintercepts,IssmDouble* monthlytrends,IssmDouble* monthlydatebreaks,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 27333) +++ /issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (revision 27334) @@ -50,21 +50,27 @@ bool isstochastic; bool istfstochastic = false; - int M,N,arorder,maorder,numbasins,numparams,numbreaks,my_rank; + int M,N,arorder,maorder,numbasins,numparams,numbreaks,nummonthbreaks,my_rank; femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); femmodel->parameters->FindParam(&numparams,FrontalForcingsNumberofParamsEnum); femmodel->parameters->FindParam(&numbreaks,FrontalForcingsNumberofBreaksEnum); + femmodel->parameters->FindParam(&nummonthbreaks,FrontalForcingsNumberofMonthBreaksEnum); femmodel->parameters->FindParam(&arorder,FrontalForcingsARMAarOrderEnum); femmodel->parameters->FindParam(&maorder,FrontalForcingsARMAmaOrderEnum); - IssmDouble* datebreaks = NULL; - IssmDouble* arlagcoefs = NULL; - IssmDouble* malagcoefs = NULL; - IssmDouble* monthlyeff = NULL; - IssmDouble* polyparams = NULL; + IssmDouble* datebreaks = NULL; + IssmDouble* arlagcoefs = NULL; + IssmDouble* malagcoefs = NULL; + IssmDouble* monthlyeff = NULL; + IssmDouble* polyparams = NULL; + IssmDouble* monthdatebreaks = NULL; + IssmDouble* monthintercepts = NULL; + IssmDouble* monthtrends = NULL; - femmodel->parameters->FindParam(&datebreaks,&M,&N,FrontalForcingsARMAdatebreaksEnum); _assert_(M==numbasins); _assert_(N==max(numbreaks,1)); - femmodel->parameters->FindParam(&polyparams,&M,&N,FrontalForcingsARMApolyparamsEnum); _assert_(M==numbasins); _assert_(N==(numbreaks+1)*numparams); - femmodel->parameters->FindParam(&arlagcoefs,&M,&N,FrontalForcingsARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(N==arorder); - femmodel->parameters->FindParam(&malagcoefs,&M,&N,FrontalForcingsARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(N==maorder); - femmodel->parameters->FindParam(&monthlyeff,&M,&N,ThermalForcingMonthlyEffectsEnum); _assert_(M==numbasins); _assert_(N==12); + femmodel->parameters->FindParam(&datebreaks,&M,&N,FrontalForcingsARMAdatebreaksEnum); _assert_(M==numbasins); _assert_(N==max(numbreaks,1)); + femmodel->parameters->FindParam(&polyparams,&M,&N,FrontalForcingsARMApolyparamsEnum); _assert_(M==numbasins); _assert_(N==(numbreaks+1)*numparams); + femmodel->parameters->FindParam(&arlagcoefs,&M,&N,FrontalForcingsARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(N==arorder); + femmodel->parameters->FindParam(&malagcoefs,&M,&N,FrontalForcingsARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(N==maorder); + femmodel->parameters->FindParam(&monthdatebreaks,&M,&N,FrontalForcingsARMAmonthdatebreaksEnum); _assert_(M==numbasins); _assert_(N==max(numbreaks,1)); + femmodel->parameters->FindParam(&monthintercepts,&M,&N,FrontalForcingsARMAmonthinterceptsEnum); _assert_(M==numbasins); _assert_(N==12*(nummonthbreaks+1)); + femmodel->parameters->FindParam(&monthtrends,&M,&N,FrontalForcingsARMAmonthtrendsEnum); _assert_(M==numbasins); _assert_(N==12*(nummonthbreaks+1)); femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); @@ -86,5 +92,5 @@ element->ArmaProcess(isstepforarma,arorder,maorder,numparams,numbreaks,tstep_arma,polyparams,arlagcoefs,malagcoefs,datebreaks,istfstochastic,FrontalForcingsRignotarmaEnum); /*Compute monthly effects*/ - element->MonthlyEffectBasin(monthlyeff,FrontalForcingsRignotarmaEnum); + element->MonthlyPiecewiseLinearEffectBasin(nummonthbreaks,monthintercepts,monthtrends,monthdatebreaks,FrontalForcingsRignotarmaEnum); } Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27333) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27334) @@ -174,5 +174,4 @@ syn keyword cConstant DebrisPackingFractionEnum syn keyword cConstant DebugProfilingEnum -syn keyword cConstant DebrisThicknessEnum syn keyword cConstant DomainDimensionEnum syn keyword cConstant DomainTypeEnum @@ -218,7 +217,11 @@ syn keyword cConstant FrontalForcingsARMAmaOrderEnum syn keyword cConstant FrontalForcingsARMAdatebreaksEnum +syn keyword cConstant FrontalForcingsARMAmonthdatebreaksEnum +syn keyword cConstant FrontalForcingsARMAmonthinterceptsEnum +syn keyword cConstant FrontalForcingsARMAmonthtrendsEnum syn keyword cConstant FrontalForcingsARMApolyparamsEnum syn keyword cConstant FrontalForcingsNumberofBasinsEnum syn keyword cConstant FrontalForcingsNumberofBreaksEnum +syn keyword cConstant FrontalForcingsNumberofMonthBreaksEnum syn keyword cConstant FrontalForcingsNumberofParamsEnum syn keyword cConstant FrontalForcingsParamEnum @@ -724,4 +727,5 @@ syn keyword cConstant DamageDbarOldEnum syn keyword cConstant DamageFEnum +syn keyword cConstant DebrisThicknessEnum syn keyword cConstant DegreeOfChannelizationEnum syn keyword cConstant DepthBelowSurfaceEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27333) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27334) @@ -211,7 +211,11 @@ FrontalForcingsARMAmaOrderEnum, FrontalForcingsARMAdatebreaksEnum, + FrontalForcingsARMAmonthdatebreaksEnum, + FrontalForcingsARMAmonthinterceptsEnum, + FrontalForcingsARMAmonthtrendsEnum, FrontalForcingsARMApolyparamsEnum, FrontalForcingsNumberofBasinsEnum, FrontalForcingsNumberofBreaksEnum, + FrontalForcingsNumberofMonthBreaksEnum, FrontalForcingsNumberofParamsEnum, FrontalForcingsParamEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27333) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27334) @@ -219,7 +219,11 @@ case FrontalForcingsARMAmaOrderEnum : return "FrontalForcingsARMAmaOrder"; case FrontalForcingsARMAdatebreaksEnum : return "FrontalForcingsARMAdatebreaks"; + case FrontalForcingsARMAmonthdatebreaksEnum : return "FrontalForcingsARMAmonthdatebreaks"; + case FrontalForcingsARMAmonthinterceptsEnum : return "FrontalForcingsARMAmonthintercepts"; + case FrontalForcingsARMAmonthtrendsEnum : return "FrontalForcingsARMAmonthtrends"; case FrontalForcingsARMApolyparamsEnum : return "FrontalForcingsARMApolyparams"; case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins"; case FrontalForcingsNumberofBreaksEnum : return "FrontalForcingsNumberofBreaks"; + case FrontalForcingsNumberofMonthBreaksEnum : return "FrontalForcingsNumberofMonthBreaks"; case FrontalForcingsNumberofParamsEnum : return "FrontalForcingsNumberofParams"; case FrontalForcingsParamEnum : return "FrontalForcingsParam"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27333) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27334) @@ -167,5 +167,4 @@ syn keyword juliaConstC DebrisPackingFractionEnum syn keyword juliaConstC DebugProfilingEnum -syn keyword juliaConstC DebrisThicknessEnum syn keyword juliaConstC DomainDimensionEnum syn keyword juliaConstC DomainTypeEnum @@ -211,7 +210,11 @@ syn keyword juliaConstC FrontalForcingsARMAmaOrderEnum syn keyword juliaConstC FrontalForcingsARMAdatebreaksEnum +syn keyword juliaConstC FrontalForcingsARMAmonthdatebreaksEnum +syn keyword juliaConstC FrontalForcingsARMAmonthinterceptsEnum +syn keyword juliaConstC FrontalForcingsARMAmonthtrendsEnum syn keyword juliaConstC FrontalForcingsARMApolyparamsEnum syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum syn keyword juliaConstC FrontalForcingsNumberofBreaksEnum +syn keyword juliaConstC FrontalForcingsNumberofMonthBreaksEnum syn keyword juliaConstC FrontalForcingsNumberofParamsEnum syn keyword juliaConstC FrontalForcingsParamEnum @@ -717,4 +720,5 @@ syn keyword juliaConstC DamageDbarOldEnum syn keyword juliaConstC DamageFEnum +syn keyword juliaConstC DebrisThicknessEnum syn keyword juliaConstC DegreeOfChannelizationEnum syn keyword juliaConstC DepthBelowSurfaceEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27333) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27334) @@ -222,7 +222,11 @@ else if (strcmp(name,"FrontalForcingsARMAmaOrder")==0) return FrontalForcingsARMAmaOrderEnum; else if (strcmp(name,"FrontalForcingsARMAdatebreaks")==0) return FrontalForcingsARMAdatebreaksEnum; + else if (strcmp(name,"FrontalForcingsARMAmonthdatebreaks")==0) return FrontalForcingsARMAmonthdatebreaksEnum; + else if (strcmp(name,"FrontalForcingsARMAmonthintercepts")==0) return FrontalForcingsARMAmonthinterceptsEnum; + else if (strcmp(name,"FrontalForcingsARMAmonthtrends")==0) return FrontalForcingsARMAmonthtrendsEnum; else if (strcmp(name,"FrontalForcingsARMApolyparams")==0) return FrontalForcingsARMApolyparamsEnum; else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum; else if (strcmp(name,"FrontalForcingsNumberofBreaks")==0) return FrontalForcingsNumberofBreaksEnum; + else if (strcmp(name,"FrontalForcingsNumberofMonthBreaks")==0) return FrontalForcingsNumberofMonthBreaksEnum; else if (strcmp(name,"FrontalForcingsNumberofParams")==0) return FrontalForcingsNumberofParamsEnum; else if (strcmp(name,"FrontalForcingsParam")==0) return FrontalForcingsParamEnum; @@ -256,12 +260,12 @@ else if (strcmp(name,"HydrologydcEplMaxThickness")==0) return HydrologydcEplMaxThicknessEnum; else if (strcmp(name,"HydrologydcEplPoreWaterMass")==0) return HydrologydcEplPoreWaterMassEnum; - else if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum; + else stage=3; + } + if(stage==3){ + if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum; else if (strcmp(name,"HydrologydcEplflipLock")==0) return HydrologydcEplflipLockEnum; else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum; else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum; - else stage=3; - } - if(stage==3){ - if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum; + else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum; else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum; else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum; @@ -379,12 +383,12 @@ else if (strcmp(name,"OceanGridNx")==0) return OceanGridNxEnum; else if (strcmp(name,"OceanGridNy")==0) return OceanGridNyEnum; - else if (strcmp(name,"OceanGridX")==0) return OceanGridXEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"OceanGridX")==0) return OceanGridXEnum; else if (strcmp(name,"OceanGridY")==0) return OceanGridYEnum; else if (strcmp(name,"OutputBufferPointer")==0) return OutputBufferPointerEnum; else if (strcmp(name,"OutputBufferSizePointer")==0) return OutputBufferSizePointerEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; + else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; @@ -502,12 +506,12 @@ else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; else if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum; - else if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum; else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum; else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum; else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum; + else if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum; else if (strcmp(name,"SmbARMAarOrder")==0) return SmbARMAarOrderEnum; else if (strcmp(name,"SmbARMAmaOrder")==0) return SmbARMAmaOrderEnum; @@ -625,12 +629,12 @@ else if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum; else if (strcmp(name,"TransientIsage")==0) return TransientIsageEnum; - else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum; else if (strcmp(name,"TransientIsdebris")==0) return TransientIsdebrisEnum; else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum; else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; + else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; @@ -748,12 +752,12 @@ else if (strcmp(name,"DeltaDsl")==0) return DeltaDslEnum; else if (strcmp(name,"DslOld")==0) return DslOldEnum; - else if (strcmp(name,"Dsl")==0) return DslEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"Dsl")==0) return DslEnum; else if (strcmp(name,"DeltaStr")==0) return DeltaStrEnum; else if (strcmp(name,"StrOld")==0) return StrOldEnum; else if (strcmp(name,"Str")==0) return StrEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; + else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; @@ -871,12 +875,12 @@ else if (strcmp(name,"LevelsetObservation")==0) return LevelsetObservationEnum; else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum; - else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum; else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum; else if (strcmp(name,"MaskOceanLevelset")==0) return MaskOceanLevelsetEnum; else if (strcmp(name,"MaskIceLevelset")==0) return MaskIceLevelsetEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"MaskIceRefLevelset")==0) return MaskIceRefLevelsetEnum; + else if (strcmp(name,"MaskIceRefLevelset")==0) return MaskIceRefLevelsetEnum; else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum; else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum; @@ -994,12 +998,12 @@ else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum; else if (strcmp(name,"SigmaVM")==0) return SigmaVMEnum; - else if (strcmp(name,"SmbAccumulatedEC")==0) return SmbAccumulatedECEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"SmbAccumulatedEC")==0) return SmbAccumulatedECEnum; else if (strcmp(name,"SmbAccumulatedMassBalance")==0) return SmbAccumulatedMassBalanceEnum; else if (strcmp(name,"SmbAccumulatedMelt")==0) return SmbAccumulatedMeltEnum; else if (strcmp(name,"SmbAccumulatedPrecipitation")==0) return SmbAccumulatedPrecipitationEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"SmbAccumulatedRain")==0) return SmbAccumulatedRainEnum; + else if (strcmp(name,"SmbAccumulatedRain")==0) return SmbAccumulatedRainEnum; else if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum; else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum; @@ -1117,12 +1121,12 @@ else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; - else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; + else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum; else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum; @@ -1240,12 +1244,12 @@ else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum; else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; - else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; + else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum; @@ -1363,12 +1367,12 @@ else if (strcmp(name,"ControlInputValues")==0) return ControlInputValuesEnum; else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum; - else if (strcmp(name,"Cuffey")==0) return CuffeyEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"Cuffey")==0) return CuffeyEnum; else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum; else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"DataSet")==0) return DataSetEnum; + else if (strcmp(name,"DataSet")==0) return DataSetEnum; else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; @@ -1486,12 +1490,12 @@ else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; - else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; + else stage=13; + } + if(stage==13){ + if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; else if (strcmp(name,"LinearFloatingMeltRatearma")==0) return LinearFloatingMeltRatearmaEnum; else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; else if (strcmp(name,"Loads")==0) return LoadsEnum; - else stage=13; - } - if(stage==13){ - if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; + else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; else if (strcmp(name,"LoveHf")==0) return LoveHfEnum; else if (strcmp(name,"LoveHt")==0) return LoveHtEnum; @@ -1609,12 +1613,12 @@ else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; - else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; + else stage=14; + } + if(stage==14){ + if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; - else stage=14; - } - if(stage==14){ - if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; + else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m (revision 27333) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m (revision 27334) @@ -6,17 +6,20 @@ classdef frontalforcingsrignotarma properties (SetAccess=public) - num_basins = 0; - num_params = 0; - num_breaks = 0; - polynomialparams = NaN; - datebreaks = NaN; - ar_order = 0; - ma_order = 0; - arma_timestep = 0; - arlag_coefs = NaN; - malag_coefs = NaN; - monthly_effects = NaN; - basin_id = NaN; - subglacial_discharge = NaN; + num_basins = 0; + num_params = 0; + num_breaks = 0; + polynomialparams = NaN; + datebreaks = NaN; + ar_order = 0; + ma_order = 0; + arma_timestep = 0; + arlag_coefs = NaN; + malag_coefs = NaN; + monthlyvals_intercepts = NaN; + monthlyvals_trends = NaN; + monthlyvals_numbreaks = 0; + monthlyvals_datebreaks = NaN; + basin_id = NaN; + subglacial_discharge = NaN; end methods @@ -52,44 +55,77 @@ end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ - %Early return - if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end - - nbas = md.frontalforcings.num_basins; - nprm = md.frontalforcings.num_params; - nbrk = md.frontalforcings.num_breaks; - md = checkfield(md,'fieldname','frontalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); - md = checkfield(md,'fieldname','frontalforcings.num_params','numel',1,'NaN',1,'Inf',1,'>',0); - md = checkfield(md,'fieldname','frontalforcings.num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0); + %Early return + if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end + + nbas = md.frontalforcings.num_basins; + nprm = md.frontalforcings.num_params; + nbrk = md.frontalforcings.num_breaks; + nMbrk = md.frontalforcings.monthlyvals_numbreaks; + md = checkfield(md,'fieldname','frontalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','frontalforcings.num_params','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','frontalforcings.num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','frontalforcings.basin_id','Inf',1,'>=',0,'<=',md.frontalforcings.num_basins,'size',[md.mesh.numberofelements 1]); - md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'NaN',1,'Inf',1,'timeseries',1); - if(nbas>1 && nbrk>=1 && nprm>1) - md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm); - elseif(nbas==1) - md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm); - elseif(nbrk==0) - md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm); - elseif(nprm==1) - md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1],'numel',nbas*(nbrk+1)*nprm); - end - md = checkfield(md,'fieldname','frontalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); + md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'NaN',1,'Inf',1,'timeseries',1); + if(nbas>1 && nbrk>=1 && nprm>1) + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm); + elseif(nbas==1) + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm); + elseif(nbrk==0) + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm); + elseif(nprm==1) + md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1],'numel',nbas*(nbrk+1)*nprm); + end + md = checkfield(md,'fieldname','frontalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','frontalforcings.ma_order','numel',1,'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','frontalforcings.arma_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %ARMA time step cannot be finer than ISSM timestep - 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(nbrk>0) - md = checkfield(md,'fieldname','frontalforcings.datebreaks','NaN',1,'Inf',1,'size',[nbas,nbrk]); - elseif(numel(md.frontalforcings.datebreaks)==0 || all(isnan(md.frontalforcings.datebreaks))) - ; - else - error('md.frontalforcings.num_breaks is 0 but md.frontalforcings.datebreaks is not empty'); - end - 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 % }}} + 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(nbrk>0) + md = checkfield(md,'fieldname','frontalforcings.datebreaks','NaN',1,'Inf',1,'size',[nbas,nbrk]); + elseif(numel(md.frontalforcings.datebreaks)==0 || all(isnan(md.frontalforcings.datebreaks))) + ; + else + error('md.frontalforcings.num_breaks is 0 but md.frontalforcings.datebreaks is not empty'); + end + + %%% Check if some monthly forcings are provided %%% + if((numel(md.frontalforcings.monthlyvals_intercepts)>1 || ~isnan(md.frontalforcings.monthlyvals_intercepts)) || (numel(md.frontalforcings.monthlyvals_trends)>1 || ~isnan(md.frontalforcings.monthlyvals_trends)) || (numel(md.frontalforcings.monthlyvals_datebreaks)>1 || ~isnan(md.frontalforcings.monthlyvals_datebreaks))) + isMonthly = true; + else + isMonthly = false; + end + if(numel(md.frontalforcings.monthlyvals_trends)>1 || ~isnan(md.frontalforcings.monthlyvals_trends)) + isMonthlyTrend = true; + else + isMonthlyTrend = false; + end + if(isMonthly) + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_numbreaks','numel',1,'NaN',1,'Inf',1,'>=',0); + if(nbas>1 && nMbrk>=1) + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nbas,12,nMbrk+1],'numel',nbas*(nMbrk+1)*12); + if(isMonthlyTrend) + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nbas,12,nMbrk+1],'numel',nbas*(nMbrk+1)*12); + end + elseif(nbas==1) + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nMbrk+1,12],'numel',nbas*(nMbrk+1)*12); + if(isMonthlyTrend) + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nMbrk+1,12],'numel',nbas*(nMbrk+1)*12); + end + elseif(nMbrk==0) + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nbas,12],'numel',nbas*(nMbrk+1)*12); + if(isMonthlyTrend) + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nbas,12],'numel',nbas*(nMbrk+1)*12); + end + end + if(nMbrk>0) + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_datebreaks','NaN',1,'Inf',1,'size',[nbas,nMbrk]); + elseif(numel(md.frontalforcings.monthlyvals_datebreaks)==0 || all(isnan(md.frontalforcings.monthlyvals_datebreaks))) + ; + else + error('md.frontalforcings.monthlyvals_numbreaks is 0 but md.frontalforcings.monthlyvals_datebreaks is not empty'); + end + end + + end % }}} function disp(self) % {{{ disp(sprintf(' Frontalforcings parameters:')); @@ -107,11 +143,15 @@ 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]'); + fielddisplay(self,'monthlyvals_intercepts','intercept of basin-specific piecewise-linear monthly values of TF added at corresponding month (default: all 0) [°C]'); + fielddisplay(self,'monthlyvals_trends','trends in basin-specific piecewise-linear monthly values of TF added at corresponding month (default: all 0) [°C/yr]'); + fielddisplay(self,'monthlyvals_numbreaks','number of breakpoints in the piecewise-linear functions of monthly values'); + fielddisplay(self,'monthlyvals_datebreaks','dates at which the breakpoints in the piecewise-linear monthly values occur (1 row per basin)'); end % }}} function marshall(self,prefix,md,fid) % {{{ yts=md.constants.yts; - nbas = md.frontalforcings.num_basins; - nprm = md.frontalforcings.num_params; - nper = md.frontalforcings.num_breaks+1; + %%% Deal with polynomial %%% + nbas = md.frontalforcings.num_basins; + nprm = md.frontalforcings.num_params; + nper = md.frontalforcings.num_breaks+1; % Scale the parameters % polyparamsScaled = md.frontalforcings.polynomialparams; @@ -150,13 +190,40 @@ polyparams2dScaled = polyparamsScaled; end - 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 if(nper==1) %a single period (no break date) dbreaks = zeros(nbas,1); %dummy else dbreaks = md.frontalforcings.datebreaks; + end + + %%% Deal with monthly effects %%% + nMper = md.frontalforcings.monthlyvals_numbreaks+1; + if((numel(md.frontalforcings.monthlyvals_intercepts)<=1 && isnan(md.frontalforcings.monthlyvals_intercepts))) + interceptsM = zeros(nbas,12); %monthly intercepts not provided, set to 0 + trendsM = zeros(nbas,12); %set monthly trends also to 0 + else + interceptsM3d = md.frontalforcings.monthlyvals_intercepts; + if((numel(md.frontalforcings.monthlyvals_trends)<=1 && isnan(md.frontalforcings.monthlyvals_trends))) + trendsM3d = 0*interceptsM3d; %monthly trends not provided, set to 0 + else + trendsM3d = md.frontalforcings.monthlyvals_trends; + end + end + % Create 2D arrays from 3D arrays if needed % + if(nMper>1 && (numel(md.frontalforcings.monthlyvals_intercepts)>1 || ~isnan(md.frontalforcings.monthlyvals_intercepts))) + interceptsM = zeros(nbas,12*nMper); + trendsM = zeros(nbas,12*nMper); + for(ii=[1:nMper]) + jj = 1+(ii-1)*12; + interceptsM(:,jj:jj+12-1) = interceptsM3d(:,:,ii); + trendsM(:,jj:jj+12-1) = trendsM3d(:,:,ii); + end + elseif(nMper==1 && (numel(md.frontalforcings.monthlyvals_intercepts)>1 || ~isnan(md.frontalforcings.monthlyvals_intercepts))) + interceptsM = interceptsM3d; + trendsM = trendsM3d; + end + if(nMper==1) %a single period (no break date) + dMbreaks = zeros(nbas,1); %dummy + else + dMbreaks = md.frontalforcings.monthlyvals_datebreaks; end @@ -174,5 +241,8 @@ WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','malag_coefs','format','DoubleMat','name','md.frontalforcings.malag_coefs','yts',yts); WriteData(fid,prefix,'data',dbreaks,'name','md.frontalforcings.datebreaks','format','DoubleMat','scale',yts); - WriteData(fid,prefix,'data',meffects,'name','md.frontalforcings.monthly_effects','format','DoubleMat'); + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','monthlyvals_numbreaks','format','Integer'); + WriteData(fid,prefix,'data',dMbreaks,'name','md.frontalforcings.monthlyvals_datebreaks','format','DoubleMat','scale',yts); + WriteData(fid,prefix,'data',interceptsM,'name','md.frontalforcings.monthlyvals_intercepts','format','DoubleMat'); + WriteData(fid,prefix,'data',trendsM,'name','md.frontalforcings.monthlyvals_trends','format','DoubleMat','scale',1/yts); end % }}} end Index: /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py =================================================================== --- /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py (revision 27333) +++ /issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py (revision 27334) @@ -25,5 +25,8 @@ self.arlag_coefs = np.nan self.malag_coefs = np.nan - self.monthly_effects = np.nan + self.monthlyvals_intercepts = np.nan + self.monthlyvals_trends = np.nan + self.monthlyvals_numbreaks = 0 + self.monthlyvals_datebreaks = np.nan self.basin_id = np.nan self.subglacial_discharge = np.nan @@ -48,5 +51,4 @@ 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 #}}} @@ -66,7 +68,8 @@ return md - nbas = md.frontalforcings.num_basins; - nprm = md.frontalforcings.num_params; - nbrk = md.frontalforcings.num_breaks; + nbas = md.frontalforcings.num_basins; + nprm = md.frontalforcings.num_params; + nbrk = md.frontalforcings.num_breaks; + nMbrk = md.frontalforcings.monthlyvals_numbreaks; md = checkfield(md, 'fieldname', 'frontalforcings.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) md = checkfield(md, 'fieldname', 'frontalforcings.num_params', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) @@ -95,8 +98,36 @@ else: raise RuntimeError('md.frontalforcings.num_breaks is 0 but md.frontalforcings.datebreaks is not empty') - 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') + + + ### Check if some monthly forcings are provided ### + if(np.all(np.isnan(md.frontalforcings.monthlyvals_intercepts))==False or np.all(np.isnan(md.frontalforcings.monthlyvals_trends))==False or np.all(np.isnan(md.frontalforcings.monthlyvals_datebreaks))==False): + isMonthly = True + else: + isMonthly = False + if(np.all(np.isnan(md.frontalforcings.monthlyvals_datebreaks))): + isMonthlyTrend = True + else: + isMonthlyTrend = False + if(isMonthly): + md = checkfield(md, 'fieldname', 'frontalforcings.monthlyvals_numbreaks', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) + if(nbas>1 and nMbrk>=1): + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nbas,12,nMbrk+1],'numel',nbas*(nMbrk+1)*12) + if(isMonthlyTrend): + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nbas,12,nMbrk+1],'numel',nbas*(nMbrk+1)*12) + elif(nbas==1): + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nMbrk+1,12],'numel',nbas*(nMbrk+1)*12) + if(isMonthlyTrend): + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nMbrk+1,12],'numel',nbas*(nMbrk+1)*12) + elif(nMbrk==0): + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nbas,12],'numel',nbas*(nMbrk+1)*12) + if(isMonthlyTrend): + md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nbas,12],'numel',nbas*(nMbrk+1)*12) + if(nMbrk>0): + md = checkfield(md, 'fieldname', 'frontalforcings.monthlyvals_datebreaks', 'NaN', 1, 'Inf', 1, 'size', [nbas,nMbrk]) + elif(np.size(md.frontalforcings.monthlyvals_datebreaks)==0 or np.all(np.isnan(md.frontalforcings.monthlyvals_datebreaks))): + pass + else: + raise RuntimeError('md.frontalforcings.monthlyvals_numbreaks is 0 but md.frontalforcings.monthlyvals_datebreaks is not empty') + return md # }}} @@ -139,12 +170,34 @@ # 2D array is already in correct format and no need for scaling # polyparams2dScaled = np.copy(polyparamsScaled) - 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 if(nper==1): dbreaks = np.zeros((nbas,1)) else: dbreaks = np.copy(md.frontalforcings.datebreaks) + + ### Deal with montly effects ### + nMper = md.frontalforcings.monthlyvals_numbreaks+1 + if(np.any(np.isnan(md.frontalforcings.monthlyvals_intercepts))): + interceptsM = np.zeros((nbas,12)) #monthly intercepts not provided, set to 0 + trendsM = np.zeros((nbas,12)) #set monthly trends also to 0 + else: + interceptsM3d = md.frontalforcings.monthlyvals_intercepts + if(np.any(np.isnan(md.frontalforcings.monthlyvals_trends))): + trendsM3d = 0*interceptsM3d #monthly trends not provided, set to 0 + else: + trendsM3d = md.frontalforcings.monthlyvals_trends + # Create 2D arrays from 3D arrays if needed # + if(nMper>1 and np.all(np.isnan(md.frontalforcings.monthlyvals_intercepts))==False): + interceptsM = np.zeros((nbas,12*nMper)) + trendsM = np.zeros((nbas,12*nMper)) + for ii in range(nMper): + interceptsM[:,ii*12:(ii+1)*12] = 1*interceptsM3d[:,:,ii] + trendsM[:,ii*12:(ii+1)*12] = 1*trendsM3d[:,:,ii] + elif(nMper==1 and np.all(np.isnan(md.frontalforcings.monthlyvals_intercepts))==False): + interceptsM = 1*interceptsM3d + trendsM = 1*trendsM3d + if(nMper==1): + dMbreaks = np.zeros((nbas,1)) + else: + dMbreaks = np.copy(md.frontalforcings.monthlyvals_datebreaks) WriteData(fid, prefix, 'name', 'md.frontalforcings.parameterization', 'data', 3, 'format', 'Integer') @@ -161,4 +214,9 @@ WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'malag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.malag_coefs', 'yts', yts) WriteData(fid, prefix, 'data', dbreaks, 'name', 'md.frontalforcings.datebreaks', 'format', 'DoubleMat','scale',yts) - WriteData(fid, prefix, 'data', meffects, 'name', 'md.frontalforcings.monthly_effects', 'format', 'DoubleMat') + WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','monthlyvals_numbreaks','format','Integer') + WriteData(fid,prefix,'data',dMbreaks,'name','md.frontalforcings.monthlyvals_datebreaks','format','DoubleMat','scale',yts) + WriteData(fid,prefix,'data',interceptsM,'name','md.frontalforcings.monthlyvals_intercepts','format','DoubleMat') + WriteData(fid,prefix,'data',trendsM,'name','md.frontalforcings.monthlyvals_trends','format','DoubleMat','scale',1/yts) # }}} + + Index: /issm/trunk-jpl/test/NightlyRun/test543.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.m (revision 27333) +++ /issm/trunk-jpl/test/NightlyRun/test543.m (revision 27334) @@ -5,6 +5,6 @@ md=setflowequation(md,'SSA','all'); md.timestepping.start_time = 0; -md.timestepping.time_step = 1; -md.timestepping.final_time = 10; +md.timestepping.time_step = 0.05; +md.timestepping.final_time = 2; %Basin separation TF @@ -52,20 +52,36 @@ trendlin = [-0.5,-0.2,0.1;0,0,0]; trendquad = [0,0.0,0;0.1,0.1,0.1]; -datebreaks = [4,7;4,7]; +datebreaks = [0.5,1.5;0.5,1.5]; polynomialparams = cat(numparams,intercept,trendlin,trendquad); +% Monthly effects params % +numbreaksM = 1; +intcpsMp0 = [-0.5,-0.5,0,0,0,0,0.5,0.5,0,0,0,0; + -1.0,-1.0,0,0,0,0,1.0,1.0,0,0,0,0]; +intcpsMp1 = [-0.25,-0.25,0,0,0,0,0.,0.,0,0,0,0; + -0.1,-0.1,0,0,0,0,0.1,0.1,0,0,0,0]; +intcpsM = cat(3,intcpsMp0,intcpsMp1); +trendsMp0 = [0,0,0,0,0,0,0.,0.0,0,0,0,0; + 0.0,0.0,0,-0.0,0,0,0.0,0.0,0,0,0,0]; +trendsMp1 = [0,-0.12,0,0,0,0,0.,0.0,0,0.0,0,0; + 0.0,-0.1,0,-0.0,0,0,0.0,0.0,0,0,0,0]; +trendsM = cat(3,trendsMp0,trendsMp1); +datebreaksM = [1;1]; -md.frontalforcings.num_basins = nb_tf; -md.frontalforcings.basin_id = idb_tf; -md.frontalforcings.num_params = numparams; %number of parameters in the polynomial -md.frontalforcings.num_breaks = numbreaks; %number of breakpoints -md.frontalforcings.subglacial_discharge = 0.01*ones(md.mesh.numberofvertices,1); -md.frontalforcings.polynomialparams = polynomialparams; -md.frontalforcings.datebreaks = datebreaks; -md.frontalforcings.ar_order = 4; -md.frontalforcings.ma_order = 2; -md.frontalforcings.arma_timestep = 2; %timestep of the ARMA model [yr] -md.frontalforcings.arlag_coefs = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters -md.frontalforcings.malag_coefs = [[0.1,0.0];[0.0,0.1]]; %moving-average parameters - +md.frontalforcings.num_basins = nb_tf; +md.frontalforcings.basin_id = idb_tf; +md.frontalforcings.num_params = numparams; %number of parameters in the polynomial +md.frontalforcings.num_breaks = numbreaks; %number of breakpoints +md.frontalforcings.subglacial_discharge = 0.01*ones(md.mesh.numberofvertices,1); +md.frontalforcings.polynomialparams = polynomialparams; +md.frontalforcings.datebreaks = datebreaks; +md.frontalforcings.ar_order = 4; +md.frontalforcings.ma_order = 2; +md.frontalforcings.arma_timestep = 2; %timestep of the ARMA model [yr] +md.frontalforcings.arlag_coefs = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters +md.frontalforcings.malag_coefs = [[0.1,0.0];[0.0,0.1]]; %moving-average parameters +md.frontalforcings.monthlyvals_numbreaks = numbreaksM; +md.frontalforcings.monthlyvals_intercepts = intcpsM; +md.frontalforcings.monthlyvals_trends = trendsM; +md.frontalforcings.monthlyvals_datebreaks = datebreaksM; % Floating Ice Melt parameters md.basalforcings.floatingice_melting_rate = 0.1*ones(md.mesh.numberofvertices,1); @@ -102,5 +118,5 @@ 1e-11,2e-11,2e-11,1e-11,1e-9,1e-10,1e-10,1e-10,... 2e-11,1e-11,1e-11,9e-11,2e-9,1e-10,1e-10,1e-10,... - 2e-10,1e-10,1e-10,1e-10,5e-9,1e-10,1e-10,1e-10,... + 2e-6,1e-6,1e-6,1e-6,5e-6,1e-6,1e-6,1e-6,... }; field_values={... @@ -113,19 +129,19 @@ (md.results.TransientSolution(1).CalvingMeltingrate),... (md.results.TransientSolution(1).BasalforcingsFloatingiceMeltingRate),... - (md.results.TransientSolution(5).Vx),... - (md.results.TransientSolution(5).Vy),... - (md.results.TransientSolution(5).Vel),... - (md.results.TransientSolution(5).Thickness),... - (md.results.TransientSolution(5).MaskIceLevelset),... - (md.results.TransientSolution(5).CalvingCalvingrate),... - (md.results.TransientSolution(5).CalvingMeltingrate),... - (md.results.TransientSolution(5).BasalforcingsFloatingiceMeltingRate),... - (md.results.TransientSolution(10).Vx),... - (md.results.TransientSolution(10).Vy),... - (md.results.TransientSolution(10).Vel),... - (md.results.TransientSolution(10).Thickness),... - (md.results.TransientSolution(10).MaskIceLevelset),... - (md.results.TransientSolution(10).CalvingCalvingrate),... - (md.results.TransientSolution(10).CalvingMeltingrate),... - (md.results.TransientSolution(10).BasalforcingsFloatingiceMeltingRate),... + (md.results.TransientSolution(20).Vx),... + (md.results.TransientSolution(20).Vy),... + (md.results.TransientSolution(20).Vel),... + (md.results.TransientSolution(20).Thickness),... + (md.results.TransientSolution(20).MaskIceLevelset),... + (md.results.TransientSolution(20).CalvingCalvingrate),... + (md.results.TransientSolution(20).CalvingMeltingrate),... + (md.results.TransientSolution(20).BasalforcingsFloatingiceMeltingRate),... + (md.results.TransientSolution(40).Vx),... + (md.results.TransientSolution(40).Vy),... + (md.results.TransientSolution(40).Vel),... + (md.results.TransientSolution(40).Thickness),... + (md.results.TransientSolution(40).MaskIceLevelset),... + (md.results.TransientSolution(40).CalvingCalvingrate),... + (md.results.TransientSolution(40).CalvingMeltingrate),... + (md.results.TransientSolution(40).BasalforcingsFloatingiceMeltingRate),... }; Index: /issm/trunk-jpl/test/NightlyRun/test543.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test543.py (revision 27333) +++ /issm/trunk-jpl/test/NightlyRun/test543.py (revision 27334) @@ -16,6 +16,6 @@ md = setflowequation(md, 'SSA', 'all') md.timestepping.start_time = 0 -md.timestepping.time_step = 1 -md.timestepping.final_time = 10 +md.timestepping.time_step = 0.05 +md.timestepping.final_time = 2 # Basin separation TF @@ -57,6 +57,19 @@ trendlin = np.array([[-0.5,-0.2,0.1],[0,0,0]]) trendquad = np.array([[0,0.0,0],[0.1,0.1,0.1]]) -datebreaks = np.array([[4,7],[4,7]]) +datebreaks = np.array([[0.5,1.5],[0.5,1.5]]) polynomialparams = np.transpose(np.stack((intercept,trendlin,trendquad)),(1,2,0)) +# Monthly effects params # +numbreaksM = 1; +intcpsMp0 = np.array([[-0.5,-0.5,0,0,0,0,0.5,0.5,0,0,0,0], + [-1.0,-1.0,0,0,0,0,1.0,1.0,0,0,0,0]]) +intcpsMp1 = np.array([[-0.25,-0.25,0,0,0,0,0.,0.,0,0,0,0], + [-0.1,-0.1,0,0,0,0,0.1,0.1,0,0,0,0]]) +intcpsM = np.transpose(np.stack((intcpsMp0,intcpsMp1)),(1,2,0)) +trendsMp0 = np.array([[0,0,0,0,0,0,0.,0.0,0,0,0,0], + [0.0,0.0,0,-0.0,0,0,0.0,0.0,0,0,0,0]]) +trendsMp1 = np.array([[0,-0.12,0,0,0,0,0.,0.0,0,0.0,0,0], + [0.0,-0.1,0,-0.0,0,0,0.0,0.0,0,0,0,0]]) +trendsM = np.transpose(np.stack((trendsMp0,trendsMp1)),(1,2,0)) +datebreaksM = np.array([[1],[1]]) md.frontalforcings.num_basins = nb_tf @@ -72,4 +85,8 @@ md.frontalforcings.arlag_coefs = np.array([[0.1, -0.1, 0.01, -0.01], [0.2, -0.2, 0.1, 0.0]]) # autoregressive parameters md.frontalforcings.malag_coefs = np.array([[0.1, 0.0], [0.0, 0.1]]) # moving-average parameters +md.frontalforcings.monthlyvals_numbreaks = numbreaksM +md.frontalforcings.monthlyvals_intercepts = intcpsM +md.frontalforcings.monthlyvals_trends = trendsM +md.frontalforcings.monthlyvals_datebreaks = datebreaksM #Floating Ice Melt parameters md.basalforcings.floatingice_melting_rate = 0.1 * np.ones((md.mesh.numberofvertices,)) @@ -119,5 +136,5 @@ 1e-11, 2e-11, 2e-11, 1e-11, 1e-9, 1e-10, 1e-10, 1e-10, 2e-11, 1e-11, 1e-11, 9e-11, 2e-9, 1e-10, 1e-10, 1e-10, - 2e-10, 1e-10, 1e-10, 1e-10, 5e-9, 1e-10, 1e-10, 1e-10] + 2e-6, 1e-6, 1e-6, 1e-6, 5e-6, 1e-6, 1e-6, 1e-6] field_values = [ md.results.TransientSolution[0].Vx, @@ -129,18 +146,18 @@ md.results.TransientSolution[0].CalvingMeltingrate, md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate, - md.results.TransientSolution[4].Vx, - md.results.TransientSolution[4].Vy, - md.results.TransientSolution[4].Vel, - md.results.TransientSolution[4].Thickness, - md.results.TransientSolution[4].MaskIceLevelset, - md.results.TransientSolution[4].CalvingCalvingrate, - md.results.TransientSolution[4].CalvingMeltingrate, - md.results.TransientSolution[4].BasalforcingsFloatingiceMeltingRate, - md.results.TransientSolution[9].Vx, - md.results.TransientSolution[9].Vy, - md.results.TransientSolution[9].Vel, - md.results.TransientSolution[9].Thickness, - md.results.TransientSolution[9].MaskIceLevelset, - md.results.TransientSolution[9].CalvingCalvingrate, - md.results.TransientSolution[9].CalvingMeltingrate, - md.results.TransientSolution[9].BasalforcingsFloatingiceMeltingRate] + md.results.TransientSolution[19].Vx, + md.results.TransientSolution[19].Vy, + md.results.TransientSolution[19].Vel, + md.results.TransientSolution[19].Thickness, + md.results.TransientSolution[19].MaskIceLevelset, + md.results.TransientSolution[19].CalvingCalvingrate, + md.results.TransientSolution[19].CalvingMeltingrate, + md.results.TransientSolution[19].BasalforcingsFloatingiceMeltingRate, + md.results.TransientSolution[39].Vx, + md.results.TransientSolution[39].Vy, + md.results.TransientSolution[39].Vel, + md.results.TransientSolution[39].Thickness, + md.results.TransientSolution[39].MaskIceLevelset, + md.results.TransientSolution[39].CalvingCalvingrate, + md.results.TransientSolution[39].CalvingMeltingrate, + md.results.TransientSolution[39].BasalforcingsFloatingiceMeltingRate]