Index: /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp	(revision 15808)
+++ /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp	(revision 15808)
@@ -0,0 +1,110 @@
+/* \file LliboutryDuval.cpp
+ * \brief figure out B of ice for a certain temperature and water fraction or enthalpy
+ */
+
+#include <math.h>
+#include "../Numerics/types.h"
+#include "../Exceptions/exceptions.h"
+
+
+/* get ice stiffness B from enthalpy, pressure and flow law exponent*/
+IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure,IssmDouble n){
+  /*Use parameterization for the rheology: Aschwanden 2012
+   *
+   *  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)
+   *  
+   *  T(H,p) = Tref + H/c_i, if H < H_s(p)
+   *         = Tpmp , if H_s(p) \le H \le H_l(p)
+   *
+   *  w(H,p) = 0, if H < H_s(p)
+   *         = (H - H_s(p))/L
+   *
+   *  H_s(p) = c_i (Tpmp - Tref)
+   *
+   *  Tpmp   = T - beta p;
+   *
+   *  A0 constant of proportionality
+   *     = 3.61 * 10^-13   if T*<263.15K
+   *     = 1.73 * 10^3     if T*>263.15K
+   *  Q  Activation energy for creep
+   *     = 6.0  * 10^4     if T*<263.15K
+   *     = 13.9 * 10^4     if T*>263.15K
+   *  R  Universal gas constant
+   *     = 8.314
+   *  
+   *  Convert A to B :  B = A^(-1/n) */
+  
+  /*Some physical constants (Aschwanden 2012)*/
+  /*TODO: get those constants from model*/
+  IssmDouble beta=7.9*pow(10.,-8.);
+  IssmDouble R=8.314;
+  IssmDouble heatcapacity=2009; // J/kg/K
+  IssmDouble Tref=253.15;
+  IssmDouble latentheat=3.34*pow(10,5.); // from Aschwanden 2012
+
+  /*Intermediaries*/
+  IssmDouble A,B,Tstar,Tpmp,H_sp,waterfraction;
+	
+  _assert_(pressure);
+  _assert_(enthalpy);
+  Tpmp=273.15-beta*pressure; 
+  H_sp=heatcapacity*(Tpmp - Tref);
+  if (enthalpy < H_sp){
+    Tstar = Tref + enthalpy/heatcapacity - beta*pressure;	
+    waterfraction = 0;
+  }
+  else{
+    Tstar=Tpmp;
+    waterfraction=(enthalpy - H_sp)/latentheat;
+  }
+
+  /*Get A*/
+  if(Tstar<263.15){
+    A=3.61*pow(10.,-13.) * exp(  -6.*pow(10.,4.)/(R*Tstar));
+  }
+  else{
+    A=1.73*pow(10.,  3.) * exp(-13.9*pow(10.,4.)/(R*Tstar));
+  }
+  A*=(1 + 181.25*waterfraction);
+
+  /*Convert to B*/
+  _assert_(n);
+  B=pow(A,-1./n);
+
+  return B;
+}
+
+/*Get stiffness from temperature, waterfraction and depth*/
+IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n){
+	/*Use parameterization for the rheology: Greve and Blatter 2009
+	 * get enthalpy from temperature and water fraction,
+	 * and use LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure,IssmDouble n) */
+  
+  IssmDouble rho_ice=910; // kg/m^3
+  IssmDouble g=9.81; //kg*m/s^2
+  IssmDouble heatcapacity=2009; // J/kg/K
+  IssmDouble Tref=253.15;
+  IssmDouble beta=7.9*pow(10.,-8.);
+  IssmDouble latentheat=3.34*pow(10,5.); // from Aschwanden 2012
+  IssmDouble Tstar, enthalpy, pressure, B;
+  _assert_(temperature);
+  _assert_(waterfraction);
+  _assert_(depth);
+
+  /*get pressure*/
+  pressure= rho_ice*g*depth;
+  Tstar=temperature-beta*pressure; // TODO: check whether plus or minus
+  /*get enthalpy*/
+  if (Tstar < 273.15){
+    enthalpy=heatcapacity*(Tstar - Tref);
+  }
+  else{
+    enthalpy=heatcapacity*(273.15 - Tref) + waterfraction*latentheat;
+  }
+	
+  B=LliboutryDuval(enthalpy, pressure, n);
+
+  return B;
+}
+