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

Index: ../trunk-jpl/test/Archives/Archive436.nc
===================================================================
Cannot display: file marked as a binary type.
svn:mime-type = application/octet-stream
Index: ../trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp
===================================================================
--- ../trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp       (revision 20091)
+++ ../trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp       (revision 20092)
@@ -8,7 +8,10 @@
 
 /* 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,45 +36,38 @@
    *     = 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;
 }