Index: /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp (revision 23155) +++ /issm/trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp (revision 23156) @@ -36,4 +36,5 @@ 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.); }/*}}}*/ @@ -115,26 +116,16 @@ 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.; } @@ -144,4 +135,6 @@ /*Clean up and return*/ xDelete(watercolumn); + xDelete(meltingrate); + xDelete(watercolumn_max); xDelete(drainagerate); }/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23155) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 23156) @@ -848,5 +848,4 @@ 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: Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23155) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 23156) @@ -710,19 +710,9 @@ element->GetInputValue(&Cc,gauss,FrictionSedimentCompressibilityCoefficientEnum); element->GetInputValue(&W,gauss,WatercolumnEnum); - element->GetInputValue(&Wmax,gauss,FrictionWatercolumnMaxEnum); - -// /*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 -// } + element->GetInputValue(&Wmax,gauss,HydrologyWatercolumnMaxEnum); + + /*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)); @@ -746,5 +736,5 @@ 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*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23155) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23156) @@ -451,5 +451,5 @@ FrictionQEnum, FrictionWaterLayerEnum, - FrictionWatercolumnMaxEnum, + HydrologyWatercolumnMaxEnum, FrictionTillFrictionAngleEnum, FrictionSedimentCompressibilityCoefficientEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23155) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23156) @@ -457,5 +457,5 @@ case FrictionQEnum : return "FrictionQ"; case FrictionWaterLayerEnum : return "FrictionWaterLayer"; - case FrictionWatercolumnMaxEnum : return "FrictionWatercolumnMax"; + case HydrologyWatercolumnMaxEnum : return "HydrologyWatercolumnMax"; case FrictionTillFrictionAngleEnum : return "FrictionTillFrictionAngle"; case FrictionSedimentCompressibilityCoefficientEnum : return "FrictionSedimentCompressibilityCoefficient"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23155) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23156) @@ -466,5 +466,5 @@ 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; Index: /issm/trunk-jpl/src/m/classes/frictionpism.m =================================================================== --- /issm/trunk-jpl/src/m/classes/frictionpism.m (revision 23155) +++ /issm/trunk-jpl/src/m/classes/frictionpism.m (revision 23156) @@ -12,5 +12,4 @@ till_friction_angle = NaN; sediment_compressibility_coefficient = NaN; - watercolumn_max = NaN; end methods @@ -49,5 +48,4 @@ 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) % {{{ @@ -59,5 +57,4 @@ 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) % {{{ @@ -71,5 +68,4 @@ 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) % {{{ Index: /issm/trunk-jpl/src/m/classes/hydrologypism.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologypism.m (revision 23155) +++ /issm/trunk-jpl/src/m/classes/hydrologypism.m (revision 23156) @@ -6,5 +6,6 @@ classdef hydrologypism properties (SetAccess=public) - drainage_rate = NaN; + drainage_rate = NaN; + watercolumn_max = NaN; end methods @@ -36,8 +37,10 @@ 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) % {{{ @@ -47,4 +50,5 @@ 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 % }}} Index: /issm/trunk-jpl/test/NightlyRun/test612.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test612.m (revision 23155) +++ /issm/trunk-jpl/test/NightlyRun/test612.m (revision 23156) @@ -7,5 +7,6 @@ %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); @@ -16,5 +17,4 @@ 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;