Index: /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp (revision 21525) +++ /issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp (revision 21526) @@ -148,4 +148,5 @@ parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum)); }/*}}}*/ @@ -173,4 +174,5 @@ ElementMatrix* Ke = element->NewElementMatrix(); IssmDouble* dbasis = xNew(2*numnodes); + IssmDouble* basis = xNew(numnodes); /*Retrieve all inputs and parameters*/ @@ -179,4 +181,9 @@ /*Get conductivity from inputs*/ IssmDouble conductivity = GetConductivity(element); + + /*Get englacial storage coefficient*/ + IssmDouble storage,dt; + element->FindParam(&storage,HydrologyStorageEnum); + element->FindParam(&dt,TimesteppingTimeStepEnum); /* Start looping on the number of gaussian points: */ @@ -187,8 +194,10 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + element->NodalFunctions(basis,gauss); for(int i=0;ivalues[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]); + Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]) + + gauss->weight*Jdet*storage/dt*basis[i]*basis[j]; } } @@ -244,4 +253,9 @@ IssmDouble conductivity = GetConductivity(element); + /*Get englacial storage coefficient*/ + IssmDouble storage,dt; + element->FindParam(&storage,HydrologyStorageEnum); + element->FindParam(&dt,TimesteppingTimeStepEnum); + /*Build friction element, needed later: */ Friction* friction=new Friction(element,2); @@ -308,4 +322,5 @@ -beta*sqrt(vx*vx+vy*vy) +ieb + +storage*head_old/dt )*basis[i]; } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21525) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21526) @@ -176,4 +176,5 @@ HydrologyRelaxationEnum, HydrologyBasalFluxEnum, + HydrologyStorageEnum, InversionControlParametersEnum, InversionControlScalingFactorsEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21525) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21526) @@ -182,4 +182,5 @@ case HydrologyRelaxationEnum : return "HydrologyRelaxation"; case HydrologyBasalFluxEnum : return "HydrologyBasalFlux"; + case HydrologyStorageEnum : return "HydrologyStorage"; case InversionControlParametersEnum : return "InversionControlParameters"; case InversionControlScalingFactorsEnum : return "InversionControlScalingFactors"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21525) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21526) @@ -185,4 +185,5 @@ else if (strcmp(name,"HydrologyRelaxation")==0) return HydrologyRelaxationEnum; else if (strcmp(name,"HydrologyBasalFlux")==0) return HydrologyBasalFluxEnum; + else if (strcmp(name,"HydrologyStorage")==0) return HydrologyStorageEnum; else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum; else if (strcmp(name,"InversionControlScalingFactors")==0) return InversionControlScalingFactorsEnum; @@ -259,9 +260,9 @@ else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum; else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum; - else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum; + if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum; + else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum; else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum; else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum; - else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; + if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; + else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum; else if (strcmp(name,"Temperature")==0) return TemperatureEnum; - else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum; + if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; + else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum; else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum; else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum; - else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum; else stage=6; } if(stage==6){ - if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum; + if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum; + else if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum; else if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum; else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum; @@ -751,9 +752,9 @@ else if (strcmp(name,"Mpi")==0) return MpiEnum; else if (strcmp(name,"Mumps")==0) return MumpsEnum; - else if (strcmp(name,"Gsl")==0) return GslEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"Cuffey")==0) return CuffeyEnum; + if (strcmp(name,"Gsl")==0) return GslEnum; + else if (strcmp(name,"Cuffey")==0) return CuffeyEnum; else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum; - else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum; + if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum; + else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum; else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum; else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum; @@ -997,9 +998,9 @@ else if (strcmp(name,"Materials")==0) return MaterialsEnum; else if (strcmp(name,"Nodes")==0) return NodesEnum; - else if (strcmp(name,"Contours")==0) return ContoursEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"Parameters")==0) return ParametersEnum; + if (strcmp(name,"Contours")==0) return ContoursEnum; + else if (strcmp(name,"Parameters")==0) return ParametersEnum; else if (strcmp(name,"Vertices")==0) return VerticesEnum; else if (strcmp(name,"Results")==0) return ResultsEnum; Index: /issm/trunk-jpl/src/m/classes/hydrologysommers.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologysommers.m (revision 21525) +++ /issm/trunk-jpl/src/m/classes/hydrologysommers.m (revision 21526) @@ -16,4 +16,5 @@ neumannflux = NaN; relaxation = 0; + storage = 0; end methods @@ -33,4 +34,5 @@ % Set under-relaxation parameter to be 1 (no under-relaxation of nonlinear iteration) self.relaxation=1; + self.storage=0; end % }}} function md = checkconsistency(self,md,solution,analyses) % {{{ @@ -51,4 +53,5 @@ md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','hydrology.relaxation','>=',0); + md = checkfield(md,'fieldname','hydrology.storage','>=',0); end % }}} function disp(self) % {{{ @@ -64,4 +67,5 @@ fielddisplay(self,'spchead','water head constraints (NaN means no constraint) (m)'); fielddisplay(self,'relaxation','under-relaxation coefficient for nonlinear iteration'); + fielddisplay(self,'storage','englacial storage coefficient (void ratio)'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -80,4 +84,5 @@ WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','spchead','format','DoubleMat','mattype',1); WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','relaxation','format','Double'); + WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','storage','format','Double'); end % }}} end