Index: /issm/trunk-jpl/src/c/analyses/enthalpy_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/enthalpy_core.cpp (revision 16358) +++ /issm/trunk-jpl/src/c/analyses/enthalpy_core.cpp (revision 16359) @@ -30,4 +30,6 @@ InputToResultx(femmodel,EnthalpyEnum); InputToResultx(femmodel,WaterfractionEnum); + InputToResultx(femmodel,WatercolumnEnum); + InputToResultx(femmodel,BasalforcingsMeltingRateEnum); } } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16358) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16359) @@ -4526,5 +4526,6 @@ IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0}; IssmDouble Jdet2d,dt; - IssmDouble rho_ice,heatcapacity,geothermalflux_value; + IssmDouble rho_ice,heatcapacity; + IssmDouble geothermalflux_value, heatflux_value; IssmDouble basalfriction,alpha2,vx,vy,vz; IssmDouble scalar,enthalpy,enthalpyup; @@ -4563,4 +4564,5 @@ gauss=new GaussPenta(0,1,2,2); gaussup=new GaussPenta(3,4,5,2); + for(int ig=gauss->begin();igend();ig++){ @@ -4575,5 +4577,5 @@ watercolumn_input->GetInputValue(&watercolumn,gauss); - if((watercolumn==0.) && (enthalpyPureIceEnthalpy(pressure))){ + if((watercolumn<=0.) && (enthalpyPureIceEnthalpy(pressure))){ /* the above check is equivalent to NOT ((watercolumn>0.) AND (enthalpyweight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice); + heatflux_value=(basalfriction+geothermalflux_value)/(rho_ice); + scalar=gauss->weight*Jdet2d*heatflux_value; if(reCast(dt)) scalar=dt*scalar; for(i=0;ivalues[i]+=scalar*basis[i]; @@ -5014,6 +5017,7 @@ bool isenthalpy,isdynamicbasalspc,setspc; int numindices, numindicesup; - IssmDouble h_pmp,pressure, pressureup; - IssmDouble enthalpy, enthalpyup; + IssmDouble pressure, pressureup; + IssmDouble h_pmp, enthalpy, enthalpyup; + IssmDouble watercolumn; int *indices = NULL, *indicesup = NULL; @@ -5035,4 +5039,5 @@ Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); //_assert_(watercolumn_input); /*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*/ @@ -5047,4 +5052,5 @@ enthalpy_input->GetInputValue(&enthalpy, gauss); pressure_input->GetInputValue(&pressure, gauss); + watercolumn_input->GetInputValue(&watercolumn,gauss); setspc = false; if (enthalpy>=matpar->PureIceEnthalpy(pressure)){ @@ -5054,5 +5060,5 @@ enthalpy_input->GetInputValue(&enthalpyup, gaussup); pressure_input->GetInputValue(&pressureup, gaussup); - setspc=(enthalpyupPureIceEnthalpy(pressureup))?true:false; + setspc=((enthalpyupPureIceEnthalpy(pressureup)) && (watercolumn>=0.))?true:false; } @@ -5096,5 +5102,6 @@ IssmDouble vx[NUMVERTICES],vy[NUMVERTICES],vz[NUMVERTICES]; IssmDouble geothermalflux[NUMVERTICES]; - IssmDouble dt,meltingrate_enthalpy; + IssmDouble dt; + IssmDouble meltingrate_enthalpy; Friction *friction = NULL; @@ -5108,5 +5115,5 @@ /*Fetch parameters and inputs */ latentheat=matpar->GetLatentHeat(); - rho_ice=this->matpar->GetRhoIce(); + rho_ice=matpar->GetRhoIce(); GetInputListOnVertices(&vx[0],VxEnum); GetInputListOnVertices(&vy[0],VyEnum); @@ -5149,8 +5156,8 @@ istemperatelayer=false; if(enthalpy[iv+3]>=matpar->PureIceEnthalpy(pressure[iv+3])) istemperatelayer=true; - if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.; + if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.; // TODO: add -k*nabla T_pmp else{ enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss); - kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure[iv]); _assert_(kappa>0.); + kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy[iv],enthalpy[iv+NUMVERTICES2D], pressure[iv],pressure[iv+NUMVERTICES2D]); _assert_(kappa>0.); for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i]; } @@ -5166,14 +5173,15 @@ matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]); - meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux[iv]))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent - } - + // -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 + } /*Update water column, basal meltingrate*/ - basalmeltingrate[iv]+=meltingrate_enthalpy; + //basalmeltingrate[iv]+=meltingrate_enthalpy; + basalmeltingrate[iv]=meltingrate_enthalpy; this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); if(reCast(dt)) - watercolumn[iv]+=dt*meltingrate_enthalpy; + watercolumn[iv]+=dt*meltingrate_enthalpy; else - watercolumn[iv]+=meltingrate_enthalpy; + watercolumn[iv]+=meltingrate_enthalpy; } /*feed updated variables back into model*/ @@ -5193,13 +5201,21 @@ bool isenthalpy; IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES]; + IssmDouble watercolumnbase[NUMVERTICES]; IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; IssmDouble latentheat, dt; + IssmDouble dw, dwc; + Penta *pentabase = NULL; /*Check wether enthalpy is activated*/ parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); if(!isenthalpy) return; - + + /*get basal element, needed for basal watercolumn*/ + pentabase=this->GetBasalElement(); + GetInputListOnVertices(&enthalpy[0],EnthalpyEnum); GetInputListOnVertices(&pressure[0],PressureEnum); + pentabase->GetInputListOnVertices(&watercolumnbase[0], WatercolumnEnum); + this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); latentheat=matpar->GetLatentHeat(); @@ -5207,13 +5223,18 @@ for(int iv=0;ivEnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]); - + dw=DrainageFunctionWaterfraction(waterfraction[iv], dt); /*drain water fraction & update enthalpy*/ - waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt); + waterfraction[iv]-=dw; matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]); + /*add drained water to watercolumn at base*/ + dwc=dw*this->IceVolume(); + watercolumnbase[iv%NUMVERTICES2D]+=dwc; } /*feed updated results back into model*/ this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum)); this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum)); - // this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum)); // temperature should not change during drainage... + pentabase->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumnbase,P1Enum)); + + delete pentabase; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16358) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16359) @@ -385,4 +385,25 @@ } /*}}}*/ +/*FUNCTION Matpar::GetEnthalpyDiffusionParameterVolume{{{*/ +IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy0, IssmDouble enthalpy1, IssmDouble pressure0, IssmDouble pressure1){ + /*returns kappa depending on distribution of enthalpy over edge of element + lambda is the barycentric coordinate that solves H0+(H1-H0)*lambda=H_pureice. + it represents the fraction of the ice column which is temperate/cold like H0. + if lambda<0 or lambda>1, then the whole ice column is cold or temperate, respectively. + */ + IssmDouble kappa, kappa0, kappa1; + if (enthalpy0!=enthalpy1){ + IssmDouble lambda=(PureIceEnthalpy(pressure0)-enthalpy0)/(enthalpy1-enthalpy0); + if ((lambda>=0.) && (lambda<=1.)){ + kappa0=GetEnthalpyDiffusionParameter(enthalpy0,pressure0); + kappa1=GetEnthalpyDiffusionParameter(enthalpy1,pressure1); + kappa=lambda*kappa0+(1.-lambda)*kappa1; + } + } + else + kappa=GetEnthalpyDiffusionParameter(enthalpy0, pressure0); + return kappa; +} +/*}}}*/ /*FUNCTION Matpar::EnthalpyToThermal {{{*/ void Matpar::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){ Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16358) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16359) @@ -126,4 +126,5 @@ IssmDouble PureIceEnthalpy(IssmDouble pressure); IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); + IssmDouble GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy0, IssmDouble enthalpy1, IssmDouble pressure0, IssmDouble pressure1); IssmDouble GetLithosphereShearModulus(); IssmDouble GetLithosphereDensity(); Index: /issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp (revision 16358) +++ /issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp (revision 16359) @@ -15,4 +15,6 @@ IssmDouble w0=0.01, w1=0.02, w2=0.03; IssmDouble Dret, D0=0, D1=0.005, D2=0.05; + IssmDouble yts=365*24*60*60; + dt/=yts; /*get drainage function value*/