source: issm/oecreview/Archive/22819-23185/ISSM-23155-23156.diff

../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -450 +450 @@
         FrictionPressureAdjustedTemperatureEnum,
         FrictionQEnum,
         FrictionWaterLayerEnum,
-        FrictionWatercolumnMaxEnum,
+        HydrologyWatercolumnMaxEnum,
         FrictionTillFrictionAngleEnum,
         FrictionSedimentCompressibilityCoefficientEnum,
         GeometryHydrostaticRatioEnum,

../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -456 +456 @@
                 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";

../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -465 +465 @@
               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;

../trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -709 +709 @@
         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 +735 @@
         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<IssmDouble>(alpha2));

../trunk-jpl/src/c/analyses/HydrologyPismAnalysis.cpp

@@ -35 +35 @@
         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 +115 @@
         IssmDouble* watercolumn  = xNew<IssmDouble>(numvertices);
         IssmDouble* drainagerate = xNew<IssmDouble>(numvertices);
         IssmDouble* meltingrate  = xNew<IssmDouble>(numvertices);
-//      IssmDouble* watercolumn_max  = xNew<IssmDouble>(numvertices);
+        IssmDouble* watercolumn_max  = xNew<IssmDouble>(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;i<numvertices;i++){
                 watercolumn[i] += (meltingrate[i]/rho_ice*rho_water-drainagerate[i])*dt;
-//              // if watercolumn height is higher than the maximum allowed height, set height to upper bound
-//              if(watercolumn[i]>watercolumn_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 +134 @@
 
         /*Clean up and return*/
         xDelete<IssmDouble>(watercolumn);
+        xDelete<IssmDouble>(meltingrate);
+        xDelete<IssmDouble>(watercolumn_max);
         xDelete<IssmDouble>(drainagerate);
 }/*}}}*/

../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

@@ -847 +847 @@
                 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");

../trunk-jpl/src/m/classes/hydrologypism.m

@@ -5 +5 @@
 
 classdef hydrologypism
         properties (SetAccess=public) 
-                drainage_rate = NaN;
+                drainage_rate   = NaN;
+                watercolumn_max = NaN;
         end
         methods
                 function self = extrude(self,md) % {{{
@@ -35 +36 @@
                         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 +49 @@
 
                         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

../trunk-jpl/src/m/classes/frictionpism.m

@@ -11 +11 @@
                 void_ratio                           = 0.;
                 till_friction_angle                  = NaN;
                 sediment_compressibility_coefficient = NaN;
-                watercolumn_max                      = NaN;
         end
         methods
                 function self = extrude(self,md) % {{{
@@ -48 +47 @@
                         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 +56 @@
                         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 +67 @@
                         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!');

../trunk-jpl/test/NightlyRun/test612.m

@@ -6 +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);
 md.basalforcings.groundedice_melting_rate = [1:md.mesh.numberofvertices]';
@@ -15 +16 @@
 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;