source: issm/oecreview/Archive/19101-20495/ISSM-20091-20092.diff

../trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp

@@ -8 +8 @@
 
 /* get ice stiffness B from enthalpy, pressure and flow law exponent*/
 IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat){
-  /*Use parameterization for the rheology: Grewe/Blatter 2009, Aschwanden 2012
+  /*Use Lliboutry & Duval's 1985 parameterization for the rheology: 
+        * see also: Grewe/Blatter 2009, Aschwanden et al. 2012
+        *
+        * ISSM uses enthalpy/temperature values that are not corrected for pressure.
    *
    *  A(H,p) = A0 exp(-Q/RT(H,p)), if H < H_s(p)
    *         = A0 exp(-Q/RTpmp) (1+181.25w(H,p)), if H_s(p) \le H < H_l(p)
@@ -33 +36 @@
    *     = 8.314
    *  
    *  Convert A to B :  B = A^(-1/n) */
-
         /*check feasibility*/
-  _assert_(pressure>=0);
-  _assert_(n>0);
-  _assert_(betaCC>=0);
-  _assert_(referencetemperature>=0);
-  _assert_(heatcapacity>0);
-  _assert_(latentheat>0);
+        _assert_(pressure+1.e-4>=0); // deal with pressure instability at ice surface
+        _assert_(n>0);
+        _assert_(betaCC>=0);
+        _assert_(referencetemperature>=0);
+        _assert_(heatcapacity>0);
+        _assert_(latentheat>0);
 
-  /*Some physical constants*/
-  IssmDouble R=8.314; 
+        /*Some physical constants*/
+        IssmDouble R=8.314; 
 
-  /*Intermediaries*/
-  IssmDouble A,B,Tstar,Tpmp,H_sp,waterfraction;
+        /*Intermediaries*/
+        IssmDouble A,B,Tstar,Tpmp,H_sp,waterfraction;
 
-  Tpmp=273.15-betaCC*pressure; 
-  H_sp=heatcapacity*(Tpmp - referencetemperature);
-  if (enthalpy < H_sp){
-    Tstar = referencetemperature + enthalpy/heatcapacity - betaCC*pressure;     
-    waterfraction = 0.;
-  }
-  else{
-    Tstar=Tpmp;
-    waterfraction=(enthalpy - H_sp)/latentheat;
-                if (waterfraction > 0.01)
-                        waterfraction = 0.01;
-  }
+        Tpmp=273.15-betaCC*pressure; //pressure melting point temperature
+        H_sp=heatcapacity*(Tpmp-referencetemperature); //pressure melting point enthalpy
+        if (enthalpy<H_sp){ //compute homologous temperature and water fraction
+                Tstar=referencetemperature+enthalpy/heatcapacity+betaCC*pressure; 
+                waterfraction=0.;
+        }
+        else{
+                Tstar=273.15;
+                waterfraction=(enthalpy-H_sp)/latentheat;
+                if (waterfraction>0.01) waterfraction=0.01; // limit softness of ice
+        }
 
-  /*Get A*/
-  if(Tstar<=263.15){
-    A=3.61e-13 * exp(  -6.e+4/(R*Tstar));
-  }
-  else{
-    A=1.73e3   * exp(-13.9e+4/(R*Tstar));
-  }
-  A*=(1. + 181.25*waterfraction);
+        /*Get A*/
+        if(Tstar<=263.15){A=3.61e-13*exp(-6.e+4/(R*Tstar));}
+        else{A=1.73e3*exp(-13.9e+4/(R*Tstar));}
+        A*=(1.+181.25*waterfraction);
 
-  /*Convert to B*/
-  B=pow(A,-1./n);
-
-  return B;
+        /*Convert to B*/
+        B=pow(A,-1./n);
+        return B;
 }