Changeset 22840 for issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

Timestamp:

06/11/18 13:42:55 (7 years ago)

Author:

schlegel

Message:

CHG: add RACMO versions to GEMB and allow dakota to run without waitonlock

File:

• issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

@@ -24 +24 @@
 const double LS = 2.8295E6;             // latent heat of sublimation (J kg-1)
 const double SB = 5.67E-8;                // Stefan-Boltzmann constant (W m-2 K-4)
-const double CA = 1005.0;                    // heat capacity of air (J kg-1 k-1)
-const double R = 8.314;                      // gas constant (mol-1 K-1)
+const double CA = 1005.0;                    // heat capacity of air (J kg-1 K-1)
+const double R = 8.314;                      // gas constant (J mol-1 K-1)
 
 void Gembx(FemModel* femmodel){  /*{{{*/
@@ -1854 +1854 @@
         IssmDouble c1=0.0;
         IssmDouble H=0.0;
+        IssmDouble M0=0.0;
+        IssmDouble M1=0.0;
+        IssmDouble c0arth=0.0;
+        IssmDouble c1arth=0.0;
 
         /*output: */
@@ -1885 +1889 @@
                 switch (denIdx){
                         case 1: // Herron and Langway (1980)
-                                c0 = (11.0 * exp(-10160.0 / (T[i] * 8.314))) * C/1000.0;
-                                c1 = (575.0 * exp(-21400.0 / (T[i]* 8.314))) * pow(C/1000.0,.5);
+                                c0 = (11.0 * exp(-10160.0 / (T[i] * R))) * C/1000.0;
+                                c1 = (575.0 * exp(-21400.0 / (T[i]* R))) * pow(C/1000.0,.5);
                                 break;
+
                         case 2: // Arthern et al. (2010) [semi-emperical]
                                 // common variable
                                 // NOTE: Ec=60000, Eg=42400 (i.e. should be in J not kJ)
-                                H = exp((-60000.0/(T[i] * 8.314)) + (42400.0/(Tmean * 8.314))) * (C * 9.81);
+                                H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
                                 c0 = 0.07 * H;
                                 c1 = 0.03 * H;
@@ -1899 +1904 @@
 
                                 // common variable
-                                H = exp((-60.0/(T[i] * 8.314))) * obp[i] / pow(re[i]/1000.0,2.0);
+                                H = exp((-60000.0/(T[i] * R))) * obp[i] / pow(re[i]/1000.0,2.0);
                                 c0 = 9.2e-9 * H;
                                 c1 = 3.7e-9 * H;
@@ -1913 +1918 @@
                                 c0 = (C/dIce) * (76.138 - 0.28965*Tmean)*8.36*pow(CtoK - T[i],-2.061);
                                 c1 = c0;
+                                break;
+
+                        case 6: // Ligtenberg and others (2011) [semi-emperical], Antarctica 
+                                // common variable
+                                H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
+                                c0arth = 0.07 * H;
+                                c1arth = 0.03 * H;
+                                M0 = fmax(1.435 - (0.151 * log(C)),0.25);
+                                M1 = fmax(2.366 - (0.293 * log(C)),0.25);
+                                c0 = M0*c0arth;
+                                c1 = M1*c1arth;
+                                break;
+
+                        case 7: // Kuipers Munneke and others (2015) [semi-emperical], Greenland
+                                // common variable
+                                H = exp((-60000.0/(T[i] * R)) + (42400.0/(T[i] * R))) * (C * 9.81);
+                                c0arth = 0.07 * H;
+                                c1arth = 0.03 * H;
+                                M0 = fmax(1.042 - (0.0916 * log(C)),0.25);
+                                M1 = fmax(1.734 - (0.2039 * log(C)),0.25);
+                                c0 = M0*c0arth;
+                                c1 = M1*c1arth;
                                 break;
                 }