Changeset 17944

Timestamp:

05/06/14 06:52:35 (11 years ago)

Author:

jbondzio

Message:

BUG: compute basal melt rate using ice density, and use relationship (1 - w)*rho_mixture = rho_ice

File:

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

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

@@ -911 +911 @@
         IssmDouble  heatflux;
         IssmDouble  vec_heatflux[dim],normal_base[dim],d1enthalpy[dim];
-        IssmDouble  temperature, waterfraction;
         IssmDouble  basalfriction,alpha2;
         IssmDouble  dt,yts;
@@ -926 +925 @@
         IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum);
         IssmDouble rho_ice    = element->GetMaterialParameter(MaterialsRhoIceEnum);
-        IssmDouble rho_water  = element->GetMaterialParameter(MaterialsRhoWaterEnum);
+        IssmDouble rho_water  = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
         Input* enthalpy_input         = element->GetInput(EnthalpyEnum);                    _assert_(enthalpy_input);
         Input* geothermalflux_input   = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
@@ -996 +995 @@
                         basalfriction=alpha2*(vx*vx + vy*vy + vz*vz);
 
-                        element->EnthalpyToThermal(&temperature,&waterfraction,enthalpy[vertexdown],pressure[vertexdown]);
                         geothermalflux_input->GetInputValue(&geothermalflux,gauss);
-                        /* -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66*/
+                        /* -Mb= Fb-(q-q_geo)/((1-w)*L*rho), and (1-w)*rho=rho_ice, cf Aschwanden 2012, eqs.1, 2, 66*/
                         heating[is]=(heatflux+basalfriction+geothermalflux);
-                        meltingrate_enthalpy[is]=heating[is]/((1-waterfraction)*latentheat*rho_water); // m/s water equivalent
+                        meltingrate_enthalpy[is]=heating[is]/(latentheat*rho_ice); // m/s water equivalent
                 }
         }