Index: /issm/trunk-jpl/m4/analyses.m4 =================================================================== --- /issm/trunk-jpl/m4/analyses.m4 (revision 27461) +++ /issm/trunk-jpl/m4/analyses.m4 (revision 27462) @@ -374,4 +374,18 @@ AC_MSG_RESULT($HAVE_HYDROLOGYSHREVE) dnl }}} +dnl with-HydrologyArmapw{{{ +AC_ARG_WITH([HydrologyArmapw], + AS_HELP_STRING([--with-HydrologyArmapw = YES], [compile with HydrologyArmapw capabilities (default is yes)]), + [HYDROLOGYARMAPW=$withval],[HYDROLOGYARMAPW=yes]) +AC_MSG_CHECKING(for HydrologyArmapw capability compilation) + +HAVE_HYDROLOGYARMAPW=no +if test "x$HYDROLOGYARMAPW" = "xyes"; then + HAVE_HYDROLOGYARMAPW=yes + AC_DEFINE([_HAVE_HYDROLOGYARMAPW_],[1],[with HydrologyArmapw capability]) +fi +AM_CONDITIONAL([HYDROLOGYARMAPW], [test x$HAVE_HYDROLOGYARMAPW = xyes]) +AC_MSG_RESULT($HAVE_HYDROLOGYARMAPW) +dnl }}} dnl with-L2ProjectionBase{{{ AC_ARG_WITH([L2ProjectionBase], Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 27461) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 27462) @@ -487,4 +487,7 @@ issm_sources += ./analyses/HydrologyDCEfficientAnalysis.cpp endif +if HYDROLOGYARMAPW +issm_sources += ./analyses/HydrologyArmapwAnalysis.cpp +endif if L2PROJECTIONEPL issm_sources += ./analyses/L2ProjectionEPLAnalysis.cpp Index: /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 27461) +++ /issm/trunk-jpl/src/c/analyses/EnumToAnalysis.cpp (revision 27462) @@ -74,4 +74,7 @@ case GLheightadvectionAnalysisEnum : return new GLheightadvectionAnalysis(); #endif + #ifdef _HAVE_HYDROLOGYARMAPW_ + case HydrologyArmapwAnalysisEnum : return new HydrologyArmapwAnalysis(); + #endif #ifdef _HAVE_HYDROLOGYDCEFFICIENT_ case HydrologyDCEfficientAnalysisEnum : return new HydrologyDCEfficientAnalysis(); Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 27461) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 27462) @@ -28,4 +28,5 @@ #include "./LoveAnalysis.h" #include "./EsaAnalysis.h" +#include "./HydrologyArmapwAnalysis.h" #include "./HydrologyDCEfficientAnalysis.h" #include "./HydrologyDCInefficientAnalysis.h" Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27461) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27462) @@ -105,4 +105,11 @@ noiseenum_type = SubglacialdischargeARMANoiseEnum; outenum_type = FrontalForcingsSubglacialDischargeEnum; + break; + case(HydrologyarmapwEnum): + arenum_type = WaterPressureValuesAutoregressionEnum; + maenum_type = WaterPressureValuesMovingaverageEnum; + basinenum_type = HydrologyBasinsIdEnum; + noiseenum_type = FrictionWaterPressureNoiseEnum; + outenum_type = WaterPressureArmaPerturbationEnum; break; } @@ -2733,4 +2740,5 @@ /*Variable declaration*/ + bool ratevariable; const int numvertices = this->GetNumberOfVertices(); int basinid,mindex,mindexnext,basinenum_type,varenum_type,indperiod; @@ -2749,5 +2757,11 @@ basinenum_type = FrontalForcingsBasinIdEnum; varenum_type = FrontalForcingsSubglacialDischargeEnum; - break; + ratevariable = true; + break; + case(HydrologyarmapwEnum): + basinenum_type = HydrologyBasinsIdEnum; + varenum_type = FrictionWaterPressureEnum; + ratevariable = false; + break; } @@ -2787,10 +2801,13 @@ for(int v=0;vAddInput(varenum_type,varlist,P1Enum); + this->AddInput(varenum_type,varlist,P1DGEnum); /*Clean-up*/ @@ -4713,15 +4730,28 @@ void Element::SubglacialWaterPressure(int output_enum){/*{{{*/ - int numvertices = this->GetNumberOfVertices(); - IssmDouble p_water[numvertices]; - Gauss* gauss=this->NewGauss(); + bool ispwHydroArma; + int M; + int numvertices = this->GetNumberOfVertices(); + IssmDouble p_water[numvertices]; + IssmDouble* perturbationvalues = xNew(numvertices); + Gauss* gauss=this->NewGauss(); Friction* friction = new Friction(this); + /*Calculate subglacial water pressure*/ for(int i=0;iGaussVertex(i); p_water[i] = friction->SubglacialWaterPressure(gauss); - } - this->AddInput(output_enum,p_water,P1DGEnum); - delete gauss; + } + /*Add perturbation in water pressure if HydrologyIsWaterPressureArmaEnum is true*/ + this->parameters->FindParam(&ispwHydroArma,HydrologyIsWaterPressureArmaEnum); + if(ispwHydroArma){ + this->GetInputListOnVertices(perturbationvalues,WaterPressureArmaPerturbationEnum); + for(int i=0;iAddInput(output_enum,p_water,P1DGEnum); + /*Clean-up*/ + delete gauss; delete friction; + xDelete(perturbationvalues); }/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 27461) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 27462) @@ -796,4 +796,7 @@ analyses_temp[numanalyses++]=L2ProjectionEPLAnalysisEnum; analyses_temp[numanalyses++]=L2ProjectionBaseAnalysisEnum; + } + if(hydrology_model==HydrologyarmapwEnum){ + analyses_temp[numanalyses++]=HydrologyArmapwAnalysisEnum; } } Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 27461) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 27462) @@ -344,19 +344,20 @@ /*Get effective pressure*/ - bool ispwStochastic; - IssmDouble Neff; - element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum); - if(ispwStochastic){ - /*Retrieve stochastic water pressure and compute ice pressure*/ - IssmDouble p_ice,p_water,Neff_limit; - element->GetInputValue(&p_water,gauss,FrictionWaterPressureEnum); - element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum); - p_ice = IcePressure(gauss); - Neff = max(Neff_limit*p_ice, p_ice - p_water); - } - else{ - /*Compute effective pressure directly*/ - Neff = EffectivePressure(gauss); - } + bool ispwHydro,ispwStochastic; + IssmDouble Neff; + element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum); + element->parameters->FindParam(&ispwHydro,HydrologyIsWaterPressureArmaEnum); + if(ispwStochastic || ispwHydro){ + /*Retrieve pre-computed water pressure and compute ice pressure*/ + IssmDouble p_ice,p_water,Neff_limit; + element->GetInputValue(&p_water,gauss,FrictionWaterPressureEnum); + element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum); + p_ice = IcePressure(gauss); + Neff = max(Neff_limit*p_ice, p_ice - p_water); + } + else{ + /*Compute effective pressure directly*/ + Neff = EffectivePressure(gauss); + } /*Get velocity magnitude*/ @@ -559,15 +560,16 @@ IssmDouble vmag = VelMag(gauss); - bool ispwStochastic; - IssmDouble Neff; - element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum); - if(ispwStochastic){ - /*Retrieve stochastic water pressure and compute ice pressure*/ - IssmDouble p_ice,p_water,Neff_limit; - element->GetInputValue(&p_water,gauss,FrictionWaterPressureEnum); + bool ispwHydro,ispwStochastic; + IssmDouble Neff; + element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum); + element->parameters->FindParam(&ispwHydro,HydrologyIsWaterPressureArmaEnum); + if(ispwStochastic || ispwHydro){ + /*Retrieve pre-computed water pressure and compute ice pressure*/ + IssmDouble p_ice,p_water,Neff_limit; + element->GetInputValue(&p_water,gauss,FrictionWaterPressureEnum); element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum); - p_ice = IcePressure(gauss); - Neff = max(Neff_limit*p_ice, p_ice - p_water); - } + p_ice = IcePressure(gauss); + Neff = max(Neff_limit*p_ice, p_ice - p_water); + } else{ /*Compute effective pressure directly*/ Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 27461) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 27462) @@ -261,4 +261,13 @@ delete analysis; } + + /*Using the armaPw hydrology model*/ + else if (hydrology_model==HydrologyarmapwEnum){ + femmodel->SetCurrentConfiguration(HydrologyArmapwAnalysisEnum); + if(VerboseSolution()) _printf0_(" updating subglacial water pressure\n"); + HydrologyArmapwAnalysis* analysis = new HydrologyArmapwAnalysis(); + analysis->UpdateSubglacialWaterPressure(femmodel); + delete analysis; + } else{ _error_("Hydrology model "<< EnumToStringx(hydrology_model) <<" not supported yet"); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27461) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27462) @@ -496,4 +496,5 @@ int hydrology_model; iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + parameters->AddObject(new BoolParam(HydrologyIsWaterPressureArmaEnum,false)); if(hydrology_model==HydrologydcEnum){ IssmDouble sedcomp, sedporo, watcomp, rhofresh, g; @@ -533,4 +534,28 @@ /*Nothing to add*/ } + else if(hydrology_model==HydrologyarmapwEnum){ + parameters->SetParam(true,HydrologyIsWaterPressureArmaEnum); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.num_basins",HydrologyNumBasinsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.num_breaks",HydrologyarmaNumBreaksEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.num_params",HydrologyarmaNumParamsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.ar_order",HydrologyarmaarOrderEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.ma_order",HydrologyarmamaOrderEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.arma_timestep",HydrologyarmaTimestepEnum)); + iomodel->FetchData(&transparam,&M,&N,"md.hydrology.datebreaks"); + parameters->AddObject(new DoubleMatParam(HydrologyarmadatebreaksEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.hydrology.polynomialparams"); + parameters->AddObject(new DoubleMatParam(HydrologyarmapolyparamsEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.hydrology.arlag_coefs"); + parameters->AddObject(new DoubleMatParam(HydrologyarmaarlagcoefsEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.hydrology.malag_coefs"); + parameters->AddObject(new DoubleMatParam(HydrologyarmamalagcoefsEnum,transparam,M,N)); + xDelete(transparam); + iomodel->FetchData(&transparam,&M,&N,"md.hydrology.monthlyfactors"); + parameters->AddObject(new DoubleMatParam(HydrologyarmaMonthlyFactorsEnum,transparam,M,N)); + xDelete(transparam); + } else{ _error_("Hydrology model "<parameters->FindParam(&starttime,TimesteppingStartTimeEnum); femmodel->parameters->FindParam(&tstep_stoch,StochasticForcingTimestepEnum); - + + /*Check if we use HydroarmaPw*/ + bool ispwHydro; + femmodel->parameters->FindParam(&ispwHydro,HydrologyIsWaterPressureArmaEnum); + #ifndef _HAVE_AD_ if((fmod(time,tstep_stoch)