Changeset 16545

Timestamp:

10/25/13 04:29:40 (11 years ago)

Author:

jbondzio

Message:

CHG: enthalpydiffusion parameter smoothed
BUG: water column stays nonnegative now, minor in drainage function

Location:

issm/trunk-jpl

Files:

• src/c/classes/Elements/Penta.cpp (modified) (10 diffs)
• src/c/classes/Materials/Matpar.cpp (modified) (1 diff)
• test/Archives/Archive120.nc (modified) ( previous)
• test/Archives/Archive121.nc (modified) ( previous)
• test/Archives/Archive122.nc (modified) ( previous)
• test/Archives/Archive325.nc (modified) ( previous)
• test/Archives/Archive327.nc (modified) ( previous)
• test/Archives/Archive431.nc (modified) ( previous)
• test/Archives/Archive432.nc (modified) ( previous)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -3824 +3824 @@
 
                 /*Conduction: */  
-                /*Need to change that depending on enthalpy value -> cold or temperate ice: */  
                 GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
 
@@ -3873 +3872 @@
                 }
 
-                /*Artifficial diffusivity*/
+                /*Artificial diffusivity*/
                 if(stabilization==1){
                         /*Build K: */
@@ -3879 +3878 @@
                         vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
                         h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
+
                         K[0][0]=h/(2*vel)*vx*vx;  K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz;
                         K[1][0]=h/(2*vel)*vy*vx;  K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz;
                         K[2][0]=h/(2*vel)*vz*vx;  K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz;
+
                         D_scalar_stab=gauss->weight*Jdet;
                         if(reCast<bool,IssmDouble>(dt)) D_scalar_stab=D_scalar_stab*dt;
@@ -4070 +4071 @@
 
                 /*Advection: */
-
                 GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); 
                 GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); 
@@ -4960 +4960 @@
         IssmDouble  vx[NUMVERTICES],vy[NUMVERTICES],vz[NUMVERTICES];
         IssmDouble  geothermalflux[NUMVERTICES];
-        IssmDouble  dt;
-        IssmDouble  meltingrate_enthalpy;
+        IssmDouble  dt, yts;
+        IssmDouble  melting_overshoot,meltingrate_enthalpy;
+        IssmDouble  lambda,heating;
         Friction   *friction  = NULL;
 
@@ -4998 +4999 @@
                 checkpositivethickness=true;
 
+                _assert_(watercolumn[iv]>=0.);
+
                 /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
                 meltingrate_enthalpy=0.;
+                heating=0.;
                 if((watercolumn[iv]>0.) && (enthalpy[iv]<matpar->PureIceEnthalpy(pressure[iv]))){
                         /*ensure that no ice is at T<Tm(p), if water layer present*/
@@ -5032 +5036 @@
                         matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]);
                         // -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66
-                        meltingrate_enthalpy=(heatflux+basalfriction+geothermalflux[iv])/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent 
-                }
+                        heating=(heatflux+basalfriction+geothermalflux[iv]);
+                        meltingrate_enthalpy=heating/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent 
+                }
+
                 /*Update water column, basal meltingrate*/
-                //basalmeltingrate[iv]+=meltingrate_enthalpy;
-                basalmeltingrate[iv]=meltingrate_enthalpy;
                 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-                if(reCast<bool,IssmDouble>(dt))
-                        watercolumn[iv]+=dt*meltingrate_enthalpy; 
-                else
+                if(reCast<bool,IssmDouble>(dt)){
+                        if(watercolumn[iv]+meltingrate_enthalpy*dt<0.){                         
+                                melting_overshoot=watercolumn[iv]+meltingrate_enthalpy*dt;
+                                lambda=melting_overshoot/(meltingrate_enthalpy*dt); _assert_(lambda>0); _assert_(lambda<1);
+                                basalmeltingrate[iv]=(1.-lambda)*meltingrate_enthalpy;
+                                watercolumn[iv]=0.;
+                                yts=365*24*60*60;
+                                enthalpy[iv]+=dt/yts*lambda*heating;
+                        }
+                        else{
+                                basalmeltingrate[iv]=meltingrate_enthalpy;
+                                watercolumn[iv]+=dt*meltingrate_enthalpy; 
+                        }
+                }
+                else{
+                        basalmeltingrate[iv]=meltingrate_enthalpy;
                         watercolumn[iv]+=meltingrate_enthalpy;
+                }         
         }  
         /*feed updated variables back into model*/
@@ -5072 +5090 @@
         pentabase=this->GetBasalElement();
         
-        GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
-        GetInputListOnVertices(&pressure[0],PressureEnum);
+        this->GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
+        this->GetInputListOnVertices(&pressure[0],PressureEnum);
         pentabase->GetInputListOnVertices(&watercolumnbase[0], WatercolumnEnum);
 
@@ -5089 +5107 @@
                 watercolumnbase[iv%NUMVERTICES2D]+=dwc;
         }
+
         /*feed updated results back into model*/
         this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
@@ -5094 +5113 @@
         pentabase->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumnbase,P1Enum));
 
-        delete pentabase;
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

@@ -377 +377 @@
 /*FUNCTION Matpar::GetEnthalpyDiffusionParameter{{{*/
 IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){
-        if(enthalpy<PureIceEnthalpy(pressure)){
+        /*if (enthalpy<=PureIceEnthalpy(pressure))
                 return thermalconductivity/heatcapacity;
-        }
-        else{
-                return 1.045*1.e-4; // K0=1.045*1e-4 from Aschwanden 2012. TODO: fetch K0 from model
-        }
+        else
+                return thermalconductivity/heatcapacity*pow(10.,-5);
+        */
+
+        IssmDouble eps=0.05*heatcapacity;
+        IssmDouble hpmp=PureIceEnthalpy(pressure);
+        IssmDouble kappa_c=thermalconductivity/heatcapacity;
+        IssmDouble kappa_t=thermalconductivity/heatcapacity*pow(10.,-1);
+        return 1./(1.+exp(-(enthalpy-(hpmp))/eps))*(kappa_t-kappa_c) + kappa_c;
 }
 /*}}}*/