Changeset 24140 for issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

Timestamp:

09/11/19 09:30:58 (6 years ago)

Author:

rueckamp

Message:

CHG: add drainiceolumn and watercolumn_upperlimit as user input; improve drainage for steadystate case

File:

• issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (modified) (6 diffs)

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

@@ -275 +275 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isenthalpy",ThermalIsenthalpyEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdynamicbasalspc",ThermalIsdynamicbasalspcEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.thermal.isdrainicecolumn",ThermalIsdrainicecolumnEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.thermal.watercolumn_upperlimit",ThermalWatercolumnUpperlimitEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
 
@@ -424 +426 @@
         element->GetInputListOnNodes(basalmeltingrates,BasalforcingsGroundediceMeltingRateEnum);
 
+        IssmDouble watercolumnupperlimit = element->FindParam(ThermalWatercolumnUpperlimitEnum);
+        
         Gauss* gauss=element->NewGauss();
         for(is=0;is<numsegments;is++){
@@ -484 +488 @@
                                 watercolumns[nodedown]+=dt*meltingrate_enthalpy[is]; 
                         }
+                        if(watercolumns[nodedown]>watercolumnupperlimit) watercolumns[nodedown]=watercolumnupperlimit;
                 }
                 else{
@@ -1168 +1173 @@
 
                 for(k=0;k<numnodes;k++){
-                        waterfractions[k]-=dt*drainage[k];
+                        if(dt==0.)
+                                waterfractions[k]-=drainage[k];
+                        else
+                                waterfractions[k]-=dt*drainage[k];
+                        
                         element->ThermalToEnthalpy(&enthalpies[k], temperatures[k], waterfractions[k], pressures[k]);
                 }
@@ -1187 +1196 @@
         IssmDouble thermalconductivity      = element->FindParam(MaterialsThermalconductivityEnum);
 
-        if(enthalpy < PureIceEnthalpy(element,pressure)){
+        if(enthalpy < PureIceEnthalpy(element,pressure))
                 return thermalconductivity/heatcapacity;
-        }
-        else{
+        else
                 return temperateiceconductivity/heatcapacity;
-        }
 }/*}}}*/
 IssmDouble     EnthalpyAnalysis::EnthalpyDiffusionParameterVolume(Element* element,int enthalpy_enum){/*{{{*/
@@ -1678 +1685 @@
         /*Intermediaries*/
         bool computebasalmeltingrates=true;
-        bool drainicecolumn=true;
+        bool isdrainicecolumn;
         IssmDouble dt;
 
         femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-
-        if(drainicecolumn && (dt>0.))   DrainWaterfraction(femmodel);
+        femmodel->parameters->FindParam(&isdrainicecolumn,ThermalIsdrainicecolumnEnum);
+
+        if(isdrainicecolumn)    DrainWaterfraction(femmodel);
         if(computebasalmeltingrates)    ComputeBasalMeltingrate(femmodel);