![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 27334
- Timestamp:
- 10/26/22 06:26:59 (2 years ago)
- Location:
- issm/trunk-jpl
- Files:
-
- 14 edited
-
src/c/analyses/LevelsetAnalysis.cpp (modified) (2 diffs)
-
src/c/classes/Elements/Element.cpp (modified) (4 diffs)
-
src/c/classes/Elements/Element.h (modified) (1 diff)
-
src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (modified) (2 diffs)
-
src/c/shared/Enum/Enum.vim (modified) (3 diffs)
-
src/c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
-
src/c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
-
src/c/shared/Enum/Enumjl.vim (modified) (3 diffs)
-
src/c/shared/Enum/StringToEnumx.cpp (modified) (13 diffs)
-
src/m/classes/frontalforcingsrignotarma.m (modified) (5 diffs)
-
src/m/classes/frontalforcingsrignotarma.py (modified) (6 diffs)
-
test/Archives/Archive543.arch (modified) ( previous)
-
test/NightlyRun/test543.m (modified) (4 diffs)
-
test/NightlyRun/test543.py (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp
r27318 r27334 248 248 parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_params",FrontalForcingsNumberofParamsEnum)); 249 249 parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_breaks",FrontalForcingsNumberofBreaksEnum)); 250 parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.monthlyvals_numbreaks",FrontalForcingsNumberofMonthBreaksEnum)); 250 251 parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ar_order",FrontalForcingsARMAarOrderEnum)); 251 252 parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.ma_order",FrontalForcingsARMAmaOrderEnum)); … … 263 264 parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAmalagcoefsEnum,transparam,M,N)); 264 265 xDelete<IssmDouble>(transparam); 265 iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.monthly_effects"); 266 parameters->AddObject(new DoubleMatParam(ThermalForcingMonthlyEffectsEnum,transparam,M,N)); 266 iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.monthlyvals_datebreaks"); 267 parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAmonthdatebreaksEnum,transparam,M,N)); 268 xDelete<IssmDouble>(transparam); 269 iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.monthlyvals_intercepts"); 270 parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAmonthinterceptsEnum,transparam,M,N)); 271 xDelete<IssmDouble>(transparam); 272 iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.monthlyvals_trends"); 273 parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAmonthtrendsEnum,transparam,M,N)); 267 274 xDelete<IssmDouble>(transparam); 268 275 /*Do not break here, generic FrontalForcingsRignot parameters still to be retrieved*/ -
issm/trunk-jpl/src/c/classes/Elements/Element.cpp
r27318 r27334 2714 2714 2715 2715 }/*}}}*/ 2716 void Element::Monthly EffectBasin(IssmDouble* monthlyeff,int enum_type){/*{{{*/2716 void Element::MonthlyPiecewiseLinearEffectBasin(int nummonthbreaks,IssmDouble* monthlyintercepts,IssmDouble* monthlytrends,IssmDouble* monthlydatebreaks,int enum_type){/*{{{*/ 2717 2717 2718 2718 /*Variable declaration*/ 2719 2719 const int numvertices = this->GetNumberOfVertices(); 2720 int basinid,mindex,basinenum_type,varenum_type; 2720 int basinid,mindex,basinenum_type,varenum_type,indperiod; 2721 int numperiods = nummonthbreaks+1; 2721 2722 IssmDouble time,fracyear,yts; 2722 2723 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}; 2723 IssmDouble* monthlyeff_b = xNew<IssmDouble>(12); 2724 IssmDouble* varlist = xNew<IssmDouble>(numvertices); 2724 IssmDouble* datebreaks_b = xNew<IssmDouble>(nummonthbreaks); 2725 IssmDouble* intercepts_b = xNew<IssmDouble>(nummonthbreaks*12); 2726 IssmDouble* trends_b = xNew<IssmDouble>(nummonthbreaks*12); 2727 IssmDouble* varlist = xNew<IssmDouble>(numvertices); 2725 2728 2726 2729 /*Get field-specific enums*/ … … 2743 2746 /*Get basin-specific parameters of the element*/ 2744 2747 this->GetInputValue(&basinid,basinenum_type); 2745 for(int ii=0;ii<12;ii++){ 2746 monthlyeff_b[ii] = monthlyeff[basinid*12+ii]; 2747 } 2748 if(nummonthbreaks>0){ 2749 for(int i=0;i<nummonthbreaks;i++) datebreaks_b[i] = monthlydatebreaks[basinid*nummonthbreaks+i]; 2750 } 2751 for(int i=0;i<numperiods;i++){ 2752 intercepts_b[i] = monthlyintercepts[basinid*12*numperiods+mindex+12*i]; 2753 trends_b[i] = monthlytrends[basinid*12*numperiods+mindex+12*i]; 2754 } 2755 2756 /*Compute piecewise-linear function*/ 2757 IssmDouble telapsed_break,piecewiselin; 2758 if(nummonthbreaks>0){ 2759 /*Find index of time compared to the breakpoints*/ 2760 indperiod = 0; 2761 for(int i=0;i<nummonthbreaks;i++){ 2762 if(time>datebreaks_b[i]) indperiod = i+1; 2763 } 2764 /*Compute intercept+trend terms with parameters of indperiod*/ 2765 if(indperiod==0) telapsed_break = time; 2766 else telapsed_break = time-datebreaks_b[indperiod]; 2767 piecewiselin = intercepts_b[indperiod]+trends_b[indperiod]*telapsed_break; 2768 } 2769 else piecewiselin = intercepts_b[indperiod]+trends_b[indperiod]*time; 2748 2770 2749 2771 /*Retrieve values non-adjusted for monthly effects*/ … … 2751 2773 2752 2774 /*Adjust values using monthly effects*/ 2753 for(int v=0;v<numvertices;v++){ 2754 varlist[v] = varlist[v]+monthlyeff_b[mindex]; 2755 } 2775 for(int v=0;v<numvertices;v++) varlist[v] = varlist[v]+piecewiselin; 2756 2776 2757 2777 /*Add input to element*/ … … 2759 2779 2760 2780 /*Clean-up*/ 2761 xDelete<IssmDouble>(monthlyeff_b); 2781 xDelete<IssmDouble>(datebreaks_b); 2782 xDelete<IssmDouble>(intercepts_b); 2783 xDelete<IssmDouble>(trends_b); 2762 2784 xDelete<IssmDouble>(varlist); 2763 }/*}}}*/ 2785 } 2786 /*}}}*/ 2764 2787 void Element::BeckmannGoosseFloatingiceMeltingRate(){/*{{{*/ 2765 2788 -
issm/trunk-jpl/src/c/classes/Elements/Element.h
r27318 r27334 164 164 void MigrateGroundingLine(IssmDouble* sheet_ungrounding); 165 165 void MismipFloatingiceMeltingRate(); 166 void Monthly EffectBasin(IssmDouble* monthlyeff,int enum_type);166 void MonthlyPiecewiseLinearEffectBasin(int nummonthbreaks,IssmDouble* monthlyintercepts,IssmDouble* monthlytrends,IssmDouble* monthlydatebreaks,int enum_type); 167 167 void BeckmannGoosseFloatingiceMeltingRate(); 168 168 void MungsmtpParameterization(void); -
issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp
r27320 r27334 50 50 bool isstochastic; 51 51 bool istfstochastic = false; 52 int M,N,arorder,maorder,numbasins,numparams,numbreaks, my_rank;52 int M,N,arorder,maorder,numbasins,numparams,numbreaks,nummonthbreaks,my_rank; 53 53 femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); 54 54 femmodel->parameters->FindParam(&numparams,FrontalForcingsNumberofParamsEnum); 55 55 femmodel->parameters->FindParam(&numbreaks,FrontalForcingsNumberofBreaksEnum); 56 femmodel->parameters->FindParam(&nummonthbreaks,FrontalForcingsNumberofMonthBreaksEnum); 56 57 femmodel->parameters->FindParam(&arorder,FrontalForcingsARMAarOrderEnum); 57 58 femmodel->parameters->FindParam(&maorder,FrontalForcingsARMAmaOrderEnum); 58 IssmDouble* datebreaks = NULL; 59 IssmDouble* arlagcoefs = NULL; 60 IssmDouble* malagcoefs = NULL; 61 IssmDouble* monthlyeff = NULL; 62 IssmDouble* polyparams = NULL; 59 IssmDouble* datebreaks = NULL; 60 IssmDouble* arlagcoefs = NULL; 61 IssmDouble* malagcoefs = NULL; 62 IssmDouble* monthlyeff = NULL; 63 IssmDouble* polyparams = NULL; 64 IssmDouble* monthdatebreaks = NULL; 65 IssmDouble* monthintercepts = NULL; 66 IssmDouble* monthtrends = NULL; 63 67 64 femmodel->parameters->FindParam(&datebreaks,&M,&N,FrontalForcingsARMAdatebreaksEnum); _assert_(M==numbasins); _assert_(N==max(numbreaks,1)); 65 femmodel->parameters->FindParam(&polyparams,&M,&N,FrontalForcingsARMApolyparamsEnum); _assert_(M==numbasins); _assert_(N==(numbreaks+1)*numparams); 66 femmodel->parameters->FindParam(&arlagcoefs,&M,&N,FrontalForcingsARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(N==arorder); 67 femmodel->parameters->FindParam(&malagcoefs,&M,&N,FrontalForcingsARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(N==maorder); 68 femmodel->parameters->FindParam(&monthlyeff,&M,&N,ThermalForcingMonthlyEffectsEnum); _assert_(M==numbasins); _assert_(N==12); 68 femmodel->parameters->FindParam(&datebreaks,&M,&N,FrontalForcingsARMAdatebreaksEnum); _assert_(M==numbasins); _assert_(N==max(numbreaks,1)); 69 femmodel->parameters->FindParam(&polyparams,&M,&N,FrontalForcingsARMApolyparamsEnum); _assert_(M==numbasins); _assert_(N==(numbreaks+1)*numparams); 70 femmodel->parameters->FindParam(&arlagcoefs,&M,&N,FrontalForcingsARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(N==arorder); 71 femmodel->parameters->FindParam(&malagcoefs,&M,&N,FrontalForcingsARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(N==maorder); 72 femmodel->parameters->FindParam(&monthdatebreaks,&M,&N,FrontalForcingsARMAmonthdatebreaksEnum); _assert_(M==numbasins); _assert_(N==max(numbreaks,1)); 73 femmodel->parameters->FindParam(&monthintercepts,&M,&N,FrontalForcingsARMAmonthinterceptsEnum); _assert_(M==numbasins); _assert_(N==12*(nummonthbreaks+1)); 74 femmodel->parameters->FindParam(&monthtrends,&M,&N,FrontalForcingsARMAmonthtrendsEnum); _assert_(M==numbasins); _assert_(N==12*(nummonthbreaks+1)); 69 75 70 76 femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); … … 86 92 element->ArmaProcess(isstepforarma,arorder,maorder,numparams,numbreaks,tstep_arma,polyparams,arlagcoefs,malagcoefs,datebreaks,istfstochastic,FrontalForcingsRignotarmaEnum); 87 93 /*Compute monthly effects*/ 88 element->Monthly EffectBasin(monthlyeff,FrontalForcingsRignotarmaEnum);94 element->MonthlyPiecewiseLinearEffectBasin(nummonthbreaks,monthintercepts,monthtrends,monthdatebreaks,FrontalForcingsRignotarmaEnum); 89 95 } 90 96 -
issm/trunk-jpl/src/c/shared/Enum/Enum.vim
r27318 r27334 174 174 syn keyword cConstant DebrisPackingFractionEnum 175 175 syn keyword cConstant DebugProfilingEnum 176 syn keyword cConstant DebrisThicknessEnum177 176 syn keyword cConstant DomainDimensionEnum 178 177 syn keyword cConstant DomainTypeEnum … … 218 217 syn keyword cConstant FrontalForcingsARMAmaOrderEnum 219 218 syn keyword cConstant FrontalForcingsARMAdatebreaksEnum 219 syn keyword cConstant FrontalForcingsARMAmonthdatebreaksEnum 220 syn keyword cConstant FrontalForcingsARMAmonthinterceptsEnum 221 syn keyword cConstant FrontalForcingsARMAmonthtrendsEnum 220 222 syn keyword cConstant FrontalForcingsARMApolyparamsEnum 221 223 syn keyword cConstant FrontalForcingsNumberofBasinsEnum 222 224 syn keyword cConstant FrontalForcingsNumberofBreaksEnum 225 syn keyword cConstant FrontalForcingsNumberofMonthBreaksEnum 223 226 syn keyword cConstant FrontalForcingsNumberofParamsEnum 224 227 syn keyword cConstant FrontalForcingsParamEnum … … 724 727 syn keyword cConstant DamageDbarOldEnum 725 728 syn keyword cConstant DamageFEnum 729 syn keyword cConstant DebrisThicknessEnum 726 730 syn keyword cConstant DegreeOfChannelizationEnum 727 731 syn keyword cConstant DepthBelowSurfaceEnum -
issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
r27323 r27334 211 211 FrontalForcingsARMAmaOrderEnum, 212 212 FrontalForcingsARMAdatebreaksEnum, 213 FrontalForcingsARMAmonthdatebreaksEnum, 214 FrontalForcingsARMAmonthinterceptsEnum, 215 FrontalForcingsARMAmonthtrendsEnum, 213 216 FrontalForcingsARMApolyparamsEnum, 214 217 FrontalForcingsNumberofBasinsEnum, 215 218 FrontalForcingsNumberofBreaksEnum, 219 FrontalForcingsNumberofMonthBreaksEnum, 216 220 FrontalForcingsNumberofParamsEnum, 217 221 FrontalForcingsParamEnum, -
issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
r27323 r27334 219 219 case FrontalForcingsARMAmaOrderEnum : return "FrontalForcingsARMAmaOrder"; 220 220 case FrontalForcingsARMAdatebreaksEnum : return "FrontalForcingsARMAdatebreaks"; 221 case FrontalForcingsARMAmonthdatebreaksEnum : return "FrontalForcingsARMAmonthdatebreaks"; 222 case FrontalForcingsARMAmonthinterceptsEnum : return "FrontalForcingsARMAmonthintercepts"; 223 case FrontalForcingsARMAmonthtrendsEnum : return "FrontalForcingsARMAmonthtrends"; 221 224 case FrontalForcingsARMApolyparamsEnum : return "FrontalForcingsARMApolyparams"; 222 225 case FrontalForcingsNumberofBasinsEnum : return "FrontalForcingsNumberofBasins"; 223 226 case FrontalForcingsNumberofBreaksEnum : return "FrontalForcingsNumberofBreaks"; 227 case FrontalForcingsNumberofMonthBreaksEnum : return "FrontalForcingsNumberofMonthBreaks"; 224 228 case FrontalForcingsNumberofParamsEnum : return "FrontalForcingsNumberofParams"; 225 229 case FrontalForcingsParamEnum : return "FrontalForcingsParam"; -
issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim
r27318 r27334 167 167 syn keyword juliaConstC DebrisPackingFractionEnum 168 168 syn keyword juliaConstC DebugProfilingEnum 169 syn keyword juliaConstC DebrisThicknessEnum170 169 syn keyword juliaConstC DomainDimensionEnum 171 170 syn keyword juliaConstC DomainTypeEnum … … 211 210 syn keyword juliaConstC FrontalForcingsARMAmaOrderEnum 212 211 syn keyword juliaConstC FrontalForcingsARMAdatebreaksEnum 212 syn keyword juliaConstC FrontalForcingsARMAmonthdatebreaksEnum 213 syn keyword juliaConstC FrontalForcingsARMAmonthinterceptsEnum 214 syn keyword juliaConstC FrontalForcingsARMAmonthtrendsEnum 213 215 syn keyword juliaConstC FrontalForcingsARMApolyparamsEnum 214 216 syn keyword juliaConstC FrontalForcingsNumberofBasinsEnum 215 217 syn keyword juliaConstC FrontalForcingsNumberofBreaksEnum 218 syn keyword juliaConstC FrontalForcingsNumberofMonthBreaksEnum 216 219 syn keyword juliaConstC FrontalForcingsNumberofParamsEnum 217 220 syn keyword juliaConstC FrontalForcingsParamEnum … … 717 720 syn keyword juliaConstC DamageDbarOldEnum 718 721 syn keyword juliaConstC DamageFEnum 722 syn keyword juliaConstC DebrisThicknessEnum 719 723 syn keyword juliaConstC DegreeOfChannelizationEnum 720 724 syn keyword juliaConstC DepthBelowSurfaceEnum -
issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
r27323 r27334 222 222 else if (strcmp(name,"FrontalForcingsARMAmaOrder")==0) return FrontalForcingsARMAmaOrderEnum; 223 223 else if (strcmp(name,"FrontalForcingsARMAdatebreaks")==0) return FrontalForcingsARMAdatebreaksEnum; 224 else if (strcmp(name,"FrontalForcingsARMAmonthdatebreaks")==0) return FrontalForcingsARMAmonthdatebreaksEnum; 225 else if (strcmp(name,"FrontalForcingsARMAmonthintercepts")==0) return FrontalForcingsARMAmonthinterceptsEnum; 226 else if (strcmp(name,"FrontalForcingsARMAmonthtrends")==0) return FrontalForcingsARMAmonthtrendsEnum; 224 227 else if (strcmp(name,"FrontalForcingsARMApolyparams")==0) return FrontalForcingsARMApolyparamsEnum; 225 228 else if (strcmp(name,"FrontalForcingsNumberofBasins")==0) return FrontalForcingsNumberofBasinsEnum; 226 229 else if (strcmp(name,"FrontalForcingsNumberofBreaks")==0) return FrontalForcingsNumberofBreaksEnum; 230 else if (strcmp(name,"FrontalForcingsNumberofMonthBreaks")==0) return FrontalForcingsNumberofMonthBreaksEnum; 227 231 else if (strcmp(name,"FrontalForcingsNumberofParams")==0) return FrontalForcingsNumberofParamsEnum; 228 232 else if (strcmp(name,"FrontalForcingsParam")==0) return FrontalForcingsParamEnum; … … 256 260 else if (strcmp(name,"HydrologydcEplMaxThickness")==0) return HydrologydcEplMaxThicknessEnum; 257 261 else if (strcmp(name,"HydrologydcEplPoreWaterMass")==0) return HydrologydcEplPoreWaterMassEnum; 258 else if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum; 262 else stage=3; 263 } 264 if(stage==3){ 265 if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum; 259 266 else if (strcmp(name,"HydrologydcEplflipLock")==0) return HydrologydcEplflipLockEnum; 260 267 else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum; 261 268 else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum; 262 else stage=3; 263 } 264 if(stage==3){ 265 if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum; 269 else if (strcmp(name,"HydrologydcMaxIter")==0) return HydrologydcMaxIterEnum; 266 270 else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum; 267 271 else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum; … … 379 383 else if (strcmp(name,"OceanGridNx")==0) return OceanGridNxEnum; 380 384 else if (strcmp(name,"OceanGridNy")==0) return OceanGridNyEnum; 381 else if (strcmp(name,"OceanGridX")==0) return OceanGridXEnum; 385 else stage=4; 386 } 387 if(stage==4){ 388 if (strcmp(name,"OceanGridX")==0) return OceanGridXEnum; 382 389 else if (strcmp(name,"OceanGridY")==0) return OceanGridYEnum; 383 390 else if (strcmp(name,"OutputBufferPointer")==0) return OutputBufferPointerEnum; 384 391 else if (strcmp(name,"OutputBufferSizePointer")==0) return OutputBufferSizePointerEnum; 385 else stage=4; 386 } 387 if(stage==4){ 388 if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; 392 else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum; 389 393 else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum; 390 394 else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum; … … 502 506 else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum; 503 507 else if (strcmp(name,"SmbASnow")==0) return SmbASnowEnum; 504 else if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum; 508 else stage=5; 509 } 510 if(stage==5){ 511 if (strcmp(name,"SmbAccualti")==0) return SmbAccualtiEnum; 505 512 else if (strcmp(name,"SmbAccugrad")==0) return SmbAccugradEnum; 506 513 else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum; 507 514 else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum; 508 else stage=5; 509 } 510 if(stage==5){ 511 if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum; 515 else if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum; 512 516 else if (strcmp(name,"SmbARMAarOrder")==0) return SmbARMAarOrderEnum; 513 517 else if (strcmp(name,"SmbARMAmaOrder")==0) return SmbARMAmaOrderEnum; … … 625 629 else if (strcmp(name,"TransientAmrFrequency")==0) return TransientAmrFrequencyEnum; 626 630 else if (strcmp(name,"TransientIsage")==0) return TransientIsageEnum; 627 else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum; 631 else stage=6; 632 } 633 if(stage==6){ 634 if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum; 628 635 else if (strcmp(name,"TransientIsdebris")==0) return TransientIsdebrisEnum; 629 636 else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum; 630 637 else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; 631 else stage=6; 632 } 633 if(stage==6){ 634 if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; 638 else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; 635 639 else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; 636 640 else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; … … 748 752 else if (strcmp(name,"DeltaDsl")==0) return DeltaDslEnum; 749 753 else if (strcmp(name,"DslOld")==0) return DslOldEnum; 750 else if (strcmp(name,"Dsl")==0) return DslEnum; 754 else stage=7; 755 } 756 if(stage==7){ 757 if (strcmp(name,"Dsl")==0) return DslEnum; 751 758 else if (strcmp(name,"DeltaStr")==0) return DeltaStrEnum; 752 759 else if (strcmp(name,"StrOld")==0) return StrOldEnum; 753 760 else if (strcmp(name,"Str")==0) return StrEnum; 754 else stage=7; 755 } 756 if(stage==7){ 757 if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; 761 else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum; 758 762 else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum; 759 763 else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum; … … 871 875 else if (strcmp(name,"LevelsetObservation")==0) return LevelsetObservationEnum; 872 876 else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum; 873 else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum; 877 else stage=8; 878 } 879 if(stage==8){ 880 if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum; 874 881 else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum; 875 882 else if (strcmp(name,"MaskOceanLevelset")==0) return MaskOceanLevelsetEnum; 876 883 else if (strcmp(name,"MaskIceLevelset")==0) return MaskIceLevelsetEnum; 877 else stage=8; 878 } 879 if(stage==8){ 880 if (strcmp(name,"MaskIceRefLevelset")==0) return MaskIceRefLevelsetEnum; 884 else if (strcmp(name,"MaskIceRefLevelset")==0) return MaskIceRefLevelsetEnum; 881 885 else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum; 882 886 else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum; … … 994 998 else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum; 995 999 else if (strcmp(name,"SigmaVM")==0) return SigmaVMEnum; 996 else if (strcmp(name,"SmbAccumulatedEC")==0) return SmbAccumulatedECEnum; 1000 else stage=9; 1001 } 1002 if(stage==9){ 1003 if (strcmp(name,"SmbAccumulatedEC")==0) return SmbAccumulatedECEnum; 997 1004 else if (strcmp(name,"SmbAccumulatedMassBalance")==0) return SmbAccumulatedMassBalanceEnum; 998 1005 else if (strcmp(name,"SmbAccumulatedMelt")==0) return SmbAccumulatedMeltEnum; 999 1006 else if (strcmp(name,"SmbAccumulatedPrecipitation")==0) return SmbAccumulatedPrecipitationEnum; 1000 else stage=9; 1001 } 1002 if(stage==9){ 1003 if (strcmp(name,"SmbAccumulatedRain")==0) return SmbAccumulatedRainEnum; 1007 else if (strcmp(name,"SmbAccumulatedRain")==0) return SmbAccumulatedRainEnum; 1004 1008 else if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum; 1005 1009 else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum; … … 1117 1121 else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum; 1118 1122 else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum; 1119 else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; 1123 else stage=10; 1124 } 1125 if(stage==10){ 1126 if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; 1120 1127 else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; 1121 1128 else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; 1122 1129 else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; 1123 else stage=10; 1124 } 1125 if(stage==10){ 1126 if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; 1130 else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; 1127 1131 else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum; 1128 1132 else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum; … … 1240 1244 else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum; 1241 1245 else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; 1242 else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; 1246 else stage=11; 1247 } 1248 if(stage==11){ 1249 if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; 1243 1250 else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; 1244 1251 else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; 1245 1252 else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; 1246 else stage=11; 1247 } 1248 if(stage==11){ 1249 if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; 1253 else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; 1250 1254 else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; 1251 1255 else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum; … … 1363 1367 else if (strcmp(name,"ControlInputValues")==0) return ControlInputValuesEnum; 1364 1368 else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum; 1365 else if (strcmp(name,"Cuffey")==0) return CuffeyEnum; 1369 else stage=12; 1370 } 1371 if(stage==12){ 1372 if (strcmp(name,"Cuffey")==0) return CuffeyEnum; 1366 1373 else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum; 1367 1374 else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum; 1368 1375 else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum; 1369 else stage=12; 1370 } 1371 if(stage==12){ 1372 if (strcmp(name,"DataSet")==0) return DataSetEnum; 1376 else if (strcmp(name,"DataSet")==0) return DataSetEnum; 1373 1377 else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum; 1374 1378 else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; … … 1486 1490 else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; 1487 1491 else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; 1488 else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; 1492 else stage=13; 1493 } 1494 if(stage==13){ 1495 if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; 1489 1496 else if (strcmp(name,"LinearFloatingMeltRatearma")==0) return LinearFloatingMeltRatearmaEnum; 1490 1497 else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; 1491 1498 else if (strcmp(name,"Loads")==0) return LoadsEnum; 1492 else stage=13; 1493 } 1494 if(stage==13){ 1495 if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; 1499 else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum; 1496 1500 else if (strcmp(name,"LoveHf")==0) return LoveHfEnum; 1497 1501 else if (strcmp(name,"LoveHt")==0) return LoveHtEnum; … … 1609 1613 else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; 1610 1614 else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; 1611 else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; 1615 else stage=14; 1616 } 1617 if(stage==14){ 1618 if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; 1612 1619 else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; 1613 1620 else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; 1614 1621 else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; 1615 else stage=14; 1616 } 1617 if(stage==14){ 1618 if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; 1622 else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; 1619 1623 else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; 1620 1624 else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; -
issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m
r27318 r27334 6 6 classdef frontalforcingsrignotarma 7 7 properties (SetAccess=public) 8 num_basins = 0; 9 num_params = 0; 10 num_breaks = 0; 11 polynomialparams = NaN; 12 datebreaks = NaN; 13 ar_order = 0; 14 ma_order = 0; 15 arma_timestep = 0; 16 arlag_coefs = NaN; 17 malag_coefs = NaN; 18 monthly_effects = NaN; 19 basin_id = NaN; 20 subglacial_discharge = NaN; 8 num_basins = 0; 9 num_params = 0; 10 num_breaks = 0; 11 polynomialparams = NaN; 12 datebreaks = NaN; 13 ar_order = 0; 14 ma_order = 0; 15 arma_timestep = 0; 16 arlag_coefs = NaN; 17 malag_coefs = NaN; 18 monthlyvals_intercepts = NaN; 19 monthlyvals_trends = NaN; 20 monthlyvals_numbreaks = 0; 21 monthlyvals_datebreaks = NaN; 22 basin_id = NaN; 23 subglacial_discharge = NaN; 21 24 end 22 25 methods … … 52 55 end % }}} 53 56 function md = checkconsistency(self,md,solution,analyses) % {{{ 54 %Early return 55 if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end 56 57 nbas = md.frontalforcings.num_basins; 58 nprm = md.frontalforcings.num_params; 59 nbrk = md.frontalforcings.num_breaks; 60 md = checkfield(md,'fieldname','frontalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); 61 md = checkfield(md,'fieldname','frontalforcings.num_params','numel',1,'NaN',1,'Inf',1,'>',0); 62 md = checkfield(md,'fieldname','frontalforcings.num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0); 57 %Early return 58 if (~strcmp(solution,'TransientSolution') | md.transient.ismovingfront==0), return; end 59 60 nbas = md.frontalforcings.num_basins; 61 nprm = md.frontalforcings.num_params; 62 nbrk = md.frontalforcings.num_breaks; 63 nMbrk = md.frontalforcings.monthlyvals_numbreaks; 64 md = checkfield(md,'fieldname','frontalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); 65 md = checkfield(md,'fieldname','frontalforcings.num_params','numel',1,'NaN',1,'Inf',1,'>',0); 66 md = checkfield(md,'fieldname','frontalforcings.num_breaks','numel',1,'NaN',1,'Inf',1,'>=',0); 63 67 md = checkfield(md,'fieldname','frontalforcings.basin_id','Inf',1,'>=',0,'<=',md.frontalforcings.num_basins,'size',[md.mesh.numberofelements 1]); 64 65 66 md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm); 67 68 69 70 71 72 73 74 68 md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'NaN',1,'Inf',1,'timeseries',1); 69 if(nbas>1 && nbrk>=1 && nprm>1) 70 md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1,nprm],'numel',nbas*(nbrk+1)*nprm); 71 elseif(nbas==1) 72 md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nprm,nbrk+1],'numel',nbas*(nbrk+1)*nprm); 73 elseif(nbrk==0) 74 md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nprm],'numel',nbas*(nbrk+1)*nprm); 75 elseif(nprm==1) 76 md = checkfield(md,'fieldname','frontalforcings.polynomialparams','NaN',1,'Inf',1,'size',[nbas,nbrk+1],'numel',nbas*(nbrk+1)*nprm); 77 end 78 md = checkfield(md,'fieldname','frontalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); 75 79 md = checkfield(md,'fieldname','frontalforcings.ma_order','numel',1,'NaN',1,'Inf',1,'>=',0); 76 80 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 77 md = checkfield(md,'fieldname','frontalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]); 78 md = checkfield(md,'fieldname','frontalforcings.malag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ma_order]); 79 if(nbrk>0) 80 md = checkfield(md,'fieldname','frontalforcings.datebreaks','NaN',1,'Inf',1,'size',[nbas,nbrk]); 81 elseif(numel(md.frontalforcings.datebreaks)==0 || all(isnan(md.frontalforcings.datebreaks))) 82 ; 83 else 84 error('md.frontalforcings.num_breaks is 0 but md.frontalforcings.datebreaks is not empty'); 85 end 86 if(any(~isnan(md.frontalforcings.monthly_effects))) 87 md = checkfield(md,'fieldname','frontalforcings.monthly_effects','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,12]); 88 if(md.timestepping.time_step>=1) 89 error('md.frontalforcings.monthly_effects are provided but md.timestepping.time_step>=1'); 90 end 91 end 92 93 end % }}} 81 md = checkfield(md,'fieldname','frontalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]); 82 md = checkfield(md,'fieldname','frontalforcings.malag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ma_order]); 83 if(nbrk>0) 84 md = checkfield(md,'fieldname','frontalforcings.datebreaks','NaN',1,'Inf',1,'size',[nbas,nbrk]); 85 elseif(numel(md.frontalforcings.datebreaks)==0 || all(isnan(md.frontalforcings.datebreaks))) 86 ; 87 else 88 error('md.frontalforcings.num_breaks is 0 but md.frontalforcings.datebreaks is not empty'); 89 end 90 91 %%% Check if some monthly forcings are provided %%% 92 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))) 93 isMonthly = true; 94 else 95 isMonthly = false; 96 end 97 if(numel(md.frontalforcings.monthlyvals_trends)>1 || ~isnan(md.frontalforcings.monthlyvals_trends)) 98 isMonthlyTrend = true; 99 else 100 isMonthlyTrend = false; 101 end 102 if(isMonthly) 103 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_numbreaks','numel',1,'NaN',1,'Inf',1,'>=',0); 104 if(nbas>1 && nMbrk>=1) 105 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nbas,12,nMbrk+1],'numel',nbas*(nMbrk+1)*12); 106 if(isMonthlyTrend) 107 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nbas,12,nMbrk+1],'numel',nbas*(nMbrk+1)*12); 108 end 109 elseif(nbas==1) 110 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nMbrk+1,12],'numel',nbas*(nMbrk+1)*12); 111 if(isMonthlyTrend) 112 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nMbrk+1,12],'numel',nbas*(nMbrk+1)*12); 113 end 114 elseif(nMbrk==0) 115 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nbas,12],'numel',nbas*(nMbrk+1)*12); 116 if(isMonthlyTrend) 117 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nbas,12],'numel',nbas*(nMbrk+1)*12); 118 end 119 end 120 if(nMbrk>0) 121 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_datebreaks','NaN',1,'Inf',1,'size',[nbas,nMbrk]); 122 elseif(numel(md.frontalforcings.monthlyvals_datebreaks)==0 || all(isnan(md.frontalforcings.monthlyvals_datebreaks))) 123 ; 124 else 125 error('md.frontalforcings.monthlyvals_numbreaks is 0 but md.frontalforcings.monthlyvals_datebreaks is not empty'); 126 end 127 end 128 129 end % }}} 94 130 function disp(self) % {{{ 95 131 disp(sprintf(' Frontalforcings parameters:')); … … 107 143 fielddisplay(self,'arlag_coefs','basin-specific vectors of AR lag coefficients [unitless]'); 108 144 fielddisplay(self,'malag_coefs','basin-specific vectors of MA lag coefficients [unitless]'); 109 fielddisplay(self,'monthly_effects','basin-specific monthly values of TF added at corresponding month (default: all 0) [°C]'); 145 fielddisplay(self,'monthlyvals_intercepts','intercept of basin-specific piecewise-linear monthly values of TF added at corresponding month (default: all 0) [°C]'); 146 fielddisplay(self,'monthlyvals_trends','trends in basin-specific piecewise-linear monthly values of TF added at corresponding month (default: all 0) [°C/yr]'); 147 fielddisplay(self,'monthlyvals_numbreaks','number of breakpoints in the piecewise-linear functions of monthly values'); 148 fielddisplay(self,'monthlyvals_datebreaks','dates at which the breakpoints in the piecewise-linear monthly values occur (1 row per basin)'); 110 149 end % }}} 111 150 function marshall(self,prefix,md,fid) % {{{ 112 151 yts=md.constants.yts; 113 nbas = md.frontalforcings.num_basins; 114 nprm = md.frontalforcings.num_params; 115 nper = md.frontalforcings.num_breaks+1; 152 %%% Deal with polynomial %%% 153 nbas = md.frontalforcings.num_basins; 154 nprm = md.frontalforcings.num_params; 155 nper = md.frontalforcings.num_breaks+1; 116 156 % Scale the parameters % 117 157 polyparamsScaled = md.frontalforcings.polynomialparams; … … 150 190 polyparams2dScaled = polyparamsScaled; 151 191 end 152 if(any(isnan(md.frontalforcings.monthly_effects))) %monthly effects not provided, set to 0153 meffects = zeros(md.frontalforcings.num_basins,12);154 else155 meffects = md.frontalforcings.monthly_effects;156 end157 192 if(nper==1) %a single period (no break date) 158 193 dbreaks = zeros(nbas,1); %dummy 159 194 else 160 195 dbreaks = md.frontalforcings.datebreaks; 196 end 197 198 %%% Deal with monthly effects %%% 199 nMper = md.frontalforcings.monthlyvals_numbreaks+1; 200 if((numel(md.frontalforcings.monthlyvals_intercepts)<=1 && isnan(md.frontalforcings.monthlyvals_intercepts))) 201 interceptsM = zeros(nbas,12); %monthly intercepts not provided, set to 0 202 trendsM = zeros(nbas,12); %set monthly trends also to 0 203 else 204 interceptsM3d = md.frontalforcings.monthlyvals_intercepts; 205 if((numel(md.frontalforcings.monthlyvals_trends)<=1 && isnan(md.frontalforcings.monthlyvals_trends))) 206 trendsM3d = 0*interceptsM3d; %monthly trends not provided, set to 0 207 else 208 trendsM3d = md.frontalforcings.monthlyvals_trends; 209 end 210 end 211 % Create 2D arrays from 3D arrays if needed % 212 if(nMper>1 && (numel(md.frontalforcings.monthlyvals_intercepts)>1 || ~isnan(md.frontalforcings.monthlyvals_intercepts))) 213 interceptsM = zeros(nbas,12*nMper); 214 trendsM = zeros(nbas,12*nMper); 215 for(ii=[1:nMper]) 216 jj = 1+(ii-1)*12; 217 interceptsM(:,jj:jj+12-1) = interceptsM3d(:,:,ii); 218 trendsM(:,jj:jj+12-1) = trendsM3d(:,:,ii); 219 end 220 elseif(nMper==1 && (numel(md.frontalforcings.monthlyvals_intercepts)>1 || ~isnan(md.frontalforcings.monthlyvals_intercepts))) 221 interceptsM = interceptsM3d; 222 trendsM = trendsM3d; 223 end 224 if(nMper==1) %a single period (no break date) 225 dMbreaks = zeros(nbas,1); %dummy 226 else 227 dMbreaks = md.frontalforcings.monthlyvals_datebreaks; 161 228 end 162 229 … … 174 241 WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','malag_coefs','format','DoubleMat','name','md.frontalforcings.malag_coefs','yts',yts); 175 242 WriteData(fid,prefix,'data',dbreaks,'name','md.frontalforcings.datebreaks','format','DoubleMat','scale',yts); 176 WriteData(fid,prefix,'data',meffects,'name','md.frontalforcings.monthly_effects','format','DoubleMat'); 243 WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','monthlyvals_numbreaks','format','Integer'); 244 WriteData(fid,prefix,'data',dMbreaks,'name','md.frontalforcings.monthlyvals_datebreaks','format','DoubleMat','scale',yts); 245 WriteData(fid,prefix,'data',interceptsM,'name','md.frontalforcings.monthlyvals_intercepts','format','DoubleMat'); 246 WriteData(fid,prefix,'data',trendsM,'name','md.frontalforcings.monthlyvals_trends','format','DoubleMat','scale',1/yts); 177 247 end % }}} 178 248 end -
issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py
r27318 r27334 25 25 self.arlag_coefs = np.nan 26 26 self.malag_coefs = np.nan 27 self.monthly_effects = np.nan 27 self.monthlyvals_intercepts = np.nan 28 self.monthlyvals_trends = np.nan 29 self.monthlyvals_numbreaks = 0 30 self.monthlyvals_datebreaks = np.nan 28 31 self.basin_id = np.nan 29 32 self.subglacial_discharge = np.nan … … 48 51 s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of AR lag coefficients [unitless]')) 49 52 s += '{}\n'.format(fielddisplay(self, 'malag_coefs', 'basin-specific vectors of MA lag coefficients [unitless]')) 50 s += '{}\n'.format(fielddisplay(self, 'monthly_effects', 'basin-specific monthly values of TF added at corresponding month (default: all 0) [°C]'))51 53 return s 52 54 #}}} … … 66 68 return md 67 69 68 nbas = md.frontalforcings.num_basins; 69 nprm = md.frontalforcings.num_params; 70 nbrk = md.frontalforcings.num_breaks; 70 nbas = md.frontalforcings.num_basins; 71 nprm = md.frontalforcings.num_params; 72 nbrk = md.frontalforcings.num_breaks; 73 nMbrk = md.frontalforcings.monthlyvals_numbreaks; 71 74 md = checkfield(md, 'fieldname', 'frontalforcings.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) 72 75 md = checkfield(md, 'fieldname', 'frontalforcings.num_params', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) … … 95 98 else: 96 99 raise RuntimeError('md.frontalforcings.num_breaks is 0 but md.frontalforcings.datebreaks is not empty') 97 if(np.any(np.isnan(md.frontalforcings.monthly_effects)==False)): 98 md = checkfield(md, 'fieldname', 'frontalforcings.monthly_effects', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, 12]) 99 if(md.timestepping.time_step>=1): 100 raise RuntimeError('md.frontalforcings.monthly_effects are provided but md.timestepping.time_step>=1') 100 101 102 ### Check if some monthly forcings are provided ### 103 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): 104 isMonthly = True 105 else: 106 isMonthly = False 107 if(np.all(np.isnan(md.frontalforcings.monthlyvals_datebreaks))): 108 isMonthlyTrend = True 109 else: 110 isMonthlyTrend = False 111 if(isMonthly): 112 md = checkfield(md, 'fieldname', 'frontalforcings.monthlyvals_numbreaks', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) 113 if(nbas>1 and nMbrk>=1): 114 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nbas,12,nMbrk+1],'numel',nbas*(nMbrk+1)*12) 115 if(isMonthlyTrend): 116 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nbas,12,nMbrk+1],'numel',nbas*(nMbrk+1)*12) 117 elif(nbas==1): 118 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nMbrk+1,12],'numel',nbas*(nMbrk+1)*12) 119 if(isMonthlyTrend): 120 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nMbrk+1,12],'numel',nbas*(nMbrk+1)*12) 121 elif(nMbrk==0): 122 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_intercepts','NaN',1,'Inf',1,'size',[nbas,12],'numel',nbas*(nMbrk+1)*12) 123 if(isMonthlyTrend): 124 md = checkfield(md,'fieldname','frontalforcings.monthlyvals_trends','NaN',1,'Inf',1,'size',[nbas,12],'numel',nbas*(nMbrk+1)*12) 125 if(nMbrk>0): 126 md = checkfield(md, 'fieldname', 'frontalforcings.monthlyvals_datebreaks', 'NaN', 1, 'Inf', 1, 'size', [nbas,nMbrk]) 127 elif(np.size(md.frontalforcings.monthlyvals_datebreaks)==0 or np.all(np.isnan(md.frontalforcings.monthlyvals_datebreaks))): 128 pass 129 else: 130 raise RuntimeError('md.frontalforcings.monthlyvals_numbreaks is 0 but md.frontalforcings.monthlyvals_datebreaks is not empty') 131 101 132 return md 102 133 # }}} … … 139 170 # 2D array is already in correct format and no need for scaling # 140 171 polyparams2dScaled = np.copy(polyparamsScaled) 141 if(np.any(np.isnan(md.frontalforcings.monthly_effects))): #monthly effects not provided, set to 0142 meffects = np.zeros((md.frontalforcings.num_basins,12))143 else:144 meffects = 1*md.frontalforcings.monthly_effects145 172 if(nper==1): 146 173 dbreaks = np.zeros((nbas,1)) 147 174 else: 148 175 dbreaks = np.copy(md.frontalforcings.datebreaks) 176 177 ### Deal with montly effects ### 178 nMper = md.frontalforcings.monthlyvals_numbreaks+1 179 if(np.any(np.isnan(md.frontalforcings.monthlyvals_intercepts))): 180 interceptsM = np.zeros((nbas,12)) #monthly intercepts not provided, set to 0 181 trendsM = np.zeros((nbas,12)) #set monthly trends also to 0 182 else: 183 interceptsM3d = md.frontalforcings.monthlyvals_intercepts 184 if(np.any(np.isnan(md.frontalforcings.monthlyvals_trends))): 185 trendsM3d = 0*interceptsM3d #monthly trends not provided, set to 0 186 else: 187 trendsM3d = md.frontalforcings.monthlyvals_trends 188 # Create 2D arrays from 3D arrays if needed # 189 if(nMper>1 and np.all(np.isnan(md.frontalforcings.monthlyvals_intercepts))==False): 190 interceptsM = np.zeros((nbas,12*nMper)) 191 trendsM = np.zeros((nbas,12*nMper)) 192 for ii in range(nMper): 193 interceptsM[:,ii*12:(ii+1)*12] = 1*interceptsM3d[:,:,ii] 194 trendsM[:,ii*12:(ii+1)*12] = 1*trendsM3d[:,:,ii] 195 elif(nMper==1 and np.all(np.isnan(md.frontalforcings.monthlyvals_intercepts))==False): 196 interceptsM = 1*interceptsM3d 197 trendsM = 1*trendsM3d 198 if(nMper==1): 199 dMbreaks = np.zeros((nbas,1)) 200 else: 201 dMbreaks = np.copy(md.frontalforcings.monthlyvals_datebreaks) 149 202 150 203 WriteData(fid, prefix, 'name', 'md.frontalforcings.parameterization', 'data', 3, 'format', 'Integer') … … 161 214 WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'malag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.malag_coefs', 'yts', yts) 162 215 WriteData(fid, prefix, 'data', dbreaks, 'name', 'md.frontalforcings.datebreaks', 'format', 'DoubleMat','scale',yts) 163 WriteData(fid, prefix, 'data', meffects, 'name', 'md.frontalforcings.monthly_effects', 'format', 'DoubleMat') 216 WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','monthlyvals_numbreaks','format','Integer') 217 WriteData(fid,prefix,'data',dMbreaks,'name','md.frontalforcings.monthlyvals_datebreaks','format','DoubleMat','scale',yts) 218 WriteData(fid,prefix,'data',interceptsM,'name','md.frontalforcings.monthlyvals_intercepts','format','DoubleMat') 219 WriteData(fid,prefix,'data',trendsM,'name','md.frontalforcings.monthlyvals_trends','format','DoubleMat','scale',1/yts) 164 220 # }}} 221 222 -
issm/trunk-jpl/test/NightlyRun/test543.m
r27318 r27334 5 5 md=setflowequation(md,'SSA','all'); 6 6 md.timestepping.start_time = 0; 7 md.timestepping.time_step = 1;8 md.timestepping.final_time = 10;7 md.timestepping.time_step = 0.05; 8 md.timestepping.final_time = 2; 9 9 10 10 %Basin separation TF … … 52 52 trendlin = [-0.5,-0.2,0.1;0,0,0]; 53 53 trendquad = [0,0.0,0;0.1,0.1,0.1]; 54 datebreaks = [ 4,7;4,7];54 datebreaks = [0.5,1.5;0.5,1.5]; 55 55 polynomialparams = cat(numparams,intercept,trendlin,trendquad); 56 % Monthly effects params % 57 numbreaksM = 1; 58 intcpsMp0 = [-0.5,-0.5,0,0,0,0,0.5,0.5,0,0,0,0; 59 -1.0,-1.0,0,0,0,0,1.0,1.0,0,0,0,0]; 60 intcpsMp1 = [-0.25,-0.25,0,0,0,0,0.,0.,0,0,0,0; 61 -0.1,-0.1,0,0,0,0,0.1,0.1,0,0,0,0]; 62 intcpsM = cat(3,intcpsMp0,intcpsMp1); 63 trendsMp0 = [0,0,0,0,0,0,0.,0.0,0,0,0,0; 64 0.0,0.0,0,-0.0,0,0,0.0,0.0,0,0,0,0]; 65 trendsMp1 = [0,-0.12,0,0,0,0,0.,0.0,0,0.0,0,0; 66 0.0,-0.1,0,-0.0,0,0,0.0,0.0,0,0,0,0]; 67 trendsM = cat(3,trendsMp0,trendsMp1); 68 datebreaksM = [1;1]; 56 69 57 md.frontalforcings.num_basins = nb_tf; 58 md.frontalforcings.basin_id = idb_tf; 59 md.frontalforcings.num_params = numparams; %number of parameters in the polynomial 60 md.frontalforcings.num_breaks = numbreaks; %number of breakpoints 61 md.frontalforcings.subglacial_discharge = 0.01*ones(md.mesh.numberofvertices,1); 62 md.frontalforcings.polynomialparams = polynomialparams; 63 md.frontalforcings.datebreaks = datebreaks; 64 md.frontalforcings.ar_order = 4; 65 md.frontalforcings.ma_order = 2; 66 md.frontalforcings.arma_timestep = 2; %timestep of the ARMA model [yr] 67 md.frontalforcings.arlag_coefs = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters 68 md.frontalforcings.malag_coefs = [[0.1,0.0];[0.0,0.1]]; %moving-average parameters 69 70 md.frontalforcings.num_basins = nb_tf; 71 md.frontalforcings.basin_id = idb_tf; 72 md.frontalforcings.num_params = numparams; %number of parameters in the polynomial 73 md.frontalforcings.num_breaks = numbreaks; %number of breakpoints 74 md.frontalforcings.subglacial_discharge = 0.01*ones(md.mesh.numberofvertices,1); 75 md.frontalforcings.polynomialparams = polynomialparams; 76 md.frontalforcings.datebreaks = datebreaks; 77 md.frontalforcings.ar_order = 4; 78 md.frontalforcings.ma_order = 2; 79 md.frontalforcings.arma_timestep = 2; %timestep of the ARMA model [yr] 80 md.frontalforcings.arlag_coefs = [[0.1,-0.1,0.01,-0.01];[0.2,-0.2,0.1,0.0]]; %autoregressive parameters 81 md.frontalforcings.malag_coefs = [[0.1,0.0];[0.0,0.1]]; %moving-average parameters 82 md.frontalforcings.monthlyvals_numbreaks = numbreaksM; 83 md.frontalforcings.monthlyvals_intercepts = intcpsM; 84 md.frontalforcings.monthlyvals_trends = trendsM; 85 md.frontalforcings.monthlyvals_datebreaks = datebreaksM; 70 86 % Floating Ice Melt parameters 71 87 md.basalforcings.floatingice_melting_rate = 0.1*ones(md.mesh.numberofvertices,1); … … 102 118 1e-11,2e-11,2e-11,1e-11,1e-9,1e-10,1e-10,1e-10,... 103 119 2e-11,1e-11,1e-11,9e-11,2e-9,1e-10,1e-10,1e-10,... 104 2e- 10,1e-10,1e-10,1e-10,5e-9,1e-10,1e-10,1e-10,...120 2e-6,1e-6,1e-6,1e-6,5e-6,1e-6,1e-6,1e-6,... 105 121 }; 106 122 field_values={... … … 113 129 (md.results.TransientSolution(1).CalvingMeltingrate),... 114 130 (md.results.TransientSolution(1).BasalforcingsFloatingiceMeltingRate),... 115 (md.results.TransientSolution( 5).Vx),...116 (md.results.TransientSolution( 5).Vy),...117 (md.results.TransientSolution( 5).Vel),...118 (md.results.TransientSolution( 5).Thickness),...119 (md.results.TransientSolution( 5).MaskIceLevelset),...120 (md.results.TransientSolution( 5).CalvingCalvingrate),...121 (md.results.TransientSolution( 5).CalvingMeltingrate),...122 (md.results.TransientSolution( 5).BasalforcingsFloatingiceMeltingRate),...123 (md.results.TransientSolution( 10).Vx),...124 (md.results.TransientSolution( 10).Vy),...125 (md.results.TransientSolution( 10).Vel),...126 (md.results.TransientSolution( 10).Thickness),...127 (md.results.TransientSolution( 10).MaskIceLevelset),...128 (md.results.TransientSolution( 10).CalvingCalvingrate),...129 (md.results.TransientSolution( 10).CalvingMeltingrate),...130 (md.results.TransientSolution( 10).BasalforcingsFloatingiceMeltingRate),...131 (md.results.TransientSolution(20).Vx),... 132 (md.results.TransientSolution(20).Vy),... 133 (md.results.TransientSolution(20).Vel),... 134 (md.results.TransientSolution(20).Thickness),... 135 (md.results.TransientSolution(20).MaskIceLevelset),... 136 (md.results.TransientSolution(20).CalvingCalvingrate),... 137 (md.results.TransientSolution(20).CalvingMeltingrate),... 138 (md.results.TransientSolution(20).BasalforcingsFloatingiceMeltingRate),... 139 (md.results.TransientSolution(40).Vx),... 140 (md.results.TransientSolution(40).Vy),... 141 (md.results.TransientSolution(40).Vel),... 142 (md.results.TransientSolution(40).Thickness),... 143 (md.results.TransientSolution(40).MaskIceLevelset),... 144 (md.results.TransientSolution(40).CalvingCalvingrate),... 145 (md.results.TransientSolution(40).CalvingMeltingrate),... 146 (md.results.TransientSolution(40).BasalforcingsFloatingiceMeltingRate),... 131 147 }; -
issm/trunk-jpl/test/NightlyRun/test543.py
r27318 r27334 16 16 md = setflowequation(md, 'SSA', 'all') 17 17 md.timestepping.start_time = 0 18 md.timestepping.time_step = 119 md.timestepping.final_time = 1018 md.timestepping.time_step = 0.05 19 md.timestepping.final_time = 2 20 20 21 21 # Basin separation TF … … 57 57 trendlin = np.array([[-0.5,-0.2,0.1],[0,0,0]]) 58 58 trendquad = np.array([[0,0.0,0],[0.1,0.1,0.1]]) 59 datebreaks = np.array([[ 4,7],[4,7]])59 datebreaks = np.array([[0.5,1.5],[0.5,1.5]]) 60 60 polynomialparams = np.transpose(np.stack((intercept,trendlin,trendquad)),(1,2,0)) 61 # Monthly effects params # 62 numbreaksM = 1; 63 intcpsMp0 = np.array([[-0.5,-0.5,0,0,0,0,0.5,0.5,0,0,0,0], 64 [-1.0,-1.0,0,0,0,0,1.0,1.0,0,0,0,0]]) 65 intcpsMp1 = np.array([[-0.25,-0.25,0,0,0,0,0.,0.,0,0,0,0], 66 [-0.1,-0.1,0,0,0,0,0.1,0.1,0,0,0,0]]) 67 intcpsM = np.transpose(np.stack((intcpsMp0,intcpsMp1)),(1,2,0)) 68 trendsMp0 = np.array([[0,0,0,0,0,0,0.,0.0,0,0,0,0], 69 [0.0,0.0,0,-0.0,0,0,0.0,0.0,0,0,0,0]]) 70 trendsMp1 = np.array([[0,-0.12,0,0,0,0,0.,0.0,0,0.0,0,0], 71 [0.0,-0.1,0,-0.0,0,0,0.0,0.0,0,0,0,0]]) 72 trendsM = np.transpose(np.stack((trendsMp0,trendsMp1)),(1,2,0)) 73 datebreaksM = np.array([[1],[1]]) 61 74 62 75 md.frontalforcings.num_basins = nb_tf … … 72 85 md.frontalforcings.arlag_coefs = np.array([[0.1, -0.1, 0.01, -0.01], [0.2, -0.2, 0.1, 0.0]]) # autoregressive parameters 73 86 md.frontalforcings.malag_coefs = np.array([[0.1, 0.0], [0.0, 0.1]]) # moving-average parameters 87 md.frontalforcings.monthlyvals_numbreaks = numbreaksM 88 md.frontalforcings.monthlyvals_intercepts = intcpsM 89 md.frontalforcings.monthlyvals_trends = trendsM 90 md.frontalforcings.monthlyvals_datebreaks = datebreaksM 74 91 #Floating Ice Melt parameters 75 92 md.basalforcings.floatingice_melting_rate = 0.1 * np.ones((md.mesh.numberofvertices,)) … … 119 136 1e-11, 2e-11, 2e-11, 1e-11, 1e-9, 1e-10, 1e-10, 1e-10, 120 137 2e-11, 1e-11, 1e-11, 9e-11, 2e-9, 1e-10, 1e-10, 1e-10, 121 2e- 10, 1e-10, 1e-10, 1e-10, 5e-9, 1e-10, 1e-10, 1e-10]138 2e-6, 1e-6, 1e-6, 1e-6, 5e-6, 1e-6, 1e-6, 1e-6] 122 139 field_values = [ 123 140 md.results.TransientSolution[0].Vx, … … 129 146 md.results.TransientSolution[0].CalvingMeltingrate, 130 147 md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate, 131 md.results.TransientSolution[ 4].Vx,132 md.results.TransientSolution[ 4].Vy,133 md.results.TransientSolution[ 4].Vel,134 md.results.TransientSolution[ 4].Thickness,135 md.results.TransientSolution[ 4].MaskIceLevelset,136 md.results.TransientSolution[ 4].CalvingCalvingrate,137 md.results.TransientSolution[ 4].CalvingMeltingrate,138 md.results.TransientSolution[ 4].BasalforcingsFloatingiceMeltingRate,139 md.results.TransientSolution[ 9].Vx,140 md.results.TransientSolution[ 9].Vy,141 md.results.TransientSolution[ 9].Vel,142 md.results.TransientSolution[ 9].Thickness,143 md.results.TransientSolution[ 9].MaskIceLevelset,144 md.results.TransientSolution[ 9].CalvingCalvingrate,145 md.results.TransientSolution[ 9].CalvingMeltingrate,146 md.results.TransientSolution[ 9].BasalforcingsFloatingiceMeltingRate]148 md.results.TransientSolution[19].Vx, 149 md.results.TransientSolution[19].Vy, 150 md.results.TransientSolution[19].Vel, 151 md.results.TransientSolution[19].Thickness, 152 md.results.TransientSolution[19].MaskIceLevelset, 153 md.results.TransientSolution[19].CalvingCalvingrate, 154 md.results.TransientSolution[19].CalvingMeltingrate, 155 md.results.TransientSolution[19].BasalforcingsFloatingiceMeltingRate, 156 md.results.TransientSolution[39].Vx, 157 md.results.TransientSolution[39].Vy, 158 md.results.TransientSolution[39].Vel, 159 md.results.TransientSolution[39].Thickness, 160 md.results.TransientSolution[39].MaskIceLevelset, 161 md.results.TransientSolution[39].CalvingCalvingrate, 162 md.results.TransientSolution[39].CalvingMeltingrate, 163 md.results.TransientSolution[39].BasalforcingsFloatingiceMeltingRate]
Note:
See TracChangeset
for help on using the changeset viewer.
