Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 26835) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 26836) @@ -219,4 +219,7 @@ iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum); break; + case AutoregressionLinearFloatingMeltRateEnum: + iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.basin_id",BasalforcingsLinearBasinIdEnum); + break; default: _error_("Basal forcing model "<GetInputValue(&basinid,SmbBasinsIdEnum); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); + if(enum_type==SMBautoregressionEnum) this->GetInputValue(&basinid,SmbBasinsIdEnum); + if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); + if(enum_type==BasalforcingsDeepwaterMeltingRateAutoregressionEnum) this->GetInputValue(&basinid,BasalforcingsLinearBasinIdEnum); for(int i=0;i(oldvarspin); } - if(enum_type==SMBautoregressionEnum) this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); - if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); + if(enum_type==SMBautoregressionEnum) this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); + if(enum_type==FrontalForcingsRignotAutoregressionEnum) this->inputs->SetArrayInput(ThermalforcingValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); + if(enum_type==BasalforcingsDeepwaterMeltingRateAutoregressionEnum) this->inputs->SetArrayInput(BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum,this->lid,varspin,numvertices*arorder); /*Cleanup and return*/ @@ -123,4 +125,10 @@ outenum_type = FrontalForcingsThermalForcingEnum; break; + case(BasalforcingsDeepwaterMeltingRateAutoregressionEnum): + arenum_type = BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum; + basinenum_type = BasalforcingsLinearBasinIdEnum; + noiseenum_type = BasalforcingsDeepwaterMeltingRateNoiseEnum; + outenum_type = BasalforcingsSpatialDeepwaterMeltingRateEnum; + break; } @@ -172,4 +180,35 @@ xDelete(varlist); xDelete(valuesautoregression); +}/*}}}*/ +void Element::BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation){/*{{{*/ + + const int NUM_VERTICES = this->GetNumberOfVertices(); + + int basinid; + IssmDouble alpha; + IssmDouble base[MAXVERTICES]; + IssmDouble values[MAXVERTICES]; + IssmDouble deepwatermelt[MAXVERTICES]; + + /*Get element-specific values*/ + this->GetInputValue(&basinid,BasalforcingsLinearBasinIdEnum); + this->GetInputListOnVertices(&base[0],BaseEnum); + this->GetInputListOnVertices(&deepwatermelt[0],BasalforcingsSpatialDeepwaterMeltingRateEnum); + + /*Compute melt rate at vertices according to basin-specific values of input arguments*/ + for(int i=0;i=upperwaterel[basinid]){ + values[i]=upperwatermelt[basinid]; + } + else if (base[i]AddInput(BasalforcingsFloatingiceMeltingRateEnum,&values[0],P1Enum); }/*}}}*/ void Element::ComputeLambdaS(){/*{{{*/ @@ -2301,5 +2340,4 @@ IssmDouble base[MAXVERTICES]; IssmDouble values[MAXVERTICES]; - IssmDouble perturbation[MAXVERTICES]; IssmDouble time; @@ -2312,5 +2350,4 @@ this->GetInputListOnVertices(&base[0],BaseEnum); - this->GetInputListOnVertices(&perturbation[0],BasalforcingsPerturbationMeltingRateEnum); for(int i=0;i=upperwaterel){ @@ -2324,6 +2361,4 @@ values[i]=deepwatermelt*alpha+(1.-alpha)*upperwatermelt; } - - values[i]+=perturbation[i]; } Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26835) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26836) @@ -70,4 +70,5 @@ void AutoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,int enum_type); void Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,bool isfieldstochastic,int enum_type); + void BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation); void ComputeLambdaS(void); void ComputeNewDamage(); Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 26835) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 26836) @@ -45,4 +45,8 @@ if(VerboseSolution())_printf0_(" call BeckmannGoosse Floating melting rate module\n"); BeckmannGoosseFloatingiceMeltingRatex(femmodel); + break; + case AutoregressionLinearFloatingMeltRateEnum: + if(VerboseSolution())_printf0_(" call Autoregression Linear Floating melting rate module\n"); + AutoregressionLinearFloatingiceMeltingRatex(femmodel); break; default: @@ -206,2 +210,99 @@ } /*}}}*/ +void AutoregressionLinearFloatingiceMeltingRateInitx(FemModel* femmodel){/*{{{*/ + + /*Initialization step of AutoregressionLinearFloatingiceMeltingRatex*/ + int M,N,Nphi,arorder,numbasins,my_rank; + IssmDouble starttime,tstep_ar,tinit_ar; + femmodel->parameters->FindParam(&numbasins,BasalforcingsLinearNumBasinsEnum); + femmodel->parameters->FindParam(&arorder,BasalforcingsAutoregressiveOrderEnum); + IssmDouble* beta0 = NULL; + IssmDouble* beta1 = NULL; + IssmDouble* phi = NULL; + femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + femmodel->parameters->FindParam(&tstep_ar,BasalforcingsAutoregressionTimestepEnum); + femmodel->parameters->FindParam(&tinit_ar,BasalforcingsAutoregressionInitialTimeEnum); + femmodel->parameters->FindParam(&beta0,&M,BasalforcingsBeta0Enum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&beta1,&M,BasalforcingsBeta1Enum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&phi,&M,&Nphi,BasalforcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); + + /*AR model spin-up with 0 noise to initialize BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum (688 = log(0.001)/log(0.99): decaying time of inluence of phi[0]=0.99 to 0.001 of beta_0*/ + int nspin = 688; + for(Object* &object:femmodel->elements->objects){ + Element* element = xDynamicCast(object); //generate element object + element->AutoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,BasalforcingsDeepwaterMeltingRateAutoregressionEnum); + } + /*Cleanup*/ + xDelete(beta0); + xDelete(beta1); + xDelete(phi); +}/*}}}*/ +void AutoregressionLinearFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/ + + /*Get time parameters*/ + IssmDouble time,dt,starttime,tstep_ar; + femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + femmodel->parameters->FindParam(&tstep_ar,BasalforcingsAutoregressionTimestepEnum); + + /*Initialize module at first time step*/ + if(time<=starttime+dt){AutoregressionLinearFloatingiceMeltingRateInitx(femmodel);} + /*Determine if this is a time step for the AR model*/ + bool isstepforar = false; + + #ifndef _HAVE_AD_ + if((fmod(time,tstep_ar)parameters->FindParam(&numbasins,BasalforcingsLinearNumBasinsEnum); + femmodel->parameters->FindParam(&arorder,BasalforcingsAutoregressiveOrderEnum); + IssmDouble tinit_ar; + IssmDouble* beta0 = NULL; + IssmDouble* beta1 = NULL; + IssmDouble* phi = NULL; + IssmDouble* deepwaterel = NULL; + IssmDouble* upperwaterel = NULL; + IssmDouble* upperwatermelt = NULL; + IssmDouble* perturbation = NULL; + + /*Get autoregressive parameters*/ + femmodel->parameters->FindParam(&tinit_ar,BasalforcingsAutoregressionInitialTimeEnum); + femmodel->parameters->FindParam(&beta0,&M,BasalforcingsBeta0Enum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&beta1,&M,BasalforcingsBeta1Enum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&phi,&M,&Nphi,BasalforcingsPhiEnum); _assert_(M==numbasins); _assert_(Nphi==arorder); + + /*Get basin-specific parameters*/ + femmodel->parameters->FindParam(&deepwaterel,&M,BasalforcingsDeepwaterElevationEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&upperwaterel,&M,BasalforcingsUpperwaterElevationEnum); _assert_(M==numbasins); + femmodel->parameters->FindParam(&upperwatermelt,&M,BasalforcingsUpperwaterMeltingRateEnum); _assert_(M==numbasins); + + /*Evaluate whether stochasticity on DeepwaterMeltingRate is requested*/ + femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); + if(isstochastic){ + int numstochasticfields; + int* stochasticfields; + femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum); + femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields); + for(int i=0;i(stochasticfields); + } + /*Time elapsed with respect to AR model initial time*/ + IssmDouble telapsed_ar = time-tinit_ar; + + /*Loop over each element to compute FloatingiceMeltingRate at vertices*/ + for(Object* &object:femmodel->elements->objects){ + Element* element = xDynamicCast(object); + /*Compute autoregression*/ + element->Autoregression(isstepforar,arorder,telapsed_ar,beta0,beta1,phi,isdeepmeltingstochastic,BasalforcingsDeepwaterMeltingRateAutoregressionEnum); + element->BasinLinearFloatingiceMeltingRate(deepwaterel,upperwatermelt,upperwaterel,perturbation); + } +}/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h (revision 26835) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h (revision 26836) @@ -16,4 +16,6 @@ void FloatingiceMeltingRateIsmip6x(FemModel* femmodel); void BeckmannGoosseFloatingiceMeltingRatex(FemModel* femmodel); +void AutoregressionLinearFloatingiceMeltingRateInitx(FemModel* femmodel); +void AutoregressionLinearFloatingiceMeltingRatex(FemModel* femmodel); #endif /* _FloatingiceMeltingRatex_H*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26835) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 26836) @@ -252,4 +252,29 @@ break; case BeckmannGoosseFloatingMeltRateEnum: + break; + case AutoregressionLinearFloatingMeltRateEnum: + /*Add parameters that are not in standard nbvertices format*/ + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsLinearNumBasinsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.ar_order",BasalforcingsAutoregressiveOrderEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.ar_initialtime",BasalforcingsAutoregressionInitialTimeEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.ar_timestep",BasalforcingsAutoregressionTimestepEnum)); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.beta0"); + parameters->AddObject(new DoubleVecParam(BasalforcingsBeta0Enum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.beta1"); + parameters->AddObject(new DoubleVecParam(BasalforcingsBeta1Enum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.phi"); + parameters->AddObject(new DoubleMatParam(BasalforcingsPhiEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.upperwater_melting_rate"); + parameters->AddObject(new DoubleVecParam(BasalforcingsUpperwaterMeltingRateEnum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.upperwater_elevation"); + parameters->AddObject(new DoubleVecParam(BasalforcingsUpperwaterElevationEnum,transparam,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.deepwater_elevation"); + parameters->AddObject(new DoubleVecParam(BasalforcingsDeepwaterElevationEnum,transparam,N)); + xDelete(transparam); break; default: Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26835) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26836) @@ -64,8 +64,14 @@ syn keyword cConstant BalancethicknessStabilizationEnum syn keyword cConstant BarystaticContributionsEnum +syn keyword cConstant BasalforcingsAutoregressionInitialTimeEnum +syn keyword cConstant BasalforcingsAutoregressionTimestepEnum +syn keyword cConstant BasalforcingsAutoregressiveOrderEnum +syn keyword cConstant BasalforcingsBeta0Enum +syn keyword cConstant BasalforcingsBeta1Enum syn keyword cConstant BasalforcingsBottomplumedepthEnum syn keyword cConstant BasalforcingsCrustthicknessEnum syn keyword cConstant BasalforcingsDeepwaterElevationEnum syn keyword cConstant BasalforcingsDeepwaterMeltingRateEnum +syn keyword cConstant BasalforcingsDeepwaterMeltingRateNoiseEnum syn keyword cConstant BasalforcingsDtbgEnum syn keyword cConstant BasalforcingsEnum @@ -77,7 +83,9 @@ syn keyword cConstant BasalforcingsIsmip6NumBasinsEnum syn keyword cConstant BasalforcingsIsmip6TfDepthsEnum +syn keyword cConstant BasalforcingsLinearNumBasinsEnum syn keyword cConstant BasalforcingsLowercrustheatEnum syn keyword cConstant BasalforcingsMantleconductivityEnum syn keyword cConstant BasalforcingsNusseltEnum +syn keyword cConstant BasalforcingsPhiEnum syn keyword cConstant BasalforcingsPicoAverageOverturningEnum syn keyword cConstant BasalforcingsPicoAverageSalinityEnum @@ -594,7 +602,10 @@ syn keyword cConstant BalancethicknessThickeningRateEnum syn keyword cConstant BasalCrevasseEnum +syn keyword cConstant BasalforcingsDeepwaterMeltingRateAutoregressionEnum +syn keyword cConstant BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum syn keyword cConstant BasalforcingsFloatingiceMeltingRateEnum syn keyword cConstant BasalforcingsGeothermalfluxEnum syn keyword cConstant BasalforcingsGroundediceMeltingRateEnum +syn keyword cConstant BasalforcingsLinearBasinIdEnum syn keyword cConstant BasalforcingsPerturbationMeltingRateEnum syn keyword cConstant BasalforcingsSpatialDeepwaterElevationEnum @@ -1212,4 +1223,5 @@ syn keyword cConstant ArrheniusEnum syn keyword cConstant AutodiffJacobianEnum +syn keyword cConstant AutoregressionLinearFloatingMeltRateEnum syn keyword cConstant Balancethickness2AnalysisEnum syn keyword cConstant Balancethickness2SolutionEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26835) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26836) @@ -58,8 +58,14 @@ BalancethicknessStabilizationEnum, BarystaticContributionsEnum, + BasalforcingsAutoregressionInitialTimeEnum, + BasalforcingsAutoregressionTimestepEnum, + BasalforcingsAutoregressiveOrderEnum, + BasalforcingsBeta0Enum, + BasalforcingsBeta1Enum, BasalforcingsBottomplumedepthEnum, BasalforcingsCrustthicknessEnum, BasalforcingsDeepwaterElevationEnum, BasalforcingsDeepwaterMeltingRateEnum, + BasalforcingsDeepwaterMeltingRateNoiseEnum, BasalforcingsDtbgEnum, BasalforcingsEnum, @@ -71,7 +77,9 @@ BasalforcingsIsmip6NumBasinsEnum, BasalforcingsIsmip6TfDepthsEnum, + BasalforcingsLinearNumBasinsEnum, BasalforcingsLowercrustheatEnum, BasalforcingsMantleconductivityEnum, BasalforcingsNusseltEnum, + BasalforcingsPhiEnum, BasalforcingsPicoAverageOverturningEnum, BasalforcingsPicoAverageSalinityEnum, @@ -590,7 +598,10 @@ BalancethicknessThickeningRateEnum, BasalCrevasseEnum, + BasalforcingsDeepwaterMeltingRateAutoregressionEnum, + BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum, BasalforcingsFloatingiceMeltingRateEnum, BasalforcingsGeothermalfluxEnum, BasalforcingsGroundediceMeltingRateEnum, + BasalforcingsLinearBasinIdEnum, BasalforcingsPerturbationMeltingRateEnum, BasalforcingsSpatialDeepwaterElevationEnum, @@ -1211,4 +1222,5 @@ ArrheniusEnum, AutodiffJacobianEnum, + AutoregressionLinearFloatingMeltRateEnum, Balancethickness2AnalysisEnum, Balancethickness2SolutionEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26835) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26836) @@ -66,8 +66,14 @@ case BalancethicknessStabilizationEnum : return "BalancethicknessStabilization"; case BarystaticContributionsEnum : return "BarystaticContributions"; + case BasalforcingsAutoregressionInitialTimeEnum : return "BasalforcingsAutoregressionInitialTime"; + case BasalforcingsAutoregressionTimestepEnum : return "BasalforcingsAutoregressionTimestep"; + case BasalforcingsAutoregressiveOrderEnum : return "BasalforcingsAutoregressiveOrder"; + case BasalforcingsBeta0Enum : return "BasalforcingsBeta0"; + case BasalforcingsBeta1Enum : return "BasalforcingsBeta1"; case BasalforcingsBottomplumedepthEnum : return "BasalforcingsBottomplumedepth"; case BasalforcingsCrustthicknessEnum : return "BasalforcingsCrustthickness"; case BasalforcingsDeepwaterElevationEnum : return "BasalforcingsDeepwaterElevation"; case BasalforcingsDeepwaterMeltingRateEnum : return "BasalforcingsDeepwaterMeltingRate"; + case BasalforcingsDeepwaterMeltingRateNoiseEnum : return "BasalforcingsDeepwaterMeltingRateNoise"; case BasalforcingsDtbgEnum : return "BasalforcingsDtbg"; case BasalforcingsEnum : return "Basalforcings"; @@ -79,7 +85,9 @@ case BasalforcingsIsmip6NumBasinsEnum : return "BasalforcingsIsmip6NumBasins"; case BasalforcingsIsmip6TfDepthsEnum : return "BasalforcingsIsmip6TfDepths"; + case BasalforcingsLinearNumBasinsEnum : return "BasalforcingsLinearNumBasins"; case BasalforcingsLowercrustheatEnum : return "BasalforcingsLowercrustheat"; case BasalforcingsMantleconductivityEnum : return "BasalforcingsMantleconductivity"; case BasalforcingsNusseltEnum : return "BasalforcingsNusselt"; + case BasalforcingsPhiEnum : return "BasalforcingsPhi"; case BasalforcingsPicoAverageOverturningEnum : return "BasalforcingsPicoAverageOverturning"; case BasalforcingsPicoAverageSalinityEnum : return "BasalforcingsPicoAverageSalinity"; @@ -596,7 +604,10 @@ case BalancethicknessThickeningRateEnum : return "BalancethicknessThickeningRate"; case BasalCrevasseEnum : return "BasalCrevasse"; + case BasalforcingsDeepwaterMeltingRateAutoregressionEnum : return "BasalforcingsDeepwaterMeltingRateAutoregression"; + case BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum : return "BasalforcingsDeepwaterMeltingRateValuesAutoregression"; case BasalforcingsFloatingiceMeltingRateEnum : return "BasalforcingsFloatingiceMeltingRate"; case BasalforcingsGeothermalfluxEnum : return "BasalforcingsGeothermalflux"; case BasalforcingsGroundediceMeltingRateEnum : return "BasalforcingsGroundediceMeltingRate"; + case BasalforcingsLinearBasinIdEnum : return "BasalforcingsLinearBasinId"; case BasalforcingsPerturbationMeltingRateEnum : return "BasalforcingsPerturbationMeltingRate"; case BasalforcingsSpatialDeepwaterElevationEnum : return "BasalforcingsSpatialDeepwaterElevation"; @@ -1214,4 +1225,5 @@ case ArrheniusEnum : return "Arrhenius"; case AutodiffJacobianEnum : return "AutodiffJacobian"; + case AutoregressionLinearFloatingMeltRateEnum : return "AutoregressionLinearFloatingMeltRate"; case Balancethickness2AnalysisEnum : return "Balancethickness2Analysis"; case Balancethickness2SolutionEnum : return "Balancethickness2Solution"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26835) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 26836) @@ -57,8 +57,14 @@ syn keyword juliaConstC BalancethicknessStabilizationEnum syn keyword juliaConstC BarystaticContributionsEnum +syn keyword juliaConstC BasalforcingsAutoregressionInitialTimeEnum +syn keyword juliaConstC BasalforcingsAutoregressionTimestepEnum +syn keyword juliaConstC BasalforcingsAutoregressiveOrderEnum +syn keyword juliaConstC BasalforcingsBeta0Enum +syn keyword juliaConstC BasalforcingsBeta1Enum syn keyword juliaConstC BasalforcingsBottomplumedepthEnum syn keyword juliaConstC BasalforcingsCrustthicknessEnum syn keyword juliaConstC BasalforcingsDeepwaterElevationEnum syn keyword juliaConstC BasalforcingsDeepwaterMeltingRateEnum +syn keyword juliaConstC BasalforcingsDeepwaterMeltingRateNoiseEnum syn keyword juliaConstC BasalforcingsDtbgEnum syn keyword juliaConstC BasalforcingsEnum @@ -70,7 +76,9 @@ syn keyword juliaConstC BasalforcingsIsmip6NumBasinsEnum syn keyword juliaConstC BasalforcingsIsmip6TfDepthsEnum +syn keyword juliaConstC BasalforcingsLinearNumBasinsEnum syn keyword juliaConstC BasalforcingsLowercrustheatEnum syn keyword juliaConstC BasalforcingsMantleconductivityEnum syn keyword juliaConstC BasalforcingsNusseltEnum +syn keyword juliaConstC BasalforcingsPhiEnum syn keyword juliaConstC BasalforcingsPicoAverageOverturningEnum syn keyword juliaConstC BasalforcingsPicoAverageSalinityEnum @@ -587,7 +595,10 @@ syn keyword juliaConstC BalancethicknessThickeningRateEnum syn keyword juliaConstC BasalCrevasseEnum +syn keyword juliaConstC BasalforcingsDeepwaterMeltingRateAutoregressionEnum +syn keyword juliaConstC BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum syn keyword juliaConstC BasalforcingsFloatingiceMeltingRateEnum syn keyword juliaConstC BasalforcingsGeothermalfluxEnum syn keyword juliaConstC BasalforcingsGroundediceMeltingRateEnum +syn keyword juliaConstC BasalforcingsLinearBasinIdEnum syn keyword juliaConstC BasalforcingsPerturbationMeltingRateEnum syn keyword juliaConstC BasalforcingsSpatialDeepwaterElevationEnum @@ -1205,4 +1216,5 @@ syn keyword juliaConstC ArrheniusEnum syn keyword juliaConstC AutodiffJacobianEnum +syn keyword juliaConstC AutoregressionLinearFloatingMeltRateEnum syn keyword juliaConstC Balancethickness2AnalysisEnum syn keyword juliaConstC Balancethickness2SolutionEnum Index: /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 26835) +++ /issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp (revision 26836) @@ -255,4 +255,5 @@ case 7: return BasalforcingsIsmip6Enum; case 8: return BeckmannGoosseFloatingMeltRateEnum; + case 9: return AutoregressionLinearFloatingMeltRateEnum; default: _error_("Marshalled Basal Forcings code \""<',0); + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','basalforcings.deepwater_elevation','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); + md = checkfield(md,'fieldname','basalforcings.upperwater_melting_rate','NaN',1,'Inf',1,'>=',0,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); + md = checkfield(md,'fieldname','basalforcings.upperwater_elevation','NaN',1,'Inf',1,'<=',0,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); + md = checkfield(md,'fieldname','basalforcings.basin_id','Inf',1,'>=',0,'<=',md.basalforcings.num_basins,'size',[md.mesh.numberofelements,1]); + md = checkfield(md,'fieldname','basalforcings.beta0','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); + md = checkfield(md,'fieldname','basalforcings.beta1','NaN',1,'Inf',1,'size',[1,md.basalforcings.num_basins],'numel',md.basalforcings.num_basins); + md = checkfield(md,'fieldname','basalforcings.ar_order','numel',1,'NaN',1,'Inf',1,'>=',0); + md = checkfield(md,'fieldname','basalforcings.ar_initialtime','numel',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','basalforcings.ar_timestep','numel',1,'NaN',1,'Inf',1,'>=',md.timestepping.time_step); %autoregression time step cannot be finer than ISSM timestep + md = checkfield(md,'fieldname','basalforcings.phi','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins,md.basalforcings.ar_order]); + end + if ismember('BalancethicknessAnalysis',analyses), + error('not implemented yet!'); + end + if ismember('ThermalAnalysis',analyses) & ~(strcmp(solution,'TransientSolution') & md.transient.isthermal==0), + error('not implemented yet!'); + end + if numel(md.basalforcings.geothermalflux)>1 + md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'timeseries',1,'>=',0); + end + end % }}} + function disp(self) % {{{ + disp(sprintf(' autoregression linear basal forcings parameters:')); + disp(sprintf(' autoregressive model is applied for deepwater_melting_rate')); + + fielddisplay(self,'num_basins','number of different basins [unitless]'); + fielddisplay(self,'basin_id','basin number assigned to each element [unitless]'); + fielddisplay(self,'beta0','basin-specific intercept values [m/yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))'); + fielddisplay(self,'beta1','basin-specific trend values [m yr^(-2)]'); + fielddisplay(self,'ar_order','order of the autoregressive model [unitless]'); + fielddisplay(self,'ar_initialtime','initial time assumed in the autoregressive model parameterization [yr]'); + fielddisplay(self,'ar_timestep','time resolution of the autoregressive model [yr]'); + fielddisplay(self,'phi','basin-specific vectors of lag coefficients [unitless]'); + fielddisplay(self,'deepwater_elevation','basin-specific elevation of ocean deepwater [m]'); + fielddisplay(self,'upperwater_melting_rate','basin-specific basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]'); + fielddisplay(self,'upperwater_elevation','basin-specific elevation of ocean upperwater [m]'); + fielddisplay(self,'groundedice_melting_rate','node-specific basal melting rate (positive if melting) [m/yr]'); + fielddisplay(self,'geothermalflux','node-specific geothermal heat flux [W/m^2]'); + + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + yts=md.constants.yts; + + WriteData(fid,prefix,'name','md.basalforcings.model','data',9,'format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','name','md.basalforcings.geothermalflux','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'fieldname','num_basins','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','ar_order','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','ar_initialtime','format','Double','scale',yts); + WriteData(fid,prefix,'object',self,'fieldname','ar_timestep','format','Double','scale',yts); + WriteData(fid,prefix,'object',self,'fieldname','basin_id','data',self.basin_id-1,'name','md.basalforcings.basin_id','format','IntMat','mattype',2); %0-indexed + WriteData(fid,prefix,'object',self,'fieldname','beta0','format','DoubleMat','name','md.basalforcings.beta0','scale',1./yts,'yts',yts); + WriteData(fid,prefix,'object',self,'fieldname','beta1','format','DoubleMat','name','md.basalforcings.beta1','scale',1./(yts^2),'yts',yts); + WriteData(fid,prefix,'object',self,'fieldname','phi','format','DoubleMat','name','md.basalforcings.phi','yts',yts); + WriteData(fid,prefix,'object',self,'fieldname','deepwater_elevation','format','DoubleMat','name','md.basalforcings.deepwater_elevation'); + WriteData(fid,prefix,'object',self,'fieldname','upperwater_melting_rate','format','DoubleMat','name','md.basalforcings.upperwater_melting_rate','scale',1./yts); + WriteData(fid,prefix,'object',self,'fieldname','upperwater_elevation','format','DoubleMat','name','md.basalforcings.upperwater_elevation'); + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.py =================================================================== --- /issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.py (revision 26836) +++ /issm/trunk-jpl/src/m/classes/autoregressionlinearbasalforcings.py (revision 26836) @@ -0,0 +1,112 @@ +import numpy as np + +from checkfield import * +from fielddisplay import fielddisplay +from project3d import * +from WriteData import * + +class autoregressionlinearbasalforcings(object): + """autoregressionlinearbasalforcings class definition + + Usage: + autoregressionlinearbasalforcings = autoregressionlinearbasalforcings() + """ + + def __init__(self, *args): # {{{ + self.num_basins = 0 + self.beta0 = np.nan + self.beta1 = np.nan + self.ar_order = 0 + self.ar_initialtime = 0 + self.ar_timestep = 0 + self.phi = np.nan + self.basin_id = np.nan + self.groundedice_melting_rate = np.nan + self.deepwater_elevation = np.nan + self.upperwater_melting_rate = np.nan + self.upperwater_elevation = np.nan + self.geothermalflux = np.nan + + nargs = len(args) + if nargs == 0: + self.setdefaultparameters() + else: + raise Exception('constructor not supported') + # }}} + + def __repr__(self): # {{{ + s = ' surface forcings parameters:\n' + s += ' autoregressive model is applied for deepwater_melting_rate\n' + s += '{}\n'.format(fielddisplay(self, 'num_basins', 'number of different basins [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'basin_id', 'basin number assigned to each element [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'beta0', 'basin-specific intercept values [m ice eq./yr] (if beta_1==0 mean=beta_0/(1-sum(phi)))')) + s += '{}\n'.format(fielddisplay(self, 'beta1', 'basin-specific trend values [m ice eq. yr^(-2)]')) + s += '{}\n'.format(fielddisplay(self, 'ar_order', 'order of the autoregressive model [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'ar_initialtime', 'initial time assumed in the autoregressive model parameterization [yr]')) + s += '{}\n'.format(fielddisplay(self, 'ar_timestep', 'time resolution of the autoregressive model [yr]')) + s += '{}\n'.format(fielddisplay(self, 'phi', 'basin-specific vectors of lag coefficients [unitless]')) + s += '{}\n'.format(fielddisplay(self, 'deepwater_elevation', 'basin-specific elevation of ocean deepwater [m]')) + s += '{}\n'.format(fielddisplay(self, 'upperwater_melting_rate', 'basin-specic basal melting rate (positive if melting applied for floating ice whith base >= upperwater_elevation) [m/yr]')) + s += '{}\n'.format(fielddisplay(self, 'upperwater_elevation', 'basin-specific elevation of ocean upperwater [m]')) + s += '{}\n'.format(fielddisplay(self, 'groundedice_melting_rate','node-specific basal melting rate (positive if melting) [m/yr]')) + s += '{}\n'.format(fielddisplay(self, 'geothermalflux','node-specific geothermal heat flux [W/m^2]')) + return s + # }}} + + def setdefaultparameters(self): #{{{ + #Nothing for now + return self + # }}} + + def extrude(self, md): # {{{ + return self # Nothing for now + # }}} + + def initialize(self, md): # {{{ + if np.all(np.isnan(self.groundedice_melting_rate)): + self.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices)) + print(" no basalforcings.groundedice_melting_rate specified: values set as zero") + return self + # }}} + + def checkconsistency(self, md, solution, analyses): # {{{ + if 'MasstransportAnalysis' in analyses: + md = checkfield(md, 'fieldname', 'basalforcings.num_basins', 'numel', 1, 'NaN', 1, 'Inf', 1, '>', 0) + md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1) + md = checkfield(md, 'fieldname', 'basalforcings.deepwater_elevation', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) + md = checkfield(md, 'fieldname', 'basalforcings.upperwater_elevation', 'NaN', 1, 'Inf', 1,'<=',0, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) + md = checkfield(md, 'fieldname', 'basalforcings.upperwater_melting_rate', 'NaN', 1, 'Inf', 1,'>=',0, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) + md = checkfield(md, 'fieldname', 'basalforcings.basin_id', 'Inf', 1, '>=', 0, '<=', md.basalforcings.num_basins, 'size', [md.mesh.numberofelements]) + md = checkfield(md, 'fieldname', 'basalforcings.beta0', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) # Scheme fails if passed as column vector + md = checkfield(md, 'fieldname', 'basalforcings.beta1', 'NaN', 1, 'Inf', 1, 'size', [1, md.basalforcings.num_basins], 'numel', md.basalforcings.num_basins) # Scheme fails if passed as column vector; NOTE: As opposed to MATLAB implementation, pass list + md = checkfield(md, 'fieldname', 'basalforcings.ar_order', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', 0) + md = checkfield(md, 'fieldname', 'basalforcings.ar_initialtime', 'numel', 1, 'NaN', 1, 'Inf', 1) + md = checkfield(md, 'fieldname', 'basalforcings.ar_timestep', 'numel', 1, 'NaN', 1, 'Inf', 1, '>=', md.timestepping.time_step) # Autoregression time step cannot be finer than ISSM timestep + md = checkfield(md, 'fieldname', 'basalforcings.phi', 'NaN', 1, 'Inf', 1, 'size', [md.basalforcings.num_basins, md.basalforcings.ar_order]) + if 'BalancethicknessAnalysis' in analyses: + raise Exception('not implemented yet!') + if 'ThermalAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isthermal: + raise Exception('not implemented yet!') + + return md + # }}} + + def marshall(self, prefix, md, fid): # {{{ + yts = md.constants.yts + + WriteData(fid, prefix, 'name', 'md.basalforcings.model', 'data', 9, 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'groundedice_melting_rate', 'name', 'md.basalforcings.groundedice_melting_rate', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'geothermalflux', 'name', 'md.basalforcings.geothermalflux', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'num_basins', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_order', 'format', 'Integer') + WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_initialtime', 'format', 'Double', 'scale', yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'ar_timestep', 'format', 'Double', 'scale', yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'basin_id', 'data', self.basin_id, 'name', 'md.basalforcings.basin_id', 'format', 'IntMat', 'mattype', 2) # 0-indexed + WriteData(fid, prefix, 'object', self, 'fieldname', 'beta0', 'format', 'DoubleMat', 'name', 'md.basalforcings.beta0', 'scale', 1 / yts, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'beta1', 'format', 'DoubleMat', 'name', 'md.basalforcings.beta1', 'scale', 1 / (yts ** 2), 'yts', yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'phi', 'format', 'DoubleMat', 'name', 'md.basalforcings.phi', 'yts', yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'deepwater_elevation', 'format', 'DoubleMat', 'name', 'md.basalforcings.deepwater_elevation') + WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_melting_rate', 'format', 'DoubleMat', 'name', 'md.basalforcings.upperwater_melting_rate', 'scale', 1 / yts, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'fieldname', 'upperwater_elevation', 'format', 'DoubleMat', 'name', 'md.basalforcings.upperwater_elevation') + + # }}}