Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23155) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23156) @@ -450,7 +450,7 @@ FrictionPressureAdjustedTemperatureEnum, FrictionQEnum, FrictionWaterLayerEnum, - FrictionWatercolumnMaxEnum, + HydrologyWatercolumnMaxEnum, FrictionTillFrictionAngleEnum, FrictionSedimentCompressibilityCoefficientEnum, GeometryHydrostaticRatioEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23155) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23156) @@ -456,7 +456,7 @@ case FrictionPressureAdjustedTemperatureEnum : return "FrictionPressureAdjustedTemperature"; case FrictionQEnum : return "FrictionQ"; case FrictionWaterLayerEnum : return "FrictionWaterLayer"; - case FrictionWatercolumnMaxEnum : return "FrictionWatercolumnMax"; + case HydrologyWatercolumnMaxEnum : return "HydrologyWatercolumnMax"; case FrictionTillFrictionAngleEnum : return "FrictionTillFrictionAngle"; case FrictionSedimentCompressibilityCoefficientEnum : return "FrictionSedimentCompressibilityCoefficient"; case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23155) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23156) @@ -465,7 +465,7 @@ else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum; else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum; else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum; - else if (strcmp(name,"FrictionWatercolumnMax")==0) return FrictionWatercolumnMaxEnum; + else if (strcmp(name,"HydrologyWatercolumnMax")==0) return HydrologyWatercolumnMaxEnum; else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum; else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum; else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum; Index: ../trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23155) +++ ../trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23156) @@ -709,21 +709,11 @@ element->parameters->FindParam(&e0,FrictionVoidRatioEnum); element->GetInputValue(&Cc,gauss,FrictionSedimentCompressibilityCoefficientEnum); element->GetInputValue(&W,gauss,WatercolumnEnum); - element->GetInputValue(&Wmax,gauss,FrictionWatercolumnMaxEnum); + element->GetInputValue(&Wmax,gauss,HydrologyWatercolumnMaxEnum); -// /*Check that water column height is within 0 and upper bound, correct if needed*/ -// // if watercolumn height is higher than the maximum allowed height, set height to upper bound -// if(W>Wmax){ -// W=Wmax; -// } -// // if watercolumn height is negative (shouldn't happen), set it to 0 -// if(W<0){ -// W=0; -// } -// // if watercolumn height is within 0 and upper bound, nothing to be done -// else{ -// //do nothing -// } + /*Check that water column height is within 0 and upper bound, correct if needed*/ + if(W>Wmax) W=Wmax; + if(W<0) W=0.; N = delta*P0*pow(10.,(e0/Cc)*(1.-W/Wmax)); @@ -745,7 +735,7 @@ IssmDouble u0,q; element->parameters->FindParam(&u0,FrictionThresholdSpeedEnum); element->parameters->FindParam(&q,FrictionPseudoplasticityExponentEnum); - IssmDouble alpha2 = tau_c/(pow(ub,1.-q)*pow(u0,q)); + IssmDouble alpha2 = tau_c/(pow(ub+1.e-10,1.-q)*pow(u0,q)); /*Final checks in debuging mode*/ _assert_(!xIsNan(alpha2)); Index: ../trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp (revision 23155) +++ ../trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp (revision 23156) @@ -35,6 +35,7 @@ iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum); iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum); iomodel->FetchDataToInput(elements,"md.hydrology.drainage_rate",HydrologyDrainageRateEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.watercolumn_max",HydrologyWatercolumnMaxEnum); iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",WatercolumnEnum,0.); }/*}}}*/ void HydrologyPismAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ @@ -114,28 +115,18 @@ IssmDouble* watercolumn = xNew(numvertices); IssmDouble* drainagerate = xNew(numvertices); IssmDouble* meltingrate = xNew(numvertices); -// IssmDouble* watercolumn_max = xNew(numvertices); + IssmDouble* watercolumn_max = xNew(numvertices); element->GetInputListOnVertices(&watercolumn[0],WatercolumnEnum); element->GetInputListOnVertices(&drainagerate[0],HydrologyDrainageRateEnum); element->GetInputListOnVertices(&meltingrate[0],BasalforcingsGroundediceMeltingRateEnum); -// element->GetInputListOnVertices(&watercolumn_max[0],FrictionWatercolumnMaxEnum); + element->GetInputListOnVertices(&watercolumn_max[0],HydrologyWatercolumnMaxEnum); /*Add water*/ -// /*Check that water column height is within 0 and upper bound, correct if needed*/ for(int i=0;iwatercolumn_max[i]){ -// watercolumn[i]=watercolumn_max[i]; -// } -// // if watercolumn height is negative (shouldn't happen), set it to 0 -// if(watercolumn[i]<0){ -// watercolumn[i]=0; -// } -// // if watercolumn height is within 0 and upper bound, nothing to be done -// else{ -// //do nothing -// } + /*Check that water column height is within 0 and upper bound, correct if needed*/ + if(watercolumn[i]>watercolumn_max[i]) watercolumn[i]=watercolumn_max[i]; + if(watercolumn[i]<0) watercolumn[i]=0.; } /* Divide by connectivity, add degree of channelization as an input */ @@ -143,5 +134,7 @@ /*Clean up and return*/ xDelete(watercolumn); + xDelete(meltingrate); + xDelete(watercolumn_max); xDelete(drainagerate); }/*}}}*/ Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23155) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23156) @@ -847,7 +847,6 @@ case 10: iomodel->FetchDataToInput(elements,"md.friction.till_friction_angle",FrictionTillFrictionAngleEnum); iomodel->FetchDataToInput(elements,"md.friction.sediment_compressibility_coefficient",FrictionSedimentCompressibilityCoefficientEnum); - iomodel->FetchDataToInput(elements,"md.friction.watercolumn_max",FrictionWatercolumnMaxEnum); break; default: _error_("friction law "<< frictionlaw <<" not supported"); Index: ../trunk-jpl/src/m/classes/hydrologypism.m =================================================================== --- ../trunk-jpl/src/m/classes/hydrologypism.m (revision 23155) +++ ../trunk-jpl/src/m/classes/hydrologypism.m (revision 23156) @@ -5,7 +5,8 @@ classdef hydrologypism properties (SetAccess=public) - drainage_rate = NaN; + drainage_rate = NaN; + watercolumn_max = NaN; end methods function self = extrude(self,md) % {{{ @@ -35,10 +36,12 @@ end md = checkfield(md,'fieldname','hydrology.drainage_rate','Inf',1,'NaN',1,'>=',0,'size',[md.mesh.numberofvertices 1]); + md = checkfield(md,'fieldname','friction.watercolumn_max','NaN',1,'Inf',1,'>',0.,'size',[md.mesh.numberofvertices 1]); end % }}} function disp(self) % {{{ disp(sprintf(' hydrologypism solution parameters:')); fielddisplay(self,'drainage_rate','fixed drainage rate [mm/yr]'); + fielddisplay(self,'watercolumn_max','maximum water column height [m], recommended default: 2 m'); end % }}} function marshall(self,prefix,md,fid) % {{{ @@ -46,6 +49,7 @@ WriteData(fid,prefix,'name','md.hydrology.model','data',4,'format','Integer'); WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','drainage_rate','format','DoubleMat','mattype',1,'scale',1./(1000.*yts)); %from mm/yr to m/s + WriteData(fid,prefix,'class','hydrology','object',self,'fieldname','watercolumn_max','format','DoubleMat','mattype',1); WriteData(fid,prefix,'data',{'Watercolumn'},'name','md.hydrology.requested_outputs','format','StringArray'); end % }}} end Index: ../trunk-jpl/src/m/classes/frictionpism.m =================================================================== --- ../trunk-jpl/src/m/classes/frictionpism.m (revision 23155) +++ ../trunk-jpl/src/m/classes/frictionpism.m (revision 23156) @@ -11,7 +11,6 @@ void_ratio = 0.; till_friction_angle = NaN; sediment_compressibility_coefficient = NaN; - watercolumn_max = NaN; end methods function self = extrude(self,md) % {{{ @@ -48,7 +47,6 @@ md = checkfield(md,'fieldname','friction.void_ratio','numel',[1],'>',0,'<',1,'NaN',1,'Inf',1); md = checkfield(md,'fieldname','friction.till_friction_angle','NaN',1,'Inf',1,'<',360.,'>',0.,'size',[md.mesh.numberofvertices 1]); %User should give angle in degrees, Matlab calculates in rad md = checkfield(md,'fieldname','friction.sediment_compressibility_coefficient','NaN',1,'Inf',1,'<',1.,'>',0.,'size',[md.mesh.numberofvertices 1]); - md = checkfield(md,'fieldname','friction.watercolumn_max','NaN',1,'Inf',1,'>',0.,'size',[md.mesh.numberofvertices 1]); end % }}} function disp(self) % {{{ disp(sprintf('Basal shear stress parameters for the PISM friction law (See Aschwanden et al. 2016 for more details)')); @@ -58,7 +56,6 @@ fielddisplay(self,'void_ratio','void ratio at a reference effective pressure [dimensionless]'); fielddisplay(self,'till_friction_angle','till friction angle [deg], recommended default: 30 deg'); fielddisplay(self,'sediment_compressibility_coefficient','coefficient of compressibility of the sediment [dimensionless], recommended default: 0.12'); - fielddisplay(self,'watercolumn_max','maximum water column height [m], recommended default: 2 m'); end % }}} function marshall(self,prefix,md,fid) % {{{ yts=md.constants.yts; @@ -70,7 +67,6 @@ WriteData(fid,prefix,'class','friction','object',self,'fieldname','void_ratio','format','Double'); WriteData(fid,prefix,'class','friction','object',self,'fieldname','till_friction_angle','format','DoubleMat','mattype',1); WriteData(fid,prefix,'class','friction','object',self,'fieldname','sediment_compressibility_coefficient','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'class','friction','object',self,'fieldname','watercolumn_max','format','DoubleMat','mattype',1); end % }}} function savemodeljs(self,fid,modelname) % {{{ error('not implemented yet!'); Index: ../trunk-jpl/test/NightlyRun/test612.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test612.m (revision 23155) +++ ../trunk-jpl/test/NightlyRun/test612.m (revision 23156) @@ -6,7 +6,8 @@ %Hydrology md.hydrology = hydrologypism(); -md.hydrology.drainage_rate = 0.001*ones(md.mesh.numberofvertices,1); +md.hydrology.drainage_rate = 0.001*ones(md.mesh.numberofvertices,1); +md.hydrology.watercolumn_max=2*ones(md.mesh.numberofvertices,1); md.initialization.pressure=zeros(md.mesh.numberofvertices,1); md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1); md.basalforcings.groundedice_melting_rate = [1:md.mesh.numberofvertices]'; @@ -15,7 +16,6 @@ md.friction=frictionpism(); md.friction.till_friction_angle=30*ones(md.mesh.numberofvertices,1); md.friction.sediment_compressibility_coefficient=0.12*ones(md.mesh.numberofvertices,1); -md.friction.watercolumn_max=2*ones(md.mesh.numberofvertices,1); md.transient.ishydrology = 1; md.transient.issmb = 0;